SimpleITK image read and write

.gitignore

.env
*.nii.gz
tmp

i.png

io_utils.py

# Helper methods for image reading and writing with SimpleITK.
# Source: https://gist.github.com/pangyuteng/7f54dbfcd67fb9d43a85f8c6818fca7b
# Author: Pangyu Teng
# Licence: https://en.wikipedia.org/wiki/MIT_License

import os
import SimpleITK as sitk

def imread(fpath):
    if fpath.endswith('.list'):
        with open(fpath,'r') as f:
            dicom_names = [x for x in f.read().split('\n') if len(x) > 0]
        if not os.path.exists(dicom_names[0]):
            dicom_names = [os.path.join(os.path.dirname(fpath),x) for x in dicom_names]
        reader = sitk.ImageSeriesReader()
        reader.SetFileNames(dicom_names)
    else:
        reader= sitk.ImageFileReader()
        reader.SetFileName(fpath)
    img = reader.Execute()
    arr = sitk.GetArrayFromImage(img)    
    spacing = img.GetSpacing()
    origin = img.GetOrigin()
    direction = img.GetDirection()    
    return arr,spacing,origin,direction

def imwrite(fpath,arr,spacing,origin,direction,use_compression=True):
    img = sitk.GetImageFromArray(arr)
    img.SetSpacing(spacing)
    img.SetOrigin(origin)
    img.SetDirection(direction)
    writer = sitk.ImageFileWriter()    
    writer.SetFileName(fpath)
    writer.SetUseCompression(use_compression)
    writer.Execute(img)

orient.py

import numpy as np
import SimpleITK as sitk
import os
import sys
import imageio.v2 as imageio

def get_feet_first_obj():
    s_png = imageio.imread("s.png")
    i_png = imageio.imread("i.png")
    arr = np.zeros((256,256,256))
    for x in range(128):
        arr[x,:,:] = i_png[:,:,0]
    for x in range(128):
        arr[-1*(x+1),:,:] = s_png[:,:,0]
    arr = arr.astype(np.int16)
    src_obj = sitk.GetImageFromArray(arr)
    src_direction = (1,0,0,0,1,0,0,0,1)
    src_obj.SetDirection(src_direction)
    src_obj.SetSpacing((1,1,1))
    src_obj.SetOrigin((0,0,0))
    return src_obj

def get_head_first_obj():
    s_png = imageio.imread("s.png")
    i_png = imageio.imread("i.png")
    arr = np.zeros((256,256,256))
    for x in range(128):
        arr[x,:,:] = s_png[:,:,0]
    for x in range(128):
        arr[-1*(x+1),:,:] = i_png[:,:,0]
    arr = arr.astype(np.int16)
    src_obj = sitk.GetImageFromArray(arr)
    src_direction = (1,0,0,0,1,0,0,0,-1)
    src_obj.SetDirection(src_direction)
    src_obj.SetSpacing((1,1,1))
    src_obj.SetOrigin((0,0,0))
    return src_obj

def get_head_first_prone_obj():
    s_png = imageio.imread("s.png")
    s_png = np.flip(s_png,axis=(0,1))
    i_png = imageio.imread("i.png")
    i_png = np.flip(i_png,axis=(0,1))
    arr = np.zeros((256,256,256))
    for x in range(128):
        arr[x,:,:] = s_png[:,:,0]
    for x in range(128):
        arr[-1*(x+1),:,:] = i_png[:,:,0]
    arr = arr.astype(np.int16)
    src_obj = sitk.GetImageFromArray(arr)
    src_direction = (-1,0,0,0,-1,0,0,0,-1)
    src_obj.SetDirection(src_direction)
    src_obj.SetSpacing((1,1,1))
    src_obj.SetOrigin((0,0,0))
    return src_obj


# iterp=sitk.sitkLinear,sitk.sitkNearestNeighbor
def orient_image(src_obj,tgt_obj,interp=sitk.sitkLinear,out_value=0):
    # special
    if isinstance(tgt_obj,str) and tgt_obj == "head_first":
            head_first_direction = (1,0,0,0,1,0,0,0,-1)

            origin_index = [None,None,None]
            if src_obj.GetDirection()[0]==head_first_direction[0]:
                origin_index[0]=0
            else:
                origin_index[0]=src_obj.GetSize()[0]-1

            if src_obj.GetDirection()[4]==head_first_direction[4]:
                origin_index[1]=0
            else:
                origin_index[1]=src_obj.GetSize()[1]-1

            if src_obj.GetDirection()[8]==head_first_direction[8]:
                origin_index[2]=0
            else:
                origin_index[2]=src_obj.GetSize()[2]-1

            origin_physical = src_obj.TransformIndexToPhysicalPoint(origin_index)
            zeros_arr = np.zeros_like(sitk.GetArrayFromImage(src_obj)).astype(np.int16)
            tgt_obj = sitk.GetImageFromArray(zeros_arr)
            tgt_obj.SetSpacing(src_obj.GetSpacing())
            tgt_obj.SetOrigin(origin_physical)
            tgt_obj.SetDirection(head_first_direction)
    else:
        raise NotImplementedError()

    transform = sitk.AffineTransform(3) # assume 3d

    return sitk.Resample(src_obj, tgt_obj,
        transform, interp,
        out_value, src_obj.GetPixelID())

feet_first_obj = get_feet_first_obj()
sitk.WriteImage(feet_first_obj,'tmp/feet_first.nii.gz')

head_first_obj = get_head_first_obj()
sitk.WriteImage(head_first_obj,'tmp/head_first.nii.gz')

head_first_prone_obj = get_head_first_prone_obj()
sitk.WriteImage(head_first_prone_obj,'tmp/head_first_prone.nii.gz')

print(feet_first_obj.GetDirection())
print(head_first_obj.GetDirection())
print(head_first_prone_obj.GetDirection())

feet_first_arr = sitk.GetArrayFromImage(head_first_obj).copy()
arr = sitk.GetArrayFromImage(head_first_obj)
feet = arr[0,:,:]
head = arr[-1,:,:]
imageio.imsave('tmp/feet-first-head.png',head.astype(np.uint8))
imageio.imsave('tmp/feet-first-feet.png',feet.astype(np.uint8))

head_first_arr = sitk.GetArrayFromImage(head_first_obj).copy()
arr = sitk.GetArrayFromImage(head_first_obj)
head = arr[0,:,:]
feet = arr[-1,:,:]
imageio.imsave('tmp/head-first-head.png',head.astype(np.uint8))
imageio.imsave('tmp/head-first-feet.png',feet.astype(np.uint8))

head_first_prone_arr = sitk.GetArrayFromImage(head_first_prone_obj).copy()
arr = sitk.GetArrayFromImage(head_first_prone_obj)
head = arr[0,:,:]
feet = arr[-1,:,:]
imageio.imsave('tmp/head-first-prone-head.png',head.astype(np.uint8))
imageio.imsave('tmp/head-first-prone-feet.png',feet.astype(np.uint8))

ro_head_first_obj = orient_image(feet_first_obj,"head_first",interp=sitk.sitkNearestNeighbor)
diff_arr = np.abs(sitk.GetArrayFromImage(feet_first_obj)-sitk.GetArrayFromImage(ro_head_first_obj))
print(np.min(diff_arr),np.max(diff_arr))
print(np.mean(diff_arr),'why no diff???')
diff_arr = feet_first_arr-sitk.GetArrayFromImage(ro_head_first_obj)
print(np.mean(diff_arr),'why no diff???')
sitk.WriteImage(ro_head_first_obj,'tmp/ro-1.nii.gz')
arr = sitk.GetArrayFromImage(ro_head_first_obj)
head = arr[0,:,:]
feet = arr[-1,:,:]
imageio.imsave('tmp/ro-1-head.png',head.astype(np.uint8))
imageio.imsave('tmp/ro-1-feet.png',feet.astype(np.uint8))

ro_head_first_obj = orient_image(head_first_obj,"head_first",interp=sitk.sitkNearestNeighbor)
diff_arr = np.abs(sitk.GetArrayFromImage(head_first_obj)-sitk.GetArrayFromImage(ro_head_first_obj))
print(np.min(diff_arr),np.max(diff_arr))
print(np.mean(diff_arr),'should be same ????')
diff_arr = head_first_arr-sitk.GetArrayFromImage(ro_head_first_obj)
print(np.mean(diff_arr),'should be same ????')
sitk.WriteImage(ro_head_first_obj,'tmp/ro-2.nii.gz')
arr = sitk.GetArrayFromImage(ro_head_first_obj)
head = arr[0,:,:]
feet = arr[-1,:,:]
imageio.imsave('tmp/ro-2-head.png',head.astype(np.uint8))
imageio.imsave('tmp/ro-2-feet.png',feet.astype(np.uint8))

ro_head_first_obj = orient_image(head_first_prone_obj,"head_first",interp=sitk.sitkNearestNeighbor)
diff_arr = np.abs(sitk.GetArrayFromImage(head_first_prone_obj)-sitk.GetArrayFromImage(ro_head_first_obj))
print(np.min(diff_arr),np.max(diff_arr))
print(np.mean(diff_arr),'why no diff???')
diff_arr = head_first_prone_arr-sitk.GetArrayFromImage(ro_head_first_obj)
print(np.mean(diff_arr),'why no diff??? need to cast to float!!!')
sitk.WriteImage(ro_head_first_obj,'tmp/ro-3.nii.gz')
arr = sitk.GetArrayFromImage(ro_head_first_obj)
head = arr[0,:,:]
feet = arr[-1,:,:]
imageio.imsave('tmp/ro-3-head.png',head.astype(np.uint8))
imageio.imsave('tmp/ro-3-feet.png',feet.astype(np.uint8))

"""

docker run -it -u $(id -u):$(id -g) \
    -w $PWD -v $PWD:$PWD pangyuteng/ml:latest bash

"""

s.png

basics:
https://simpleitk.readthedocs.io/en/v1.2.4/Documentation/docs/source/fundamentalConcepts.html

resampling method a (better than b IMO)
https://gist.github.com/mrajchl/ccbd5ed12eb68e0c1afc5da116af614a

def resample_img(itk_image, out_spacing=[2.0, 2.0, 2.0], is_label=False):
    
    # Resample images to 2mm spacing with SimpleITK
    original_spacing = itk_image.GetSpacing()
    original_size = itk_image.GetSize()

    out_size = [
        int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0]))),
        int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1]))),
        int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])))]

    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(out_spacing)
    resample.SetSize(out_size)
    resample.SetOutputDirection(itk_image.GetDirection())
    resample.SetOutputOrigin(itk_image.GetOrigin())
    resample.SetTransform(sitk.Transform())
    resample.SetDefaultPixelValue(itk_image.GetPixelIDValue())

    if is_label:
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        resample.SetInterpolator(sitk.sitkBSpline)

    return resample.Execute(itk_image)

# Assume to have some sitk image (itk_image) and label (itk_label)
resampled_sitk_img = resample_img(itk_image, out_spacing=[2.0, 2.0, 2.0], is_label=False)
resampled_sitk_lbl = resample_img(itk_label, out_spacing=[2.0, 2.0, 2.0], is_label=True)

--

**resampling method b:** 

# Resample using the three interpolators and combine the results so that we view the
# images one after the other, highlights the boundary differences
resampled_segmentation_nn = sitk.Resample(
    segmentation,
    new_size,
    sitk.Transform(),
    sitk.sitkNearestNeighbor,
    segmentation.GetOrigin(),
    new_spacing,
    segmentation.GetDirection(),
    0,
    segmentation.GetPixelID(),
)

# if you need GetArrayFromImage to respect  SetDirection:
sitk.Resample(mask_obj, img_obj, sitk.AffineTransform(3), sitk.sitkNearestNeighbor, 0, mask_obj.GetPixelID())

pretty much the same amount of work as method a....

import sys
import numpy as np
import SimpleITK as sitk


def resample(src_obj,out_size,method=sitk.sitkNearestNeighbor):
    src_size = np.array(src_obj.GetSize())
    src_spacing = np.array(src_obj.GetSpacing())
    tgt_size = np.array(out_size)
    tgt_spacing = src_size*src_spacing/tgt_size

    ref_obj = sitk.Image(out_size, sitk.sitkInt16)
    ref_obj.SetDirection(src_obj.GetDirection())
    ref_obj.SetOrigin(src_obj.GetOrigin())
    ref_obj.SetSpacing(tgt_spacing)

    tgt_obj = sitk.Resample(
        src_obj, ref_obj, sitk.Transform(),
        method, 0, src_obj.GetPixelID()
    )
    return tgt_obj

if __name__ == "__main__":
    img_file = sys.argv[1]
    out_file = sys.argv[2]
    out_size = (128,128,128)
    tgt_obj = sitk.ReadImage(img_file)
    out_obj = resample(tgt_obj,out_size)
    sitk.WriteImage(out_obj,out_file)

read region from nifti

https://stackoverflow.com/questions/63189077/simpleitk-reading-a-slice-of-an-image
https://simpleitk.readthedocs.io/en/master/link_AdvancedImageReading_docs.html

dicom file creation with pydicom

pydicom/pydicom#953
https://stackoverflow.com/questions/14350675/create-pydicom-file-from-numpy-array

rtstruct

https://github.com/qurit/rt-utils

ct_img = sitk.ReadImage(ct_scan_path)
nda = sitk.GetArrayFromImage(ct_img)
img_preprocessed = sitk.GetImageFromArray(nda_preprocessed)
img_preprocessed.CopyInformation(ct_img)
...
sitk.WriteImage(image, outputImageFileName)

https://itk.org/pipermail/community/2017-November/013783.html

cropping:

# Get the bounding box of the anatomy
label_shape_filter = sitk.LabelShapeStatisticsImageFilter()    
label_shape_filter.Execute(label_obj)
...
bounding_box = label_shape_filter.GetBoundingBox(roi_idx)
# The bounding box's first "dim" entries are the starting index and last "dim" entries the size
cropped_size = bounding_box[int(len(bounding_box)/2):]
cropped_start = bounding_box[0:int(len(bounding_box)/2)]
cropped_img_obj = sitk.RegionOfInterest(img_obj,cropped_size,cropped_start)
...
# ref. https://simpleitk.org/SPIE2019_COURSE/03_data_augmentation.html

pasting:

blank = np.zeros_like(sitk.GetArrayFromImage(img_obj)).astype(np.uint8)
mask_obj = sitk.GetImageFromArray(blank)
mask_obj.CopyInformation(img_obj)

lobes_obj = sitk.Paste(
    destinationImage=mask_obj,
    sourceImage=cropped_obj,
    sourceSize=cropped_size,
    sourceIndex=[0,0,0],
    destinationIndex=cropped_start)

label, regionprops
https://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/Python_html/35_Segmentation_Shape_Analysis.html

# skimage

from skimage.measure import label, regionprops

label_image = label(arr>0)
region_list = regionprops(label_image)
region_list = sorted(region_list,key=lambda x: x.area,reverse=True)

# sitk

shape_stats = sitk.LabelShapeStatisticsImageFilter()
shape_stats.ComputeOrientedBoundingBoxOn()
shape_stats.Execute(non_border_seg)
...
shape_stats.GetPhysicalSize(i),
shape_stats.GetElongation(i),
shape_stats.GetFlatness(i),
shape_stats.GetOrientedBoundingBoxSize(i),
intensity_stats.GetMean(i),
intensity_stats.GetStandardDeviation(i),
intensity_stats.GetSkewness(i),

more copy pastas

head first

HEAD_FIRST_ORIENTATION = (1,0,0,0,1,0,0,0,-1)
...
ref_space_obj = sitk.Image(source_obj.GetSize(), source_obj.GetPixelID())
ref_space_obj.SetOrigin(origin_physical)
ref_space_obj.SetDirection(HEAD_FIRST_ORIENTATION)
ref_space_obj.SetSpacing(source_obj.GetSpacing())

transform = sitk.AffineTransform(3)

return sitk.Resample(source_obj, ref_space_obj, 
  transform, sitk.sitkNearestNeighbor, 
  out_value, source_obj.GetPixelID())

rotation

https://stackoverflow.com/a/69098613/868736

def get_center(img_obj):
    width, height, depth = img_obj.GetSize()
    return img_obj.TransformIndexToPhysicalPoint((int(np.ceil(width/2)),
                                              int(np.ceil(height/2)),
                                              int(np.ceil(depth/2))))

origin_physical = source_obj.TransformIndexToPhysicalPoint(origin_index)
# transform = sitk.Transform()
...
rotation_center = get_center(source_obj)
theta_x = np.deg2rad(180) if source_orient_vec[0] != HEAD_FIRST_ORIENTATION[0] else 0
theta_y = np.deg2rad(180) if source_orient_vec[4] != HEAD_FIRST_ORIENTATION[4] else 0
theta_z = np.deg2rad(180) if source_orient_vec[8] != HEAD_FIRST_ORIENTATION[8] else 0
print(theta_x,theta_y,theta_z)
translation = (0,0,0)
transform = sitk.Euler3DTransform(
    rotation_center, theta_x, theta_y, theta_z, translation
)
return sitk.Resample(source_obj, ref_space_obj, 
    transform, sitk.sitkNearestNeighbor, 
    out_value, source_obj.GetPixelID())

flip

https://stackoverflow.com/questions/63873094/simple-itk-is-not-behaving-as-expected-while-flip-using-affine-transform

    center = get_center(source_obj)
    iden_x = -1 if source_orient_vec[0] != HEAD_FIRST_ORIENTATION[0] else 1
    iden_y = 1 if source_orient_vec[4] != HEAD_FIRST_ORIENTATION[4] else 1
    iden_z = 1 if source_orient_vec[8] != HEAD_FIRST_ORIENTATION[8] else 1
    transform = sitk.AffineTransform(3)
    transform.SetMatrix([iden_x,0,0,0,iden_y,0,0,0,iden_z])
    transform.SetCenter(center)
    print(transform.GetMatrix())

    return sitk.Resample(source_obj, ref_space_obj, 
        transform, sitk.sitkNearestNeighbor, 
        out_value, source_obj.GetPixelID())