[SOLVED] Difficulty getting registration to work with mask on moving image

Beginner Questions
,
1

Update 6 February: I got this to work. I confirm the mask on the moving image works.

Hello! I am having difficulty to get a registration with a mask on the moving image to work. I have attached the full code below, as well as the two images that I’m trying to register. One of the images, the moving image, has “bad” values that I want the registration to ignore.

I am trying to register with a translation transform. The code works fine when I do not apply any mask. I can also get it to work when applying a mask to the fixed image, with m.SetMetricFixedMask. But when I try to apply a mask to the moving image, with m.SetMetricMovingMask, the registration refuses to run. It tells me that “All samples map outside moving image buffer”, but I do not understand why. I checked the mask, and it has 0 values exactly where the bad values are, and 1 elsewhere.

Could someone give me a hint about what I’m not doing right?

fix.tif (842.7 KB)
mov.tif (837.1 KB)
mov_mask.png
metric_moving_mask.py (1.4 KB)

Here are the contents of metric_moving_mask.py:

import SimpleITK as sitk  
  
# unregistered images, in a type suitable for direct ingestion into the registration  
fix = sitk.ReadImage("fix.tif")  
mov = sitk.ReadImage("mov.tif")                 
mov_mask = sitk.ReadImage("mov_mask.png")   
                                        
# set up registration               
m = sitk.ImageRegistrationMethod()  
m.SetMetricMovingMask(mov_mask)                                                                               
m.SetMetricAsMattesMutualInformation()  
m.SetMetricSamplingStrategy(m.NONE)  
m.SetOptimizerAsGradientDescentLineSearch(2, 50, 1e-6)  
m.SetInterpolator(sitk.sitkLinear)   
m.AddCommand(  
    sitk.sitkIterationEvent,  
    lambda: print(  
        f"{m.GetOptimizerIteration():3}, {m.GetMetricValue():10.5f}, {m.GetOptimizerPosition()}"    
    ),  
)  
  
# execute registration   
t = sitk.TranslationTransform(2)  
m.SetInitialTransform(t)  
final_transform = m.Execute(fix, mov)  
print(f"final transform: {final_transform}")                                    
print(f"Optimizer stop condition: {m.GetOptimizerStopConditionDescription()}")  
                                                                     
# diagnostic image                                                   
resampler = sitk.ResampleImageFilter()                               
resampler.SetReferenceImage(fix)                                     
resampler.SetInterpolator(sitk.sitkLinear)  
resampler.SetTransform(final_transform)  
mov_r = resampler.Execute(mov)  
fix_8 = sitk.Cast(sitk.RescaleIntensity(fix), sitk.sitkUInt8)  
mov_r_8 = sitk.Cast(sitk.RescaleIntensity(sitk.Mask(mov_r, mov_mask)), sitk.sitkUInt8)                                       
im_registered = sitk.Compose(fix_8, mov_r_8, fix_8)                                      
sitk.WriteImage(im_registered, "diagnostic.png")
2

Samples are randomly taken in the fixed image (within the fixed mask). Then they are transformed to the moving image. If they fall outside of the image (or within bad region), they are discarded. If the moving mask’s “good” region is small, this situation frequently happens.

If you don’t have a fixed mask, you can swap fixed/moving designation and later invert the transform. If you have both fixed and moving mask, try setting the image with more restrictive mask as the “fixed”. This trick is harder to pull of with transforms that do not have a closed-form inverse.

3

Hi Dženan, thank you for your reply. Your reply made me try out a few more things, and I finally got it to work. I will document it here for future reference :slight_smile:

Somewhat surprisingly (to me, at least :smiley: ), when I change the optimizer method, I can get the code shown above to work on the images attached. For anyone interested in the solution, change the line that reads

m.SetOptimizerAsGradientDescentLineSearch(2, 50, 1e-6)

into

m.SetOptimizerAsPowell()

With this change, the code now converges nicely on the right solution (you can inspect the diagnostic image created at the end). You will also need to change the line

mov_r_8 = sitk.Cast(sitk.RescaleIntensity(sitk.Mask(mov_r, mov_mask)), sitk.sitkUInt8)

into

mask_r = resampler.Execute(mov_mask)
mov_r_8 = sitk.Cast(sitk.RescaleIntensity(sitk.Mask(mov_r, mask_r)), sitk.sitkUInt8)

Of course, this begs the question why the gradient descent method failed. I suspect that it is trying extreme values of the parameters, but without extra information from the optimizer, I cannot confirm that.

Anyway, I consider this problem solved. Thank you for pointing me in the right direction!

1 Like
4

The gradient based methods are sensitive to the scales. They can over step the image bounds when not set correctional.

I don’t see any “SetOptimizerScales…” methods being used? You can give those methods a try.

1 Like
5

That is correct. For a pure translation, I didn’t think settings the scales (which I assume are 1 by default) would provide any benefit. Am I mistaken?

6

This is a question best answered by an experiment. The auto scaling scale takes into consideration the pixel spacing, and the intensity values, and the particular metric, which all can effect the gradient.

1 Like