Simple ITK (SITK) + ITK hacks (python)

conversion.py

import itk
import SimpleITK as sitk

def convert_sitk_to_itk(sitk_img):
    
    img_array = sitk.GetArrayFromImage(sitk_img)
    itk_img = itk.GetImageFromArray(img_array)
    itk_img.SetOrigin(tuple(sitk_img.GetOrigin()))
    itk_img.SetSpacing(tuple(sitk_img.GetSpacing()))
    itk_img.SetDirection(itk.GetMatrixFromArray(np.reshape(np.array(sitk_img.GetDirection()), [sitk_img.GetDimension()]*2)))

    return itk_img
 
def convert_itk_to_sitk(image_itk, dtype):
    # cast = itk.UC, itk.F ??
        
    img_array = itk.GetArrayFromImage(image_itk)
    if dtype in ['short', 'int16']:    
        img_array = np.array(img_array, dtype=np.int16)
    elif dtype in ['unsigned int', 'uint8']:
        img_array = np.array(img_array, dtype=np.uint8)
        
    image_sitk = sitk.GetImageFromArray(img_array, isVector=image_itk.GetNumberOfComponentsPerPixel()>1)
    image_sitk.SetOrigin(tuple(image_itk.GetOrigin()))
    image_sitk.SetSpacing(tuple(image_itk.GetSpacing()))
    image_sitk.SetDirection(itk.GetArrayFromMatrix(image_itk.GetDirection()).flatten())
    return image_sitk

itk.py

import itk

# 1 - Basics
img = itk.imread('') # read some 3D volume
## itk.imread('', itk.F) float
## itk.imread('', itk.UC) # unsigned-char (e.g. for segmentation masks)
print (' - dimension: ', img.GetImageDimension())
print (' - size     : ', img.GetLargestPossibleRegion().GetSize())
print (' - spacing  : ', img.GetSpacing())
print (' - origin   : ', img.GetOrigin())
print (' - direction: ', itk.GetArrayFromMatrix(img.GetDirection()))
print (' - pixel_val: ', itk.GetPixel((h,w,d))) # (h,w,d) - voxel index
img_array = itk.GetArrayFromImage(img)
itk.imwrite(img, 'filename.mha')

sitk.py

import SimpleITK as sitk

img = sitk.ReadImage('') # some 3D volume
print (' - dimension: ', img.GetDimension())
print (' - size     : ', img.GetSize())
print (' - spacing  : ', img.GetSpacing())
print (' - origin   : ', img.GetOrigin())
print (' - direction: ', img.GetDirection())
print (' - pixel_val: ', itk.GetPixel((h,w,d))) # (h,w,d) - voxel index
print (' - type     : ', img.GetPixelIDTypeAsString())
img_array = sitk.GetArrayFromImage(img)
sitk.WriteImage(img, 'filename.mha', useCompression=True)


"""
- References
-- [Migration Guide](https://simpleitk.readthedocs.io/en/master/migrationGuide2.0.html)
"""

vol_to_dicom_for_ct.py

def vol_to_dicom_for_ct(path_img_ct, patient_name, patient_id, path_dicom):

    """
    Converts a .nrrd/.mha/.nifti file into its .dcm files

    Params
    ------
    path_img_ct: str, the path of the .nrrd/.mha/.nifti file
    patient_name: str
    patient_id: str
    path_dicom: str, the final output directory
    
    Note: Verify the output with dciodvfy
        - Ref 1: https://www.dclunie.com/dicom3tools/workinprogress/index.html
        - Ref 2: https://manpages.debian.org/unstable/dicom3tools/dciodvfy.1.en.html
        - Ref 3: # Motivation: https://stackoverflow.com/questions/14350675/create-pydicom-file-from-numpy-array
    """

    try:
        
        import sys
        import copy
        import random
        import shutil
        import subprocess
        import numpy as np

        if Path(path_img_ct).exists():
            
            try:
                import pydicom
                import pydicom._storage_sopclass_uids
            except:
                subprocess.check_call([sys.executable, '-m', 'pip', 'install', '--user', 'pydicom'])
                import pydicom  

            try:
                import SimpleITK as sitk
            except:
                subprocess.check_call([sys.executable, '-m', 'pip', 'install', '--user', 'SimpleITK']) # 2.1.1
                import SimpleITK as sitk
            
            try:
                import matplotlib.pyplot as plt
            except:
                subprocess.check_call([sys.executable, '-m', 'pip', 'install', '--user', 'matplotlib']) # 2.1.1
                import matplotlib.pyplot as plt

            # Step 0 - Create save directory
            if Path(path_dicom).exists():
                shutil.rmtree(path_dicom)
            Path(path_dicom).mkdir(exist_ok=True, parents=True)
            
            # Step 1 - Get volume params
            img_ct      = sitk.ReadImage(str(path_img_ct))
            img_spacing = tuple(img_ct.GetSpacing())
            img_origin  = tuple(img_ct.GetOrigin()) # --> dicom.ImagePositionPatient
            img_array   = sitk.GetArrayFromImage(img_ct).astype(np.int16) # [D,H,W]

            # Step 2 - Create dicom dataset
            ds                     = pydicom.dataset.Dataset()
            ds.FrameOfReferenceUID = pydicom.uid.generate_uid() # this will stay the same for all .dcm files of a volume
            
            # Step 2.1 - Modality details
            ds.SOPClassUID         = pydicom._storage_sopclass_uids.CTImageStorage
            ds.Modality            = 'CT'
            ds.ImageType           = ['ORIGINAL', 'PRIMARY', 'AXIAL']

            # Step 2.2 - Image Details
            ds.PixelSpacing               = [float(img_spacing[0]), float(img_spacing[1])]
            ds.SliceThickness             = str(img_spacing[-1])
            ds.Rows                       = img_array.shape[1]
            ds.Columns                    = img_array.shape[2]
            
            ds.PatientPosition            = 'HFS'
            ds.ImageOrientationPatient    = [1, 0, 0, 0, 1, 0]
            ds.PositionReferenceIndicator = 'SN'
            
            ds.SamplesPerPixel            = 1
            ds.PhotometricInterpretation  = 'MONOCHROME2'
            ds.BitsAllocated              = 16
            ds.BitsStored                 = 16
            ds.HighBit                    = 15
            ds.PixelRepresentation        = 1
            
            ds.RescaleIntercept           = "0.0"
            ds.RescaleSlope               = "1.0"
            ds.RescaleType                = 'HU'
            
            # Step 3.1 - Metadata
            fileMeta = pydicom.Dataset()
            fileMeta.MediaStorageSOPClassUID = pydicom._storage_sopclass_uids.CTImageStorage
            fileMeta.MediaStorageSOPInstanceUID = pydicom.uid.generate_uid() # this will change for each .dcm file of a volume
            fileMeta.TransferSyntaxUID = pydicom.uid.ExplicitVRLittleEndian
            ds.file_meta = fileMeta

            # Step 3.2 - Include study details
            ds.StudyInstanceUID    = pydicom.uid.generate_uid()
            ds.StudyDescription    = ''
            ds.StudyDate           = '19000101'                   # needed to create DICOMDIR
            ds.StudyID             = str(random.randint(0,1000))  # needed to create DICOMDIR

            # Step 3.3 - Include series details
            ds.SeriesInstanceUID   = pydicom.uid.generate_uid()
            ds.SeriesDescription   = ''
            ds.SeriesNumber        = str(random.randint(0,1000)) # needed to create DICOMDIR

            # Step 3.4 - Include patient details
            ds.PatientName      = patient_name
            ds.PatientID        = patient_id

            # Step 3.5 - Manufacturer details
            ds.Manufacturer           = 'MICCAI2015'
            ds.ReferringPhysicianName = 'Mody'                   # needed for identification in RayStation
            ds.ManufacturerModelName  = 'test_offsite'

            # Step 4 - Make slices
            for slice_id in range(img_array.shape[0]):
                
                # Step 4.1 - Slice identifier
                random_uuid = pydicom.uid.generate_uid()
                ds.file_meta.MediaStorageSOPInstanceUID = random_uuid
                ds.SOPInstanceUID = random_uuid
                ds.InstanceNumber = str(slice_id+1)
                
                vol_origin_tmp          = list(copy.deepcopy(img_origin))
                vol_origin_tmp[-1]      += img_spacing[-1]*slice_id
                ds.ImagePositionPatient = vol_origin_tmp

                # Step 4.2 - Slice data
                img_slice               = img_array[slice_id,:,:]
                # plt.imshow(img_slice); plt.savefig(str(Path(path_dicom, '{}.png'.format(slice_id)))); plt.close()
                ds.PixelData            = img_slice.tobytes()

                save_path               = Path(path_dicom).joinpath(str(ds.file_meta.MediaStorageSOPInstanceUID) + '.dcm')
                ds.save_as(str(save_path), write_like_original=False)

            return ds.StudyInstanceUID, ds.SeriesInstanceUID

        else:
            print (' - [ERROR][vol_to_dicom_for_ct()] Error in path: path_img_ct: ', path_img_ct)
            return None, None

    except:
        traceback.print_exc()

weird_itk_sitk_write.py

"""
Shows that when changing shape of a raw numpy array, sitk handles it, but itk does not. Weird :p
- Ref: https://discourse.itk.org/t/importing-image-from-array-and-axis-reorder/1192/10
"""
import itk
import SimpleITK as sitk
import numpy as np

array = np.random.random((240,240,80))
spacing=(0.8,0.8,2.5)
origin=(1,1,1)

img_itk = itk.GetImageFromArray(array)
img_itk.SetOrigin(tuple(origin))
img_itk.SetSpacing(tuple(spacing))
itk.imwrite(img_itk, 'itk_sample1.mha')

img_sitk = sitk.GetImageFromArray(array)
img_sitk.SetOrigin(tuple(origin))
img_sitk.SetSpacing(tuple(spacing))
sitk.WriteImage(img_sitk, 'sitk_sample1.mha')

array_new = np.moveaxis(array, [0,1,2], [2,1,0]).copy() # [NOTE] Need to use .copy() here
img_itk = itk.GetImageFromArray(array_new)
img_itk.SetOrigin(tuple(origin))
img_itk.SetSpacing(tuple(spacing))
itk.imwrite(img_itk, 'itk_sample2.mha')

array_new = np.moveaxis(array, [0,1,2], [2,1,0]) # [NOTE] No need to use .copy() here
img_sitk = sitk.GetImageFromArray(array_new)
img_sitk.SetOrigin(tuple(origin))
img_sitk.SetSpacing(tuple(spacing))
sitk.WriteImage(img_sitk, 'sitk_sample2.mha')

print (' - itk.imread(itk_sample1.mha) shape:', itk.GetArrayFromImage(itk.imread('itk_sample1.mha')).shape)
print (' - sitk.ReadImage(sitk_sample1.mha) shape:', sitk.GetArrayFromImage(sitk.ReadImage('sitk_sample1.mha')).shape)
print (' - itk.imread(itk_sample2.mha) shape:', itk.GetArrayFromImage(itk.imread('itk_sample2.mha')).shape)
print (' - sitk.ReadImage(sitk_sample2.mha) shape:', sitk.GetArrayFromImage(sitk.ReadImage('sitk_sample2.mha')).shape)