I am trying to implement an atlas based segmentation model; however, when I resample after registration, I always get an image full of zeros.
Registration function (input images are non empty images of size 240,240,155 and 256, 256, 261 in numpy order:
def register_test3(fixed, moving):
initialTx = sitk.CenteredTransformInitializer(fixed, moving,
sitk.AffineTransform(
fixed.GetDimension()))
R = sitk.ImageRegistrationMethod()
R.SetShrinkFactorsPerLevel([3, 2, 1])
R.SetSmoothingSigmasPerLevel([2, 1, 1])
R.SetMetricAsJointHistogramMutualInformation(20)
R.MetricUseFixedImageGradientFilterOff()
R.SetOptimizerAsGradientDescent(learningRate=1.0,
numberOfIterations=300,
estimateLearningRate=R.EachIteration)
R.SetOptimizerScalesFromPhysicalShift()
R.SetInitialTransform(initialTx)
R.SetInterpolator(sitk.sitkLinear)
R.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(R))
R.AddCommand(sitk.sitkMultiResolutionIterationEvent,
lambda: command_multiresolution_iteration(R))
outTx1 = R.Execute(fixed, moving)
print("-------")
print(outTx1)
print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
print(f" Iteration: {R.GetOptimizerIteration()}")
print(f" Metric value: {R.GetMetricValue()}")
displacementField = sitk.Image(fixed.GetSize(), sitk.sitkVectorFloat64)
displacementField.CopyInformation(fixed)
displacementTx = sitk.DisplacementFieldTransform(displacementField)
del displacementField
displacementTx.SetSmoothingGaussianOnUpdate(varianceForUpdateField=0.0,
varianceForTotalField=1.5)
R.SetMovingInitialTransform(outTx1)
R.SetInitialTransform(displacementTx, inPlace=True)
R.SetMetricAsANTSNeighborhoodCorrelation(4)
R.MetricUseFixedImageGradientFilterOff()
R.SetShrinkFactorsPerLevel([3, 2, 1])
R.SetSmoothingSigmasPerLevel([2, 1, 1])
R.SetOptimizerScalesFromPhysicalShift()
R.SetOptimizerAsGradientDescent(learningRate=1,
numberOfIterations=300,
estimateLearningRate=R.EachIteration)
R.Execute(fixed, moving)
compositeTx = sitk.CompositeTransform([outTx1, displacementTx])
""" print("-------")
print(displacementTx)
print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
print(f" Iteration: {R.GetOptimizerIteration()}")
print(f" Metric value: {R.GetMetricValue()}")
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(fixed)
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
resampler.SetDefaultPixelValue(100)
resampler.SetTransform(compositeTx)
out = resampler.Execute(mask)"""
return compositeTx #, out
And I resample the it to the mask as follows:
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(moving_image_mask)
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
resampler.SetDefaultPixelValue(100)
resampler.SetTransform(q)
out = resampler.Execute(fixed_image)
No matter why I do, it does not seem to work (switching atlas and the target image, resampling with linear interpolator does not work either). The transform q in the resampler is the output of the registration function (compositeTx). What am I doing wrong here?