below hola.py with itk(simpleitk) not working AT ALL.
TODO: please replace with Simple Elasix
ref:
https://simpleelastix.readthedocs.io/AffineRegistration.html
please do also checkout VoxelMorph:
https://github.com/voxelmorph/voxelmorph
sample script
| ''' | |
| ref | |
| https://simpleitk.readthedocs.io/en/v1.2.0/Examples/ImageRegistrationMethod1/Documentation.html | |
| https://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/Python_html/60_Registration_Introduction.html | |
| ''' | |
| import SimpleITK as sitk | |
| import sys | |
| import os | |
| raise valueError("not great, go use Elastix!") | |
| fixed_file = sys.argv[1] | |
| moving_file = sys.argv[2] | |
| transform_file = sys.argv[3] | |
| moved_file = sys.argv[4] | |
| kind = sys.argv[5] | |
| fixed = sitk.ReadImage(fixed_file, sitk.sitkFloat32) | |
| moving = sitk.ReadImage(moving_file, sitk.sitkFloat32) | |
| R = sitk.ImageRegistrationMethod() | |
| R.SetMetricAsMeanSquares() | |
| R.SetOptimizerAsRegularStepGradientDescent(4.0, .01, 200 ) | |
| if kind == 'affine': | |
| R.SetInitialTransform(sitk.AffineTransform(fixed.GetDimension())) | |
| else: | |
| R.SetInitialTransform(sitk.TranslationTransform(fixed.GetDimension())) | |
| R.SetInterpolator(sitk.sitkLinear) | |
| outTx = R.Execute(fixed, moving) | |
| sitk.WriteTransform(outTx, transform_file) | |
| resampler = sitk.ResampleImageFilter() | |
| resampler.SetReferenceImage(fixed) | |
| resampler.SetInterpolator(sitk.sitkLinear) | |
| resampler.SetDefaultPixelValue(100) | |
| resampler.SetTransform(outTx) | |
| out = resampler.Execute(moving) | |
| sitk.WriteImage(out,moved_file) | |
| #simg1 = sitk.Cast(sitk.RescaleIntensity(fixed), sitk.sitkUInt8) | |
| #simg2 = sitk.Cast(sitk.RescaleIntensity(out), sitk.sitkUInt8) | |
| #cimg = sitk.Compose(simg1, simg2, simg1//2.+simg2//2.) | |
| #sitk.Show( cimg, "ImageRegistration1 Composition" ) | |
| """ | |
| python hola.py \ | |
| /dingo_data/cvib-airflow/RESEARCH/10128/stp-suv/10128_ILD_004_1/voxelmorph/cropped_hrct.nii.gz \ | |
| /dingo_data/cvib-airflow/RESEARCH/10128/stp-suv/10128_ILD_004_1/voxelmorph/cropped_wb_ct.nii.gz \ | |
| /dingo_data/cvib-airflow/RESEARCH/10128/stp-suv/10128_ILD_004_1/voxelmorph/T.tfm \ | |
| /dingo_data/cvib-airflow/RESEARCH/10128/stp-suv/10128_ILD_004_1/voxelmorph/sitk_moved_wb_ct.nii.gz | |
| """ |
| def elastix_register_and_transform(fixed_file,moving_file, | |
| moved_file=None,moving_list=[]): | |
| fixed = sitk.ReadImage(fixed_file) | |
| moving = sitk.ReadImage(moving_file) | |
| elastixImageFilter = sitk.ElastixImageFilter() | |
| elastixImageFilter.SetOutputDirectory('/tmp') | |
| elastixImageFilter.SetFixedImage(fixed) | |
| elastixImageFilter.SetMovingImage(moving) | |
| elastixImageFilter.SetParameterMap(sitk.GetDefaultParameterMap("affine")) | |
| elastixImageFilter.LogToConsoleOn() | |
| elastixImageFilter.Execute() | |
| if moved_file: | |
| sitk.WriteImage(elastixImageFilter.GetResultImage(),moved_file) | |
| for moving_file in moving_list: | |
| if moving_file.endswith('ct.nii.gz'): | |
| out_pixel_value = -1000 | |
| ismask=False | |
| elif moving_file.endswith('lung.nii.gz'): | |
| out_pixel_value = 0 | |
| ismask=True | |
| elif moving_file.endswith('suv.nii.gz'): | |
| out_pixel_value = 0.0 | |
| ismask=False | |
| else: | |
| raise NotImplementedError() | |
| og_obj = sitk.ReadImage(moving_file) | |
| transformixImageFilter = sitk.TransformixImageFilter() | |
| transformixImageFilter.SetOutputDirectory("/tmp") | |
| transformixImageFilter.SetMovingImage(og_obj) | |
| transform_tuple = elastixImageFilter.GetTransformParameterMap() | |
| transform = list(transform_tuple)[0] | |
| transform['DefaultPixelValue']=[str(out_pixel_value)] | |
| if ismask: | |
| transform['FinalBSplineInterpolationOrder']=["0"] | |
| transform["ResultImagePixelType"] = ["int"] | |
| transform_tuple = (transform,) | |
| transformixImageFilter.SetTransformParameterMap(transform_tuple) | |
| transformixImageFilter.LogToConsoleOn() | |
| transformixImageFilter.Execute() | |
| moved = transformixImageFilter.GetResultImage() | |
| moved = sitk.Cast(moved,og_obj.GetPixelID()) | |
| sitk.WriteImage(moved,moving_file) |