Partial Volume and (linear) interpolation

[romainVala]

Hi there

I just posted the same question on the itk forum, but I copy past it there, in case you have an other insight
(and because resample_from_to (like itk resample) produce weird results)

I am trying to construct a ground truth partial volume map (of brain region)
For instance I start with a high resolution (0.25 mm³) binnary mask of Gray Matter (GM) and I want to compute the partial volume at 1mm^3 resolution
Since this is a multiple of the initial resolution I can compute the exact partial volume, simply as the average of the original binary mask (average over the 444 voxel at 0.25 mm resolution)

I computed this partial volume using torch AveragePooling

Then I compared with the resampling method as implemented by torchio (but it rely on itk ResampleImageFilter from simple itk) with the linear interpolation I obtain a very different results

Actually it is kind of similar, but all values are only in the set ( [0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000])

Sorry if is is a bug on my side, but I double checked and I do start with a nifti f32 bit float volume. Since I average over 4^3 voxels I expect to have 64 possible values (| 0 0.0156, 0.0312 ect …]

I compared different implementation of interpolation:
_ resample_from_to from nibable give the same result as itk
_ mrgrid (from mrtrix) does give the exact same result a AveragePooling,
_ flirt gives me a similar results as mrgrid (but more partial volume !)

I am very surprise to have so much different results. How can it be that a simple linear interpolation leads to different implementation ?

My initial objective was to quantify the error when estimating the Partial volume with voxel size not being a multiple of 0.25, but I can not even get a clear answer for the simple case …

Many thanks for your insight

[effigies]

If your input values are all binary, then I would expect linear interpolation at voxel corners to give you the values you saw ($\{n/8 : n\in[0..8]\}$ ). The way that resampling works is that you take the location of the output voxel and map it into coordinates for the input space, then use your interpolation rule to determine how to find a value that does not line up exactly with one of your input voxels. In the case of linear, you will take a weighted average of the eight nearest voxels; on the corners, these weights will be equal, and so you can only end up with the values you found.

To reword it, with linear interpolation, you're averaging neighboring values at a point, not averaging all values in a volume. The resampling algorithm does not have any notion of volumes in the output space. It's not the same operation as average pooling unless you happen to be downsampling by a factor of two, with each output voxel containing exactly 8 whole input voxels. (So you could perform this twice, with 0.5mm^3 and 1mm^3, if you like.)

mrgrid's default oversampling seems to be what's doing the trick here. I would guess that it's increasing the number of voxels being averaged by the linear interpolation step, but I'm not positive how you would do that with the tools we use. It very well could be that it does a double-interpolation in this case. @Lestropie might be able to explain more.

I would guess that FLIRT is performing a Gaussian blur before resampling.

[romainVala]

Thanks for your input !
I do now better understand what interpolation is.

Indeed you are right if I compare now the interpolation from 0.25 to 0.5, I do find the same output for itk (nibable) and mrgrid. so performing it twice would lead me to the mrgrid result

Kind of weird though: our physical interpretation of dowsampling, is really to take an average within the new voxel value. doing an interpolation (downsampling) without taking into account the target voxel dimension seems quite strange to me. I know that mrtrix devs come from MRI physics, may be this is why, they implement a downsampling that is closer to what you expect with MRI (ie partial volume).

But I find it quite strange that two different tool perform very different output just for a linear interpolation ( is mrgrid performing an interpolation ?, or do they have a different interpretation of interpolation. ?). I'll ask to the mrtrix forum to comment here.

My new question is now : how can you obtain (best) estimate of partial volume at different voxel resolution. (from 0.25 to 0.9 for instance)
I first though interpolation would give me an answer, but I know understand this is quite different from partial volume information. (except may be with mrgrid ...)

thanks

[Lestropie]

When it comes to generating new intensity samples on a grid that is lower in resolution than the input, some form of antialiasing is required, whether it be smoothing of the input image prior to resampling or oversampling. So it's not "just an interpolation" in that scenario. mrgrid with the -info flag I would expect to provide some feedback on what is being performed here.

I do however want to explicitly pull the discussion away from interpolation. Because your question was not "how do i resample a high resolution image onto a low resolution voxel grid"; it was "how do I compute partial volumes". That's actually a different question, and one I've thought about a fair bit both historically and recently.

I can think of three cases worth considering.

1. Let's say for instance that your 1mm voxel grid aligns perfectly with your 0.25mm voxel grid, such that every 1mm voxel contains 64 0.25mm voxels. Perhaps what you actually want is, for every 1mm voxel, what fraction of the 64 encapsulated 0.25mm voxels are true in the input image. That's a different mathematical operation. It's also one for which I need to implement an MRtrix3 function at some point in the near future, though you could do it without having to write any C++ code.

2. Alternatively, consider the case where there's no simple relationship between input and output voxel grid. One thing I've been thinking about is how one might compute, for each output voxel, the volume of intersection with every input voxel. From this you could compute a weighted average. I like the idea, but I've not pursued it.

3. Yet another option exists in the scenario where you are not interested in the partial volumes of the input voxels and how these propagate to the output voxels; instead, you expect that some fraction of the output voxel is entirely one tissue, and the other fraction is entirely not that tissue. In this case you can take the notion of "partial volume" very literally. You can compute from the input binary mask a closed surface that covers the interface between true and false voxels, and then for the output image compute the fraction of the volume of that voxel that lies within the closed surface. MRtrix3 commands voxel2mesh and mesh2voxel do exactly this.

[romainVala]

thank you for the explanation.
I do not understand why resampling from high to low resolution and computing partial volumes are 2 different question. I guess our difference comes from what we understand by "resampling" ...
When I want to resample a high resolved volume to a lower one, I would like to compute a MRI volume that would be equivalent to an acquisition done at lower resolution. MRI physic answer is then simple : I need to average all high resolved voxel contained in my big voxel (if both grid are alligned). The reason why this is true is because the MRI signal within a given voxel is the weighted sum of signal from different tissue present in the voxel (the weights being the volume fraction)
does it make sense ?
It is probably specific to MRI, but my point is that we do not want to interpolate but we want to compute a weighted average, exactly as you explained in the point 2 .

About your first point, yes this is what I want, but I do not understand, why you need a new code. From what I tested mrgrid is already perfoming this operation. (that was the point of this issue: mrgrid is computing an average whereas nibabel and itk are just computing the interpolation). No need of C++ code to check .. I simply used torch AveragePolling and I did obtain the exact same result as with mrgrid. ... what do I missed ?
I tested with the -info flag, and the extra information I get is only (I used option -voxel 1)
mrgrid: [INFO] using oversampling factors [ 4 4 4 ]
does it mean you oversample the output image (1mm) to 0.25 mm and then do the average ?
if yes it does the job of point 1 no ?

you point 2 is where it becomes interesting, and more difficult. If I now run mrgrid with -voxel 0.9 I get the same info message (oversampling factors [4 4 4] but I do not understand what computation is done.
Since mrgrid seems to already perform the right operation (when the grid are allinged) I wonder how good it performs in case where the grids are not aligned ... ?
Note here that the weighted average will not be the exact result. Consider the case of 1mm grid, and two adjacent voxel. voxel 1 is containing 50% of white matter (Swm=1) and 50% of gray matter (Sgm=0.5) so its Signal value is 0.75, voxel 2 (neighbor in the x direction) is containing 100% of csf (wiht Scsf = 0.1) it signal is then 0.1. I now you resample the volume after a shift of 0.5 mm in the x direction. For the voxel in the new grid that is exactly containing half of voxel 1 and half of voxel 2 you will compute a value by taking the average (weights are 0.5) of the 2 voxels so (0.75+0.1)/2. But since this voxel is now containing 50% of gm and 50% of csf its true value is (0.5+0.1)/2 .
I guess this is a no solution issue : if you do not have the high resolution information, you can not perform a perfect interpolation.

I also like you point 3, my understanding is that it is very similar to whant you want to achieve in point 2. I used to play with it a few years ago, and I do not remember why I was not happy with voxel2mesh solution and I went for toblerone implementation. (if you are interested, I can try again and explain further)

Finally there may be a forth option, which is going with fourier transform. When dealing with translation and rotation perfoming it in the image domain (and thus using interpolation via resampling) is less optimal than computing it in the fourier domaine.
performing the resampling several times make a clear point that fourier interpolation does not get blured

Personally I am still impress by this results which make me thing we should go with fourier resampling ...

[Lestropie]

I do not understand why resampling from high to low resolution and computing partial volumes are 2 different question.

It's different under a specific interpretation of "partial volume", where one is interested specifically in the proportions of the voxel volume that are inside or outside of isocontours. This image conveys calculation of voxel partial volume based on an explicit surface representation, but bear in mind that it's possible to get from an image to a surface representation using the Marching Cubes method.

From what I tested mrgrid is already perfoming this operation.

Sorry, I likely misread this bit. Certainly some time in the past we relied on smoothing of the input image; it was perhaps changed to oversampling at some point.

If I now run mrgrid with -voxel 0.9 I get the same info message (oversampling factors [4 4 4] but I do not understand what computation is done.
Since mrgrid seems to already perform the right operation (when the grid are allinged) I wonder how good it performs in case where the grids are not aligned ... ?

My expectation is that it is taking the volume occupied by the output voxel, drawing a lattice grid of points within that volume in scanner space, and drawing a sample at each of those points using interpolation. In the case where each voxel in the low resolution grid perfectly encapsulates 64 high-resolution voxels, that lattice of points will lie perfectly at the centres of those 64 voxels, and so what you get is exactly the same as doing a 64-voxel aggregate; "interpolation" has no effect. But the exact same process is applied regardless of the voxel grid alignment; so for unaligned grids, it's taking 4x4x4 lattice-distributed sample points.

https://github.com/MRtrix3/mrtrix3/blob/3.0.4/core/adapter/reslice.h#L196-L212

So in that case, for my own application I was thinking of, you're correct in that I don't need to write any new C++ code. 🍾

I do not remember why I was not happy with voxel2mesh solution and I went for toblerone implementation

Thanks for raising it; found the paper, had not seen it previously. This is exactly what I had wanted to create when I first implemented mesh2voxel for ACT way back in ~2012, but I opted for the brute force sampling approach as it was good enough at the time.

Finally there may be a forth option, which is going with fourier transform.

Some of your conclusions may be based on the esoterics of the data you're dealing with, which look highly piecewise constant. In that scenario using oversampling but with a nearest neighbour interpolator might get you closer to what you want.

[romainVala]

Ok, just to summarize my understanding
Mrgrid resample interpolation is different from the one use in nibable or itk because mrgrid first oversample the target grid (by a multiple of 2) to be a resolution closest to the initial one, then it perform interpolation, and average to upsample back to the desire resolution.
@effigies don't you think this is a better way to interpolate

I probably miss something, but it seems to me that mrgrid interpolation is better at least when interpoling from high to low resolution

About case 2)

But the exact same process is applied regardless of the voxel grid alignment; so for unaligned grids, it's taking 4x4x4 lattice-distributed sample points.

https://github.com/MRtrix3/mrtrix3/blob/3.0.4/core/adapter/reslice.h#L196-L212

So in that case, for my own application I was thinking of, you're correct in that I don't need to write any new C++ code. 🍾

"my own application" refer here to your case2 (in you first post). But I don't think that the current implementation is doing the exact same thing. In other word, is the interpolation of "4x4x4" lattice distributes sample points equivalent to computing the volume intersection between the two grids ?
it is not equivalent since interpolation is base on the image intensities. May be it becomes equivalent if one apply it to a binary mask. (but then which interpolation method would give the exact result (partial volume) linear cubic sinc ?)

I do believe your case number 2 is worth to implement. My previous numerical example was a bit naive ( of course it won't be perfect for real images) but at least, if one compute it from a binary mask then this solution should give the exact partial volume. (and this is why it is worth to implement)

It's different under a specific interpretation of "partial volume", where one is interested specifically in the proportions of the voxel volume that are inside or outside of isocontours.

I fully agree I just try to work with a know case, so I suppose I have a "voxelique" isocontour at high resolution, (binary mask at 0.25 mm) and I want to compute the exact partial volume map at different resolution (and grid position)

Some of your conclusions may be based on the esoterics of the data you're dealing with, which look highly piecewise constant. In that scenario using oversampling but with a nearest neighbour interpolator might get you closer to what you want.

Indeed I show here example on synthetic images (ideal image of piecewise constant tissue), but I do not see why this bias the results ... ? ... I checked, I do obtain the same with a real data. The blurring is obvious with resampling in image domain whereas it is not in fourier domain. Unfortunately, The draw back is that I then get issue with gibbs artefact due to discret fourier
Applying five consecutive rotation (of 1 degree) with nearest neighbour does avoid bluring (by definition) but it induces large discontinuities in the image (funny results, I can put screenshot if interested)

[effigies]

@effigies don't you think this is a better way to interpolate

No. Interpolation is a well-defined operation, and adding oversampling and low-pass filtering (averaging) to it unconditionally is going to violate somebody's expectations. It is a reasonable thing for a CLI tool that is solving a particular problem to do, but not an interpolation function.

Just on the scope of nibabel: the focus is on IO and making structures amenable to manipulation with numpy/scipy, not on implementing neuroimaging tasks. There are projects like nipy, nilearn and dipy for more opinionated and targeted operations. nibabel.processing is a very thin layer for wrapping scipy.ndimage resampling methods, but for anything more complicated I would expect users to use those as templates and write their own.

Now, that isn't to say that we can't make it easier to do these things. I could imagine providing something like:

def oversample(img: ImgT, factor: int | tuple[int, int, int]) -> ImgT: ...

def apply(func: Callable[[np.ndarray], np.ndarray], img: imgT) -> ImgT:
    return img.__class__(func(img.get_fdata()), img.affine, img.header)

This would let you do something like:

from functools import partial
from scipy.ndimage import convolve

def myresampler(img):
    oversampled = oversample(img, ...)
    filtered = apply(partial(convolve, ...), oversampled)
    return resample_from_to(filtered, ...)

[romainVala]

sorry I can't resist to show with real data (and nearest neighbour)
montage.pdf

When performing nearest neighbour recursively, the error accumulate, and we have a funny mixt of fixed voxel and "moving" voxel, ...

blurring induced by interpolation is a well known fact, and this is why all preprocessing pipeline try to perform the resampling operation only once (combining affine and non linear deformation for instance). but I do not understand why fourier interpolation is not implemented in any software, despite its obvious advantage ... (and no black nor white magic here ! ;-)

[effigies]

I do not understand why fourier interpolation is not implemented in any software,

• 3drotate
• 3dvolreg

For Python, scipy has some Fourier filters. I would first see if the utility you want exists in scipy or scikit-image. If it does not, we could implement it here (feel free to open an issue to scope out the proposal, or open a PR if it's quick), but I think in general it would be good to get these tools into domain-agnostic packages where they can be easily applied to nibabel images.

[romainVala]

thanks @effigies and sorry if this post goes a bit beyond nibabel main focus,

I understand your point about nibabel scope, but I do not understand why fsl or mrtrix resample function would try to solve a particular problem : Resampling definition should not depend on the particular problem those tool are dedicated to ...

thanks to both explanations, I know better understand but I am still surprise that such low level operation (linear interpolation) is implemented differently in different software.

So mrgrid is doing both at the same time oversampling and interpolation. May be the important question now is : is it the correct solution when performing downsampling ?
If you agree that yes, then may be there should be an implementation in nibabel too ...

(Note that in current neuro-imaging pipeline we usually perform a resampling toward the higher resolution, so may be the exactness of dowsampling does not matter, but for training deep learning, it is common to use downsampling for data augmentation, and I would prefer to perform a more physical operation )

About fourier interpolation, thanks for pointing afni tools, but note that this is used for rotation (and translation) only, not for "general" resampling toward a given reference. But I could verify my implementation (which is also currently only for rotation too) : I do obtain the same results as afni, and I also agree with their documentation :
-Fourier = Use a Fourier method (the default: most accurate; slowest).

[Lestropie]

mrgrid first oversample the target grid (by a multiple of 2)

Doesn't have to be a multiple of 2; essentially a ceiling operation on the ratio of voxel spacings.

In other word, is the interpolation of "4x4x4" lattice distributes sample points equivalent to computing the volume intersection between the two grids?

In the case where the sampling lattice aligns perfectly with the input voxel grid, and you use an interpolator that yields the exact voxel value when sampling precisely at the centre of the voxel, yes it will be the same. As soon as that is violated, I presume no, but I've not attempted the math.

Unfortunately, The draw back is that I then get issue with gibbs artefact due to discret fourier

I would have thought that sinc interpolation with a reasonable kernel window would achieve similar but with less ringing?

Applying five consecutive rotation (of 1 degree) with nearest neighbour does avoid bluring (by definition) but it induces large discontinuities in the image

Yeah, that's expected; that's pretty much worst case scenario for that suggestion. But if the piecewise constant image you're demonstrating on is not actually a suitable exemplar for your intended application it's probably a moot point.