sitk.Resample returns empty image

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?

Earlier cross posting and answer on SimpleITK github issues.

1 Like