Additional Python Examples

Overview

Code

Python

#!/usr/bin/env python

from __future__ import print_function

import sys
import SimpleITK as sitk
import os

# verify that we have the correct number of arguments
if ( len( sys.argv ) != 5 ):
    sys.stderr.write( "Usage: prog inputFile outputFile replaceValue upperThreshold\n" )
    exit( 1 )
    
# copy the arguments in to variables
inputFileName = sys.argv[1]
outputFileName = sys.argv[2]
replaceValue = int( sys.argv[3] )
upperThreshold = float( sys.argv[4] )

# Read the file into an sitkImage
image = sitk.ReadImage( inputFileName )

# Threshold the value [0,2), results in values inside the range 1, 0 otherwise
boundary = sitk.BinaryThreshold( image, 0, upperThreshold, 1, 0 )

boundary = sitk.BinaryMorphologicalClosing( boundary, 1 )

# Remove any label pixel not connected to the boarder
boundary = sitk.BinaryGrindPeak( boundary )



boundary = sitk.Cast( boundary, image.GetPixelID() )

# Multiply, the input image by not the boarder.
# This will multiply the image by 0 or 1, where 0 is the
# boarder. Making the board 0
image = image * ~boundary

# add the replace value to the pixel on the board
image = image + ( boundary * replaceValue )


if ( not "SITK_NOSHOW" in os.environ ):
    sitk.Show( image, "Boarder Segmentation" )

#!/usr/bin/env python

from __future__ import print_function

import SimpleITK as sitk
import sys
import os

if len ( sys.argv ) < 2:
    print( "Usage: %s <input>" % ( sys.argv[0] ) )
    sys.exit ( 1 )


image = sitk.Cast( sitk.ReadImage( sys.argv[1] ), sitk.sitkFloat32 ) 

edges = sitk.CannyEdgeDetection( image, lowerThreshold=200, upperThreshold=400, variance=[4]*3 )

stats = sitk.StatisticsImageFilter()
stats.Execute( image )

if ( not "SITK_NOSHOW" in os.environ ):
    sitk.Show( sitk.Maximum( image*0.5, edges*stats.GetMaximum()*.5) )

'''=========================================================================
from __future__ import print_function

import SimpleITK as sitk
import sys
import os

#
# Check Command Line
#
if len( sys.argv ) < 7:
  print("Usage: ConnectedThresholdImageFilter inputImage outputImage lowerThreshold upperThreshold seedX seedY [seed2X seed2Y ... ]");
  sys.exit( 1 )


#
# Read the image
#
reader = sitk.ImageFileReader()
reader.SetFileName( sys.argv[1] )
image = reader.Execute();

#
# Blur using CurvatureFlowImageFilter
#
blurFilter = sitk.CurvatureFlowImageFilter()
blurFilter.SetNumberOfIterations( 5 )
blurFilter.SetTimeStep( 0.125 )
image = blurFilter.Execute( image )

#
# Set up ConnectedThresholdImageFilter for segmentation
#
segmentationFilter = sitk.ConnectedThresholdImageFilter()
segmentationFilter.SetLower( float(sys.argv[3]) )
segmentationFilter.SetUpper( float(sys.argv[4]) )
segmentationFilter.SetReplaceValue( 255 )

for i in range( 5, len(sys.argv)-1, 2 ):
  seed = [ int(sys.argv[i]), int(sys.argv[i+1]) ]
  segmentationFilter.AddSeed( seed )
  print( "Adding seed at: ", seed, " with intensity: ", image.GetPixel(*seed) )

# Run the segmentation filter
image = segmentationFilter.Execute( image )
image[seed] = 255

#
# Write out the result
#
writer = sitk.ImageFileWriter()
writer.SetFileName( sys.argv[2] )
writer.Execute( image )


if ( not "SITK_NOSHOW" in os.environ ):
  sitk.Show( image, "ConntectedThreshold" )

#=========================================================================
from __future__ import print_function

import SimpleITK as sitk

import sys, time

if len( sys.argv ) < 3:
    print( "Usage: python " + __file__ + " <input_image> <output_image>" )
    sys.exit ( 1 )

# Read the image, single 3D DICOM slice
image = sitk.ReadImage( sys.argv[1] )

# Modify the image (mean)
mean_image = sitk.BoxMean( image, [3,3,1] )

# Save the modified image:
# We do not provide a method that copies all keys blindly. Our intent is
# to force us to remember that we always need to modify the meta-data dicitonary
# keys when we save processed images in DICOM format.
# In case of the DICOM format, amongst other keys we need to set the Image Type (0008,0008),
# the Series Date (0008,0021), and Series Time (0008,0021). Obviously we need to set a
# different series number (0020,0011), series instance UID (0020,000E), etc. - we don't do
# that here.
# Please see the DICOM standard (http://dicom.nema.org/standard.html) for complete details on
# how to create valid DICOM images.

all_keys = image.GetMetaDataKeys()
for key in all_keys:
    mean_image.SetMetaData( key, image.GetMetaData( key ) )
mean_image.SetMetaData( "0008|0008", "DERIVED\SECONDARY" )
modification_time = time.strftime("%H%M%S")
modification_date = time.strftime("%Y%m%d")
mean_image.SetMetaData( "0008|0031", modification_time )
mean_image.SetMetaData( "0008|0021", modification_date )

sitk.WriteImage( mean_image, sys.argv[2] )

# Finally, read the image back and see that changes were made
# Note that the image type (0008|0008) can contain additional spaces. The requirement is that
# the string have an even length (ftp://dicom.nema.org/medical/DICOM/2013/output/chtml/part05/sect_6.2.html).
# Where spaces are added is not specified and thus may vary ("DERIVED\SECONDARY ", "DERIVED\ SECONDARY" are
# equivalent).
modified_image = sitk.ReadImage( sys.argv[2] )
if modified_image.GetMetaData( "0008|0008" ).replace(" ","") != "DERIVED\SECONDARY" or \
   modified_image.GetMetaData( "0008|0031" ) != modification_time or \
   modified_image.GetMetaData( "0008|0021" ) != modification_date:
    sys.exit(1)

sys.exit( 0 )

#!/usr/bin/env python
#=========================================================================


#
#  This example shows how to read a specific series from a Dicom directory that
#  may contain more than one series.  The script scans for all series.  If an
#  output name is given, it writes out the requested series.  If no specific
#  series name is given, the first series found is written.
#

from __future__ import print_function

import sys, getopt
import SimpleITK as sitk

target_series = ""
output_image = ""

def usage():
    print (  "\nUsage: %s [-s series_name] input_directory [output_image]\n" % (sys.argv[0]) )

# Parse command line options
try:
    opts, args = getopt.getopt(sys.argv[1:], "s:", [ "series" ] )
except getopt.GetoptError as err:
    usage()
    sys.exit(1)

for o, a in opts:
    if o in ("-s", "--series"):
        target_series = a
    else:
         assert False, "unhandled options"


# Get input/output names
if len(args) < 1:
    print( args )
    usage()
    sys.exit(1)

input_directory = args[0]
if len(args)>1:
    output_image = args[1]

# Find the Dicom series
reader = sitk.ImageSeriesReader()
written = False

series_found = reader.GetGDCMSeriesIDs(input_directory)

# Process each Dicom series
if len(series_found):

    for serie in series_found:

        print( "\nSeries:", serie )

        # Get the Dicom filename corresponding to the current series
        dicom_names = reader.GetGDCMSeriesFileNames(input_directory, serie)

        print( "\nFiles in series: ", dicom_names )

        if len(dicom_names):
            reader.SetFileNames(dicom_names)
            image = reader.Execute()
            print( "\nImage size: ", image.GetSize() )

            if (output_image != "") and not written:
                if (target_series == "" or target_series == serie):

                    print( "\nWriting", output_image )
                    sitk.WriteImage(image, output_image)
                    written = True
else:
    sys.exit(1)

print ()

#!/usr/bin/env python

from __future__ import print_function

import SimpleITK as sitk
import sys
import os

if len ( sys.argv ) != 4:
    print( "Usage: %s inputImage sliceNumber outputImage" % ( sys.argv[0] ) )
    sys.exit ( 1 )

zslice = int( sys.argv[2] )

inputImage = sitk.ReadImage( str(sys.argv[1]) )

size = list( inputImage.GetSize() )
size[2] = 0

index = [ 0, 0, zslice  ]

Extractor = sitk.ExtractImageFilter()
Extractor.SetSize( size )
Extractor.SetIndex( index )

sitk.WriteImage( Extractor.Execute( inputImage ), str(sys.argv[3]) )


if ( not "SITK_NOSHOW" in os.environ ):
    sitk.Show( Extractor.Execute( inputImage ) )

#!/usr/bin/env python

from __future__ import print_function
from __future__ import division

import SimpleITK as sitk
import sys
import os
import math

if len ( sys.argv ) < 4:
    print( "Usage: FFTConvolution <input> <kernel> <output>" )
    sys.exit ( 1 )


inputFileName = sys.argv[1]
kernelFileName = sys.argv[2]
outputFileName = sys.argv[3]

### Input Image ###
# read the input image
img = sitk.ReadImage( inputFileName )

# save the type of pixel the input is, so that we can cast the result
# back out to the same type  
pixelID = img.GetPixelID()

# pad the image 
img = sitk.MirrorPad( img, [128] *2, [128]*2 )
size = img.GetSize();

# perform the FFT
fftimg = sitk.ForwardFFT( sitk.Cast( img, sitk.sitkFloat32 ) )


### Kernel Image ###
# Read the kernel image file
kernel = sitk.ReadImage( kernelFileName )

# flip kernel about all axis
kernel = sitk.Flip( kernel, [True]*2 )


# normalize the kernel to sum to ~1
stats = sitk.StatisticsImageFilter();
stats.Execute( kernel )
kernel = sitk.Cast( kernel / stats.GetSum(), sitk.sitkFloat32 )

upadding = [0]*2
upadding[0] = int( math.floor( (size[0] - kernel.GetSize()[0])/2.0 ) )
upadding[1] = int( math.floor( (size[1] - kernel.GetSize()[1])/2.0 ) )

lpadding = [0]*2
lpadding[0] = int( math.ceil( (size[0] - kernel.GetSize()[0])/2.0 ) )
lpadding[1] = int( math.ceil( (size[1] - kernel.GetSize()[1])/2.0 ) )

# pad the kernel to prevent edge artifacts
kernel = sitk.ConstantPad( kernel, upadding, lpadding, 0.0 )

# perform FFT on kernel
fftkernel = sitk.ForwardFFT( sitk.FFTShift( kernel ) )

# meta-data must match for multiplication
fftkernel.SetSpacing( fftimg.GetSpacing() )
fftkernel.SetOrigin( fftimg.GetOrigin() )
fftkernel.SetDirection( fftimg.GetDirection() )


### Convolution ###
# Finally perform the convolution in Fourier space by multiplication
img =  sitk.InverseFFT( fftimg*fftkernel )

# remove the padding
img = sitk.Crop( img, [128]*2, [128]*2 )

### Writing ###
# write the output image the same type as the input
sitk.WriteImage( sitk.Cast( img, pixelID ), outputFileName )



if ( not "SITK_NOSHOW" in os.environ ):
    sitk.Show( sitk.ReadImage( inputFileName ), "original" )
    sitk.Show( sitk.ReadImage( kernelFileName ), "kernel" )
    sitk.Show( sitk.Cast( img, pixelID ), "FFT_Convolution" )

#=========================================================================
from __future__ import print_function

import SimpleITK as sitk
import os
import sys

if len ( sys.argv ) < 8:
    print( "Usage:", sys.argv[0], " <inputImage> <outputImage> <sigma> <InitialDistance> <PropagationScaling> <seedX> <seedY> <?seedZ>" )
    sys.exit ( 1 )

inputImage = sys.argv[1]
outputImage = sys.argv[2]
sigma = float(sys.argv[3])
initialDistance = float(sys.argv[4])
propagationScaling = float(sys.argv[5])
seed = [float(sys.argv[6]), float(sys.argv[7]) ]

if len( sys.argv ) > 8:
    seed.append( float(sys.argv[8]) )


reader = sitk.ImageFileReader()
reader.SetFileName ( inputImage )
image = reader.Execute()

gradientMagnitude = sitk.GradientMagnitudeRecursiveGaussianImageFilter()
gradientMagnitude.SetSigma( sigma )

featureImage = sitk.BoundedReciprocal( gradientMagnitude.Execute( image ) )

seedImage = sitk.Image( image.GetSize()[0], image.GetSize()[1], sitk.sitkUInt8 )
seedImage.SetSpacing( image.GetSpacing() )
seedImage.SetOrigin( image.GetOrigin() )
seedImage.SetDirection( image.GetDirection() )
seedImage[ seedImage.TransformPhysicalPointToIndex(seed) ] = 1

distance = sitk.SignedMaurerDistanceMapImageFilter()
distance.InsideIsPositiveOff()
distance.UseImageSpacingOn()

initialImage = sitk.BinaryThreshold( distance.Execute( seedImage ), -1000, 10 )
initialImage = sitk.Cast( initialImage, featureImage.GetPixelID() ) * -1 + 0.5


geodesicActiveContour = sitk.GeodesicActiveContourLevelSetImageFilter()
geodesicActiveContour.SetPropagationScaling( propagationScaling )
geodesicActiveContour.SetCurvatureScaling( .5 )
geodesicActiveContour.SetAdvectionScaling( 1.0 )
geodesicActiveContour.SetMaximumRMSError( 0.01 )
geodesicActiveContour.SetNumberOfIterations( 1000 )

levelset = geodesicActiveContour.Execute( initialImage, featureImage )

print( "RMS Change: ", geodesicActiveContour.GetRMSChange() )
print( "Elapsed Iterations: ", geodesicActiveContour.GetElapsedIterations() )

contour = sitk.BinaryContour( sitk.BinaryThreshold( levelset, -1000, 0 ) )

if ( not "SITK_NOSHOW" in os.environ ):
    sitk.Show( sitk.LabelOverlay( image, contour ), "Levelset Countour" )

#!/usr/bin/env python

from __future__ import print_function

import SimpleITK as sitk
import os

xImg = sitk.Image( 256, 256, sitk.sitkFloat32 )
yImg = sitk.Image( 256, 256, sitk.sitkFloat32 )

for y in range( 0, xImg.GetSize()[1] ):
    for x in range( 0, xImg.GetSize()[0] ):
        xImg.SetPixel( x, y, x )
        yImg[x, y] = y


sigma = 50

xImg = sitk.Subtract( xImg,  xImg.GetSize()[0] / 2 )
yImg = yImg - yImg.GetSize()[1] / 2

gaussianImg = sitk.Exp( -1 * (xImg**2 + yImg**2) / (2.0 * sigma**2) )


if ( not "SITK_NOSHOW" in os.environ ):
    sitk.Show( gaussianImg, "Gaussian Blob" )

'''=========================================================================
from __future__ import print_function

import SimpleITK as sitk
import sys
import os


#
# Check Command Line
#
if len( sys.argv ) < 7:
  print( "Usage: NeighborhoodConnectedImageFilter inputImage outputImage lowerThreshold upperThreshold seedX seedY [seed2X seed2Y ... ]")
  sys.exit( 1 )


#
# Read the image
#
reader = sitk.ImageFileReader()
reader.SetFileName( sys.argv[1] )
image = reader.Execute();


#
# Blur using CurvatureFlowImageFilter
#
blurFilter = sitk.CurvatureFlowImageFilter()
blurFilter.SetNumberOfIterations( 5 )
blurFilter.SetTimeStep( 0.125 )
image = blurFilter.Execute( image )

#
# Set up NeighborhoodConnectedImageFilter for segmentation
#
segmentationFilter = sitk.NeighborhoodConnectedImageFilter()
segmentationFilter.SetLower( float(sys.argv[3]) )
segmentationFilter.SetUpper( float(sys.argv[4]) )
segmentationFilter.SetReplaceValue( 255 )

radius = [2,2]
segmentationFilter.SetRadius( radius )

for i in range( 5, len(sys.argv)-1, 2 ):
  seed = [int(sys.argv[i]), int(sys.argv[i+1])]
  segmentationFilter.AddSeed( seed )
  print( "Adding seed at: ", seed, " with intensity: ", image.GetPixel(*seed) )

# Run the segmentation filter
image = segmentationFilter.Execute( image )

#
# Write out the result
#
writer = sitk.ImageFileWriter()
writer.SetFileName( sys.argv[2] )
writer.Execute( image )



if ( not "SITK_NOSHOW" in os.environ ):
  sitk.Show( image, "NeighborhoodConnectedThreshold" )

#!/usr/bin/env python

from __future__ import print_function

import SimpleITK as sitk
import sys
import os

if len ( sys.argv ) != 2:
    print( "Usage: %s inputImage" % ( sys.argv[0] ) )
    sys.exit ( 1 )

inputImage = sitk.ReadImage( sys.argv[1] )


if ( not "SITK_NOSHOW" in os.environ ):
    sitk.Show( inputImage )

#!/usr/bin/env python

from __future__ import print_function

import SimpleITK as sitk
import sys
import os

if len ( sys.argv ) < 4:
    print( "Usage: %s <input> <sigma> <output>" % ( sys.argv[0] ) )
    sys.exit ( 1 )


image = sitk.ReadImage( sys.argv[1] )

pixelID = image.GetPixelID()

image  = sitk.SmoothingRecursiveGaussian( image,  float( sys.argv[2] ) )

sitk.WriteImage( sitk.Cast( image, pixelID ), sys.argv[3] )


if ( not "SITK_NOSHOW" in os.environ ):
    sitk.Show( sitk.Cast( image, pixelID ), "Simple Gaussian Procedural" )