# Advice on registering noisy head CT on template

**URL:** <https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849>\
**Category:** Uncategorized\
**Created:** [March 18, 2020, 6:47pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849 "2020-03-18T18:47:03Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![trypag](https://discourse.itk.org/user_avatar/discourse.itk.org/trypag/32/1271_2.png) [@trypag](https://discourse.itk.org/u/trypag)\
**Post date:** [March 18, 2020, 6:47pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/1 "2020-03-18T18:47:03Z")

</div>

Hello,

I am trying to register a head CT (with headset) on a head CT template. The presence of a headset is challenging the whole registration step (I think). It can eventually succeed if I wait long enough.  
I am looking for advice on how to make the process more robust and faster, based on my current configuration.

Fixed image config :  
Size: (512, 512, 55)  
Spacing: (0.3547742962837219, 0.3547742962837219, 2.438863754272461)  
Origin: (92.43216705322266, 83.35346984863281, -12.921518325805664)  
Direction: (-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0)

Moving image config :  
Size: (512, 512, 55)  
Spacing: (0.474609375, 0.474609375, 5.0)  
Origin: (-121.5, 14.5, -4.4000244140625)  
Direction: (1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)

 ![Screenshot from 2020-03-18 19-47-35](https://discourse.itk.org/uploads/default/original/2X/5/5ea94cf33d627dedada66deea94912b4af8a9b64.png)

```python
initial_transform = sitk.CenteredTransformInitializer(fixed_image, 
                                                      moving_image, 
                                                      sitk.Euler3DTransform(), 
                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)
# registration
registration_method = sitk.ImageRegistrationMethod()
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.1)
registration_method.SetInterpolator(sitk.sitkLinear)
registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=200)
registration_method.SetOptimizerScalesFromPhysicalShift()

# Set the initial moving and optimized transforms.
optimized_transform = sitk.Euler3DTransform()
initial_transform = sitk.CenteredTransformInitializer(sitk.Cast(fixed_image,moving_image.GetPixelID()), 
                                                      moving_image, 
                                                      sitk.Euler3DTransform(), 
                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)
registration_method.SetMovingInitialTransform(initial_transform)
registration_method.SetInitialTransform(optimized_transform)

```

The optimization process looks good to me, though in my opinion it seems pretty long regarding the simple task I want to solve.  
 ![Screenshot from 2020-03-18 19-41-07](https://discourse.itk.org/uploads/default/original/2X/9/93f8f5ee5133e75a3bbe0f22b39cc169e2afccb1.png)

I noticed that the fixed image is not in the same orientation as the moving image, I guess that having the same orientation will help a bit. I am not sure if the complexity to register comes from the parameters or the difference in the images.  
I would really appreciate some guidance 🙂

Thanks

---

<div class="post-metadata">

**Author:** ![blowekamp](https://discourse.itk.org/user_avatar/discourse.itk.org/blowekamp/32/79_2.png) [@blowekamp](https://discourse.itk.org/u/blowekamp)\
**Post date:** [March 18, 2020, 7:22pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/2 "2020-03-18T19:22:44Z")

</div>

Many additional features and ways to optimized and improve a registration are demonstrated in the SimpleITK Notebooks:

> **[InsightSoftwareConsortium/SimpleITK-Notebooks](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks)**
>
> Jupyter notebooks for learning how to use SimpleITK - InsightSoftwareConsortium/SimpleITK-Notebooks

Please look through Notebooks in the 60’s for registration tips.

---

<div class="post-metadata">

**Author:** ![trypag](https://discourse.itk.org/user_avatar/discourse.itk.org/trypag/32/1271_2.png) [@trypag](https://discourse.itk.org/u/trypag)\
**Post date:** [March 19, 2020, 12:58pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/3 "2020-03-19T12:58:03Z")

</div>

Thanks for the tutorial, adding Multi Resolution really helped to find a fast and good estimate of the final transform, love it !

I am now quite satisfied with the registration process, it seems accurate.  
Though the resampled moving images have very noticeable artifacts, as you can see mostly on the skull contours

 ![Screenshot from 2020-03-19 13-50-37](https://discourse.itk.org/uploads/default/original/2X/9/963decbdc0cf40dd0e778115759968d06550b839.png) ![Screenshot from 2020-03-19 13-48-47](https://discourse.itk.org/uploads/default/original/2X/9/9f8b435648c55942150e689e1b9debd2ca4884d6.png) ![Screenshot from 2020-03-19 13-45-13](https://discourse.itk.org/uploads/default/original/2X/7/7238ccb9cc7aa7ed8cebe55982af33101505ea30.png)

I guess the fixed image I am using has a higher Z resolution (more than twice actually) compared to the moving image. What would you advice in such situations ? resampling the fixed image to a lower depth resolution ? any ideas would be appreciated, thanks.

---

<div class="post-metadata">

**Author:** ![blowekamp](https://discourse.itk.org/user_avatar/discourse.itk.org/blowekamp/32/79_2.png) [@blowekamp](https://discourse.itk.org/u/blowekamp)\
**Post date:** [March 19, 2020, 1:10pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/4 "2020-03-19T13:10:06Z")

</div>

Your moving image has a 1:10 x/y:z spacing ratio. The data limits the quality. Using a higher order interpolator may help somewhat.

The best way to proceed will depend on your particular application, purpose and quantified goals. For example if you wanted to proceed with additional registration then you would not resample, but use the produced transform as the “MovingInitialTransform” for the next phase of registration.

---

<div class="post-metadata">

**Author:** ![trypag](https://discourse.itk.org/user_avatar/discourse.itk.org/trypag/32/1271_2.png) [@trypag](https://discourse.itk.org/u/trypag)\
**Post date:** [March 19, 2020, 1:24pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/5 "2020-03-19T13:24:22Z")

</div>

The only goal is to normalize all the CTs to have same head orientation, centering and scaling. It’s mostly a pre-processing step before deep learning algorithms and also to ease visual interpretation for the MD. I would rather use the original moving image spacing while preserving the ROI obtained from the registration.  
Then I would only keep the slices from the moving image that are contained in the fixed image physical space, without resampling to the fixed image resolution.

Does it makes sense with regard to the intended goal ?  
Thank you

---

<div class="post-metadata">

**Author:** ![zivy](https://discourse.itk.org/user_avatar/discourse.itk.org/zivy/32/1726_2.png) [@zivy](https://discourse.itk.org/u/zivy)\
**Post date:** [March 19, 2020, 1:46pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/6 "2020-03-19T13:46:48Z")

</div>

Hello @trypag,

Given your goal, deep learning preprocessing, I am not sure you need to register your data using the standard registration approach. I would first try a simpler route:

1. Identify bounding head bounding box in each of the scans (Otsu filter works reasonably well for this).
2. Compute the affine transformation between template bounding box (reference image) and each of the other boxes (center difference gives translation, size differences give scaling and orientation is assumed to be the same, but can be computed from axes of oriented bounding box if needed).
3. Resample all of the bounding box images onto the reference image.

This was the approach taken in [this Jupyter notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/master/Python/70_Data_Augmentation.ipynb) which dealt with data augmentation. Possibly also relevant for your work.

---

<div class="post-metadata">

**Author:** ![trypag](https://discourse.itk.org/user_avatar/discourse.itk.org/trypag/32/1271_2.png) [@trypag](https://discourse.itk.org/u/trypag)\
**Post date:** [March 19, 2020, 2:25pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/7 "2020-03-19T14:25:12Z")

</div>

Hello @zivy,

That’s right, deep learning is not so sensitive to translation and scaling, so I could go for a simpler way. I have already tried some parts of the data augmentation notebook and I might be going back to it in a few days, it’s really great and so under rated by the deep learning community !

Though affine registration does a pretty good job at focusing on what I am mostly interested in, the brain. Otsu filter would not help me to estimate this fine grained ROI. Some of my CT scans contain neck and shoulder parts, some not.  
What I would like to achieve first is to get rid of these parts so that my ROI is centered on the brain. Hope this makes sense again 🙂

---

<div class="post-metadata">

**Author:** ![zivy](https://discourse.itk.org/user_avatar/discourse.itk.org/zivy/32/1726_2.png) [@zivy](https://discourse.itk.org/u/zivy)\
**Post date:** [March 19, 2020, 2:29pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/8 "2020-03-19T14:29:24Z")

</div>

Hello @trypag,  
Given your explantation, data characteristics, the registration approach makes sense. We always try to use the simplest model, but life is complicated…

---

<div class="post-metadata">

**Author:** ![trypag](https://discourse.itk.org/user_avatar/discourse.itk.org/trypag/32/1271_2.png) [@trypag](https://discourse.itk.org/u/trypag)\
**Post date:** [March 19, 2020, 9:11pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/9 "2020-03-19T21:11:48Z")

</div>

I am puzzled by this last resampling step …

```python
moving_resampled = sitk.Resample(moving_image_true, fixed_image,
                                 final_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())

```

What I would like is the following :  
1/ Register the moving image on the fixed image in the physical space (DONE)  
2/ Use the fixed image ROI as the final ROI  
3/ Resample X:Y with the fixed image resolution AND leave the resampled image Z resolution to the original value

A first starting point would be to apply the transform without changing the original resolution (only use the fixed image ROI), then move to a solution where X and Y axis are sampled.

For 2/  
I could use `.TransformIndexToPhysicalPoint()` to find the fixed image extend in the physical domain, then transform these positions to the moving image domain based on `final_transform`. Am I wrong ?

For 3/  
No single idea about how to select which axis can be resampled or not, guidance would be appreciated again, thanks 🙂

---

<div class="post-metadata">

**Author:** ![zivy](https://discourse.itk.org/user_avatar/discourse.itk.org/zivy/32/1726_2.png) [@zivy](https://discourse.itk.org/u/zivy)\
**Post date:** [March 19, 2020, 10:06pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/10 "2020-03-19T22:06:39Z")

</div>

Hello @trypag,

This is readily achievable if you understand the SimpleITK/ITK tenet that an image is a physical grid in space. You just need to figure out how to describe various grids that occupy the same physical location. See [this notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/master/Python/21_Transforms_and_Resampling.ipynb) for a discussion on resampling.

Below is code that should do what you want (haven’t run it so it may have errors, but the idea is there):

```auto
fixed_spacing = fixed_image.GetSpacing()
moving_spacing = moving_image.GetSpacing()
           #spacing x-y from fixed, z from moving
new_spacing = fixed_spacing[0:2] + moving_spacing[2:3] 
fixed_size = fixed_image.GetSize()
          # compute the new image size using the new spacing so that it occupies the same
          # region in space as the fixed image.
new_size = fixed_size[0:2] + (int(0.5 + fixed_spacing[2]/moving_spacing[2]*fixed_size[2]) ,)

moving_resampled = sitk.Resample(moving_image_true, new_size, final_transform, sitk.sitkLinear, fixed_image.GetOrigin(), new_spacing, fixed_image.GetDirection(), 0, moving_image.GetPixelID())

```

---

<div class="post-metadata">

**Author:** ![trypag](https://discourse.itk.org/user_avatar/discourse.itk.org/trypag/32/1271_2.png) [@trypag](https://discourse.itk.org/u/trypag)\
**Post date:** [March 20, 2020, 6:43pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/11 "2020-03-20T18:43:49Z")

</div>

Hello @zivy,

Thank you so much for your help !  
I took some time to read the tutorial again and then your solution. Everything makes sense except there :

```auto
new_size = fixed_size[0:2] + (int(0.5 + fixed_spacing[2]/moving_spacing[2]*fixed_size[2]) ,)

```

I would rewrite it :

```auto
new_size = fixed_size[0:2] + (int(0.5 + (fixed_spacing[2]*fixed_size[2])/moving_spacing[2] ,))

```

Could you explain what does `0.5` stands for ?  
I understand that with `(fixed_spacing[2]*fixed_size[2])/moving_spacing[2]` we find the size of the moving image based on the size in the physical space, but I don’t get this 0.5 value.

Thanks !

---

<div class="post-metadata">

**Author:** ![zivy](https://discourse.itk.org/user_avatar/discourse.itk.org/zivy/32/1726_2.png) [@zivy](https://discourse.itk.org/u/zivy)\
**Post date:** [March 20, 2020, 7:00pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/12 "2020-03-20T19:00:13Z")

</div>

You’re welcome.

The 0.5 is in combination with the cast, we want the nearest integral value to a positive value `x` so I did `int(0.5+x)`. Likely the more Pythonic way is to do `round(x)`.

---

<div class="post-metadata">

**Author:** ![trypag](https://discourse.itk.org/user_avatar/discourse.itk.org/trypag/32/1271_2.png) [@trypag](https://discourse.itk.org/u/trypag)\
**Post date:** [March 20, 2020, 7:02pm UTC](https://discourse.itk.org/t/advice-on-registering-noisy-head-ct-on-template/2849/13 "2020-03-20T19:02:48Z")

</div>

Alright, got it 🙂  
Thanks again, the case is solved !