Resampling to Isotropic: Signal Processing Theory


(Bryce Besler) #1

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


(Pablo Hernandez-Cerdan) #2

I think so.

Have a look at: https://github.com/InsightSoftwareConsortium/LesionSizingToolkit/blob/master/include/itkIsotropicResamplerImageFilter.hxx


(Bryce Besler) #3

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

Thank You!


(Simon Warfield) #4

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


(Pablo Hernandez-Cerdan) #5

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.


(Bryce Besler) #6

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.