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?
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?
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)?
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?
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.
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?
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_81PDF 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?
looks great.