Resampling to Isotropic: Signal Processing Theory

Hi ITK Community!

Sanity check here about Nyquist-Shannon sampling theorem and resampling images to be isotropic. I have read through the online example [1]. I have also read the sinc approximation paper [2].

Setup I have CT data that is captured at 0.5 x 0.5 x 2.0 mm resolution. I want my data in 1.0 x 1.0 x 1.0 mm resolution for quantification purposes (FEA, actually).

Theory Following classsic DSP theory, I would do the following.

I apply a low pass filter (say Gaussian) in the X and Y directions at a spatial frequency of 0.5 mm^-1 (Sampling frequency is 1 mm^-1 so Nyquist frequency is 0.5 mm^-1).

Then, I would drop samples in X and Y and zero order hold in Z.

Then, I would apply a low pass filter in Z at 0.5 mm^-1.

This all makes sense from a classic DSP approach. I avoid aliasing in X and Y with prefiltering and I get a reconstructed filter with postfiltering in Z.

Interpolation Filter The problem arrises with the concept of interpolation in the ITK Resample filter.

In the online example [1], X and Y are blurred with a Gaussian filter and a linear interpolator is used for resampling. Letā€™s say I use a Hamming Windowed Sinc interpolator. I map a point from the isotropic image to the fixed image, sample by convolving the Hamming Windowed Sinc, and store that value. This effectively causes blurring in my input image, similar to the effect of prefiltering.

Question If I use a proper interpolation filter (3rd order BSpline, windowed sinc) I do not have to prefilter my image. I will have minimal aliasing because I am effectively prefiltering my image in the interpolator. Is this correct?

[1] https://itk.org/Doxygen/html/Examples_2Filtering_2ResampleVolumesToBeIsotropic_8cxx-example.html
[2] https://www.researchgate.net/profile/WJ_Niessen/publication/221401597_Quantitative_Comparison_of_Sinc-Approximating_Kernels_for_Medical_Image_Interpolation/links/0deec525e649cb45a6000000/Quantitative-Comparison-of-Sinc-Approximating-Kernels-for-Medical-Image-Interpolation.pdf

I think so.

Have a look at: LesionSizingToolkit/include/itkIsotropicResamplerImageFilter.hxx at master Ā· InsightSoftwareConsortium/LesionSizingToolkit Ā· GitHub

Thank you for the example! Here theyā€™re doing 3rd order BSpline interpolation, no prefiltering. This agrees with my understanding.

Thank You!

1 Like

A good explanation of the difference between prefiltering and interpolating, and what has been called ā€˜generalized interpolationā€™ is available here:
http://bigwww.epfl.ch/publications/thevenaz9901.pdf

ITK implements a version of the code that is here:
http://bigwww.epfl.ch/thevenaz/interpolation/

and is also explained in this paper:
Unser et al. https://doi-org/10.1109/78.193220

1 Like

What a good and complete description http://bigwww.epfl.ch/publications/thevenaz9901.pdf
Thanks @warfield.

We should create a sticky post holding high-quality papers/reviews like that on important topics.

Wow, thank you @warfield. That was a journey.

Iā€™m going to check my understanding here.

Interpolation Represent the continuous signal as the convolution of the discrete signal with an interpolation function.

Generalized Interpolation Represent the continuous signal as the convolution of a coefficient image and basis functions. The distinction is to calculate the coefficients from the image before interpolating (donā€™t assume the coefficients are the discrete signal).

Aliasing
From what Iā€™m reading, antialias depends on your implementation. If decimation is implemented (Unser et al. https://doi-org/10.1109/78.193220, section D. Least Squares Fitting), then this acts like the antialiasing filter.

So, in general, antialias filtering is needed. Interpolation provides a continuous representation for the discrete signal. Since you want to store that signal with fewer samples, you need to remove frequency content above the Nyquest sample frequency before building a representation (i.e. before resampling).

For the implementation in ITK, antialiasing filter is needed as per the example.

This raises concerns about how Iā€™ve been doing my registrations! I suppose you should always register to a high resolution image so you donā€™t need antialias filtering, among other reasons (fidelity, etc.)

P.S. I was mixing up the terms prefilteringand antialias filtering. Antialias filtering means to remove frequency content for discrete representation. Prefiltering, in the context of interpolation, means to determine the coefficients.

I was about to post a very similar topic when I found this discussion.
Many thanks for all the informative feedback and the pointers for further reading.
This leaves me with 3 remaining questions:

Regarding

For the implementation in ITK, antialiasing filter is needed as per the example.

does that mean that even when using IsotropicResamplerImageFilter instead of the ResampleImageFilter (as used in ResampleVolumesToBeIsotropic) it is still essential to apply low pass filtering (ā€œremove frequency content above the Nyquest sample frequencyā€) for the spatial directions to be down-sampled? As far as I understand the code of IsotropicResamplerImageFilter it mainly only differs by using BSpline interpolation instead of linear interpolation as in ResampleVolumesToBeIsotropic. Is that correct or am I missing something here?

So in other words it would NOT be appropriate to remove the low pass filtering from ResampleVolumesToBeIsotropic when replacing the linear interpolation with BSpline interpolation?

ITK implements a version of the code that is here: http://bigwww.epfl.ch/thevenaz/interpolation/

The ITK implementation referred to here is IsotropicResamplerImageFilter or ResampleVolumesToBeIsotropic?

And my original question:
What would be the correct way to generalize e.g. ResampleVolumesToBeIsotropic to resample arbitrary voxel spacing relations? For example if the final intended voxel spacing in x is close to the input voxel size in x I would think that the blurring should depend on the difference of the two. IsotropicResamplerImageFilter only seems to handle the case in which the desired output spacing equals the input spacing. E.g. in the case of a 1:10 voxel size ratio and a desired iso output spacing equaling to the z-spacing: What low pass filter should be used (Gauss, Hamming, etc) for xy and how to calculate appropriate parameters for the filter (e.g. sigma)?

As IsotropicResamplerImageFilter is in a remote module, ā€œITK implementsā€ probably does not refer to it.

Looking at the code of IsotropicResamplerImageFilter, it does not include smoothing so it needs to be applied separately.

1 Like

Thanks @dzenanz for your reply.
So I would conclude that ā€œITK implementsā€ refers to the BSpline interpolation implementation of ITK and not to ResampleVolumesToBeIsotropic (which seems to originate from @luisibanez from 2005 and was ever since using linear interpolation).
@warfield is that what you meant?

To summarize on a high level:

  • An anti-aliasing filter, e.g. Gaussian smoothing, is required to remove high frequencies that will cause artifacts when downsampling.
  • An accurate interpolator, e.g. windowed-sinc interpolator, is required to avoid artifacts when upsampling.
3 Likes

I meant that the BSpline interpolation uses an algorithmic strategy and code that is based on what Thevanez originally published.

1 Like

Many thanks for these clearifications @matt.mccormick and @warfield

So what I wonder then is how to choose the correct sigma for Gaussian smoothing when downsampling. ResampleVolumesToBeIsotropic uses the new pixel size directly:

However, I would think that is not a very appropriate choice in case the original pixel size is very similar, e.g. for downsampling from 0.9mm to 1mm pixel size. Would it be better to choose the difference of the original and new pixel size for sigma or would that be too small? Or should it be a function that depends on the difference?

Despite the argument that sophistication is not appropriate in the case of the example

and therefore linear interpolation is used

it might be better to change the example to use BSpline or windowed-sinc interpolator to be a better role model in general, as it seems this example is uses as a base for more general cases.

1 Like

This sounds much better to me than the target spacing. But ff some of the spacings are already close to your target, maybe it is better not to resample along that dimension. E.g. instead of resampling 0.9x0.9x1.2 into 1x1x1, maybe resample it into 0.9x0.9x0.9 or 1.2x1.2x1.2, or something of that sort.

Windowed-sinc interpolator is quite compute intensive, and is maybe not appropriate for on-the-fly execution. But it is probably the best choice for one-off resampling of all kinds. And as computation power has grown much more than image resolutions since that example was originally written, perhaps it is time to update it to use BSpline or windowed-sinc interpolator. Do you want to do it @grothausmann.roman?

1 Like

Ideally, sigma should be specified in each direction based on the physical sampling distance in the output image.

Thanks for all your feedback, perhaps we can discuss this easier when meeting at WBIR2020.

1 Like

Since WBIR2020 is postponed, continuing the discussion here.
I found this article from Cardoso M.J., Modat M., Vercauteren T., Ourselin S. (2015) Scale Factor Point Spread Function Matching: Beyond Aliasing in Image Resampling ā€“ MICCAI 2015 https://doi.org/10.1007/978-3-319-24571-3_81 PDF which provides the formula sigma= sqrt((isoSpacing^2 - inputSpacing[0]^2)/(2*sqrt(2*ln(2)))^2), which is not a function of the difference but based on difference and also leads to sigma= 0 in case input and output spacing are equal.
@dzenanz @matt.mccormick What do you think? Is its application appropriate in the discussed case, i.e. is its use in my resample-iso correct?

2 Likes

@grothausmann.roman :+1: looks great.

David Gobbi implemented an image resize filter in VTK that automatically applies anti-aliasing and interpolation as needed (it computes optimal anti-aliasing filter parameters). It would be nice if ITK would offer a similar feature. Ideally, it would be just a flag in itk::ResampleImageFilter to enable/disable anti-aliasing (to choose between higher-quality vs. faster resize).

2 Likes

I still have this on my todo list. But it is now likely I will get to this in the next few months.

2 Likes