Image Grid Manipulation

Overview

There are numerous SimpleITK filters that have similar functions, but very important differences. The filters that will be compared are:

• JoinSeriesImageFilter() - Joins multiple N-D images into an (N+1)-D image

• ComposeImageFilter() - Combines several scalar images into a multi-component vector image

• VectorIndexSelectionCastImageFilter() - Extracts the selected index of the input pixel type vector (the input image pixel type must be a vector and the output a scalar). Additionally, this filter can cast the output pixel type (SetOutputPixelType method).

• ExtractImageFilter() - Crops an image to the selected region bounds using vectors; Collapses dimensions unless dimension is two

• CropImageFilter() - Similar to ExtractImageFilter(), but crops an image by an itk::Size at upper and lower bounds

• Image Slicing Operator - Uses slicing (img[i:j, k:l]) to extract a subregion of an image

All of these operations will maintain the physical location of the pixels, instead modifying the image’s metadata.

Comparisons

Composition Filters

While JoinSeriesImageFilter() merges multiple images of the same pixel type in N dimensions into an image in N+1 dimensions, ComposeImageFilter() will combine scalar images into a vector image of the same dimension. The former is useful for connecting a series of contiguous images while the latter is more useful for merging multiple channels of the same object into one image (such as RGB).

Extraction Filters

VectorIndexSelectionCastImageFilter() will isolate a single channel in a vector image and return a scalar image. On the other hand, ExtractImageFilter() and CropImageFilter() will extract and return a subregion of an image, using an ExtractionRegion size and index and itk::Size’s respectively. However, note that only the ExtractImageFilter() collapses dimensions. The image slicing operator can also serve the same purpose.

Code

PythonR

#!/usr/bin/env python
import sys
import SimpleITK as sitk

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

# Two vector images of same pixel type and dimension expected
image_1 = sitk.ReadImage(sys.argv[1])
image_2 = sitk.ReadImage(sys.argv[2])

# Join two N-D Vector images to form an (N+1)-D image
join = sitk.JoinSeriesImageFilter()
joined_image = join.Execute(image_1, image_2)

# Extract first three channels of joined image (assuming RGB)
select = sitk.VectorIndexSelectionCastImageFilter()
channel1_image = select.Execute(joined_image, 0, sitk.sitkUInt8)
channel2_image = select.Execute(joined_image, 1, sitk.sitkUInt8)
channel3_image = select.Execute(joined_image, 2, sitk.sitkUInt8)

# Recompose image (should be same as joined_image)
compose = sitk.ComposeImageFilter()
composed_image = compose.Execute(channel1_image, channel2_image, channel3_image)

# Select same subregion using image slicing operator
sliced_image = composed_image[100:400, 100:400, 0]

# Select same subregion using ExtractImageFilter
extract = sitk.ExtractImageFilter()
extract.SetSize([300, 300, 0])
extract.SetIndex([100, 100, 0])
extracted_image = extract.Execute(composed_image)

# Select same sub-region using CropImageFilter (NOTE: CropImageFilter cannot
# reduce dimensions unlike ExtractImageFilter, so cropped_image is a three
# dimensional image with depth of 1)
crop = sitk.CropImageFilter()
crop.SetLowerBoundaryCropSize([100, 100, 0])
crop.SetUpperBoundaryCropSize(
    [composed_image.GetWidth() - 400, composed_image.GetHeight() - 400, 1]
)
cropped_image = crop.Execute(composed_image)