I’m a new user of simple ITK and I’m trying to register CT images and high resolution images from a camera.
When I run the registration with the following code, I get the error that joint pdf summed to zero:
registration_method = sitk.ImageRegistrationMethod()
# Similarity metric settings.
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=100)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.1)
registration_method.SetInterpolator(sitk.sitkLinear)
# Optimizer settings.
registration_method.SetOptimizerAsGradientDescent(learningRate=registration_method.EachIteration, numberOfIterations=100, convergenceMinimumValue=1e-6, convergenceWindowSize=10)
registration_method.SetOptimizerScalesFromPhysicalShift()
# Setup for the multi-resolution framework.
registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])
registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,0])
registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
# Don't optimize in-place, we would possibly like to run this cell multiple times.
registration_method.SetInitialTransform(initial_transform, inPlace=False)
# Connect all of the observers so that we can perform plotting during registration.
registration_method.AddCommand(sitk.sitkStartEvent, start_plot)
registration_method.AddCommand(sitk.sitkEndEvent, end_plot)
registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, update_multires_iterations)
registration_method.AddCommand(sitk.sitkIterationEvent, lambda: plot_values(registration_method))
final_transform = registration_method.Execute(fixed_image, moving_image)
By commenting the smoothing and shrinking factors, I don’t get the error, but the registration results are not really good.
Can someone explain me what are those factors and how I should choose them?
Before I answer your specific question, we need to understand how you initialize your registration. How is the value of initial_transform obtained, and can you visually validate that after this initialization the images have a reasonable overlap. Once you confirm that your initialization is reasonable, you can move to the second step which is actually doing the registration. If initialization is significantly off: (1) no overlap between image content or (2) contents overlap but significant rotation (e.g. 180 degrees) the registration will likely fail. See this notebook for a discussion of registration initialization.
Now that you have a good initialization, the shrinking and smoothing parameters define an image pyramid. The original image is first smoothed, deal with aliasing by removing high frequency components, and then resampled based on the shrinking factors. The pyramid scheme makes the registration more robust to initial misalignments (one voxel in the highest level of the pyramid is much larger than one voxel in the lowest level of the pyramid), and it often speeds up the process.
Finally, please note that plotting the metric values during registration is useful for debugging purposes but significantly slows down the process due to the time it takes to plot. In your final code you can either omit it (remove all the calls to AddCommand) or have them under a conditional flag that is set beforehand:
if debug_on:
registration_method.AddCommand(...)
...
The images have a small rotation over the z axis and also a small rotation over the x axis. Nonetheless, when I perform the registration as mentioned in the previous comment, no changes are apparently made:
The “joint PDF summed to zero” error usually means that there is no overlap between the images (registration is failing - taking a step that is too large in parameter space). Possibly at the high level of the pyramid the image is very small (~100x100) and this allows for such issues. This is a combination of image size, image content and sampling strategy.
The image content appears to be challenging, to me this looks like a cylinder with very small protrusions in three places. This makes registration a challenging task as rotation around the cylinder’s axis would barely have an effect on the similarity metric (the small protrusions are out-weighted by the cylinder when doing random sampling).
Based on your image, you may not need to use the multi-resolution approach, possibly try the single resolution approach (comment out the three lines dealing with the multi-resolution approach). Additionally if my understanding of the image content is correct you may want to consider a sampling strategy of None, change the code to:
I suspected that the protrusions might not be enought for a significant change in the metric value, but I don’t know if it is actually the case.
I ran the algorithm without the multi-resolution framework and no sampling strategy (I imagine the same as 100% sampling) and got the same results as before.
Besides that, my objects have a number that identify then in the top layers of the images (total of layers: 233, 217 - 233 are the number). I croped the image to contain only the number and conducted the same algorithm, getting the following results:
The results look reasonable to me, the registration is going in the right direction but it is likely stopped prematurely. The stopping condition is reported as “Maximum number of iterations (100) exceeded.” . The next thing to try is a larger number of iterations. The number of iterations stopping criterion is used to prevent extremely long registration times with minimal improvement in accuracy. When the registration is fast enough then you don’t need to set this to a tight bound, just give it a large number of iterations and it will likely stop before reaching this bound because it converges.
With respect to a sampling of NONE vs. RANDOM+100% sampling rate, these are not equivalent. The sampling of NONE uses all image voxels. The random sampling with 100% samples the image N times (N is the number of voxels), this will be different from using all voxels as it is sampling with replacement. In theory you could have a sample where the same voxel was repeated N times (the probability for this event is equivalent to the probability of independently choosing any other N samples, which isn’t zero).
Thanks very much for your help. It indeed worked for a larger number of iterations when I was considering just the number. Considering the whole component, it did not work. Probably because of the big overlapping area of the cylinder compared to the small protrusions.
Now the difference of no sampling and random sampling with sampling rate of 100% is also clearer to me.
