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

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())