NF: Conformation function and CLI tool to apply shape, orientation and zooms

[kaczmarj]

closes #851
related to neuronets/kwyk#5

This PR adds a function to resample a volume to arbitrary shape and voxel sizes. This is meant to be a close replica of freesurfer's mri_convert --conform command. In freesurfer, the image is reoriented to LIA, but in this command, the image is reoriented to RAS.

So this command:

• resamples to a specified output shape
• resamples to a specified voxel size
• casts to uint8 (while also changing the range of data to [0, 255]).
• reorients to RAS

The heavy lifting is done by nibabel.processing.resample_from_to().

I have not added tests yet.

[codecov]

Codecov Report

Merging #853 into master will increase coverage by 0.03%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #853      +/-   ##
==========================================
+ Coverage   91.70%   91.73%   +0.03%     
==========================================
  Files          96       97       +1     
  Lines       12311    12359      +48     
  Branches     2173     2177       +4     
==========================================
+ Hits        11290    11338      +48     
  Misses        684      684              
  Partials      337      337              

Impacted Files	Coverage Δ
nibabel/affines.py	100.00% <100.00%> (ø)
nibabel/cmdline/conform.py	100.00% <100.00%> (ø)
nibabel/processing.py	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e2f50b4...2177a59. Read the comment docs.

[kaczmarj]

ready for review!

[yarikoptic]

minor comments.
Is actual (command line) "command" interface coming in this PR or later?

[yarikoptic]

should it become resilient to x_max == x_min when it would either set it to new_min or blow some better than ZeroDivisonError exception?

[kaczmarj]

thank you @yarikoptic - i will account for x_max == x_min and i will add a test for this case.

[yarikoptic]

is that what freesurfer does? I thought it might be nice to have this one applicable to 4D (or 5D - AFNI etc) series, to "conform" the first 3 dimensions to desired resolution/orientation etc

[kaczmarj]

i don't think it is what freesurfer does, but for simplicity and due to my lack of knowledge, i made it this way. how should the output shape and zooms be changed to support 4- and 5-D images?

[effigies]

The output shape and zooms would not be affected by rolling the first three axes.

In any event, I'm okay with getting this in and adding 4/5D image support later, but if we're going to have a test, then we should mark it as one that we intend to relax. Otherwise it may look like a desired constraint to future developers.

[kaczmarj]

Is actual (command line) "command" interface coming in this PR or later?

@yarikoptic - it could come in this pr. is there a preferred entrypoint? nib-conform?

[effigies]

So you're already not doing exactly what FreeSurfer does, so I don't think you need to stick too close to what it does.

I think you obviously want to be able to specify dimensions, zooms and orientation from the CLI. I would suggest that the dtype and rescaling might also be useful parameters. For rescaling, I would make it a boolean flag to indicate that the data should be rescaled to the full range of the dtype (True) or preserve values and adjusting scale parameters as needed (False).

Anyway, I haven't reviewed the code yet, but figured I might chime in about capabilities.

[kaczmarj]

For rescaling, I would make it a boolean flag to indicate that the data should be rescaled to the full range of the dtype (True) or preserve values and adjusting scale parameters as needed (False).

this makes sense for integer types, but i don't think it should be done for float types.

or perhaps we could make things much easier and not cast the data.

[effigies]

I agree that rescaling into the range for floats doesn't make sense.

[effigies]

Hi @kaczmarj just checking in on this. Do you want to aim for a 3.1 release (soon-ish...) or 3.2?

[kaczmarj]

hi @effigies - thanks for the ping. i will work on this tonight and will update. i aim for 3.1.

[pep8speaks]

Hello @kaczmarj, Thank you for updating!

Cheers! There are no style issues detected in this Pull Request. 🍻 To test for issues locally, pip install flake8 and then run flake8 nibabel.

Comment last updated at 2020-04-15 20:35:57 UTC

[kaczmarj]

@effigies - i added the command-line program nib-conform (including tests). i also narrowed the conform function. now, it conforms to a requested shape and voxel size. it does not change data type or transform data. how does that sound?

[kaczmarj]

Thanks @effigies. As far as I'm concerned, this PR should be good for final review. When merging, can you squash and merge? There are some commits in here are irrelevant and include incorrect code.

[effigies]

Will review tomorrow or Friday. Happy to squash, if that's what you prefer.

[effigies]

LGTM. @yarikoptic What do you think?

[effigies]

Bump @yarikoptic, in case this notification escaped, too.

[effigies]

@kaczmarj Sorry, forgot to squash. I think we'll survive a little mess.

[leej3]

@kaczmarj just to flag a gotcha for the freesurfer conform code... they seem to do image histogram normalization (shift values to center white matter somewhere around 100) and I believe they truncate to uint8 unless told not to do so. We got tripped up by that in attempting an identical but all python implementation. @patrick-mcclure may have some input on details of how to resolve it (though we may have just stuck with the freesurfer output?). It might be worth taking into account if this is for neural nets.

[patrick-mcclure]

@kaczmarj @leej3 For our successful projects, we, indeed, did stick to volumes pre-processed using Freesurfer's mri_convert.

[leej3]

Ok, thanks, good to know. And was the diagnosis that it was primarily the histogram normalization or were there other hypotheses?

[kaczmarj]

Thanks @leej3 and @patrick-mcclure. You are right that this does not perform histogram normalization or truncating to uint8. I think not truncating to uint8 should be fine, because that truncation will lose precision anyway. I'd be curious to hear if you think histogram normalization had a beneficial effect on the model performance.

[leej3]

Well you certainly can't mix and match:train on mri_convert output and predict on output that was conformed without normalization. My memory was that the normalization performed by freesurfer was important for heterogeneous datasets, i.e. it was essential for any broadly usable tool. Correct me if I'm wrong @patrick-mcclure.

[satra]

they seem to do image histogram normalization

@leej3 - is this done during mri_convert conforming step or later in the freesurfer process? my recollection is that this is done later because that requires getting at what white matter is. mri_convert is generally quite simple.

[patrick-mcclure]

From what I remember, training using normalized volumes and then testing on non-normalized volumes lead to a significant decrease in test performance compared to testing using the same volumes with normalization.

…

On Mon, Apr 20, 2020 at 9:51 AM john lee ***@***.***> wrote: Well you certainly can't mix and match:train on mri_convert output and predict on output that was conformed without normalization. My memory was that the normalization performed by freesurfer was important for heterogeneous datasets, i.e. it was essential for any broadly usable tool. Correct me if I'm wrong @patrick-mcclure <https://github.com/patrick-mcclure>. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#853 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACPK4HSNKQ7MZK7D54LAUR3RNRHPDANCNFSM4J2AI2LA> .

[kaczmarj]

@satra @leej3 - it looks like the histogram normalization happens during the conversion to uint8: https://github.com/freesurfer/freesurfer/blob/dev/mri_convert/mri_make_uchar.cpp

from the comments at the top of that file:

 * uses the Tal xform to find a ball of voxels that are mostly brain.
 * the top of the intensity histogram in this ball will then be white matter,
 * which allows us to center it at the desired value approximately (110)

[satra]

@patrick-mcclure and @leej3 - i think there are two things getting mixed up here.

1. what the conform step in freesurfer does, which is what this should implement
2. what is used for model training and what the model can be applied to.

let's keep 2 out of this discussion, since that relates to preprocessing for models that nibabel could not care less about.

[satra]

@kaczmarj - i would still like to know how mri_convert can determine where this ball is without doing any registration. i can conform a blank mri scan if i want to.

[leej3]

The mri_make_uchar tool compiled from the code @kaczmarj referenced is in fact run later, I believe here in recon-all.

It's coming back to me now. I think I had gotten confused by the output naming and assumed orig.mgz was the output of the first step of the pipeline with just a basic conform performed but it is in fact following the uint8 conversion (which does that normalization).

So summary is that mri_convert is indeed just performing a format conversion and with the -conform option it additionally performs operations equivalent to the ones here.

Agreed, 2. is another problem

[kaczmarj]

to add onto that @leej3 , mri_make_uchar is not used by mri_convert --conform. sorry for any confusion. it does look like it's used by recon-all as @leej3 noted.

mri_convert --conform will do the following things (i went through the code again):

1. read in volume
2. create an mri template with 1^3 mm voxels, size 256^3, and LIA orientation.
   • the header of the template indicates that the volume is uint8 (but the values have not been cast yet).
3. change data type to uint8 with MRIchangeType. this function does not depend on an atlas. it creates a 1000-bin histogram of the mri values and then scales the values to fit the limits of uint8 (it supports other integer types, too)

[leej3]

Thanks for detailing that.

[satra]

@kaczmarj - perhaps why it needs to create a histogram is relevant to implementing 3. since the basic scaling can be implemented just with nibabel's changing datatype function, which does the scaling on the range of the data.

[kaczmarj]

@satra - i am still trying to figure that out. i might need to translate the c++ code to python/numpy to give a good answer.

The linear transformation is:

val = dest_min + scale * (val - src_min);

where dest_min is the minimum of the output data type (0 for uint8), and val is the value being scaled.

other values are as follows:

src_min, src_max = limits(mri);
src_min = (float)bin * bin_size + src_min;
src_max = (float)bin * bin_size + src_min;
scale = (dest_max - dest_min) / (src_max - src_min);

the bin values in the calculation of src_min and src_max are different and come from for loops. i'm having trouble understanding what's going on there.

[kaczmarj]

i implemented relevant parts of MRIchangeType below for easier debugging. this is copied as closely as possible to the freesurfer code.

the histogram is used to modify the denominator of the scale function. more specifically, the histogram is used to update src_min and src_max. src_max is also updated according to the number of nonzero values.

import numpy as np

def change_type(x, f_low=0.0, f_high=0.999):
    # The defaults come from `mri_convert.c`
    x = np.asarray(x)
    N_HIST_BINS = 1000
    hist_bins = [0] * N_HIST_BINS
    src_min, src_max = x.min(), x.max()
    dest_min, dest_max = 0, 255  # for uint8

    bin_size = (src_max - src_min) / float(N_HIST_BINS)

    bin_threshold = N_HIST_BINS / 5.0
    # for in vivo
    bin_threshold = 10

    DZERO = lambda x: abs(float(x)) < 1e-15
    nonzero = 0
    for val in x.ravel():
        if not DZERO(val):
            nonzero += 1
        bin_ = int((val - src_min) / bin_size)
        if bin_ < 0:
            bin_ = 0
        if bin_ > N_HIST_BINS:
            bin_ = N_HIST_BINS - 1
        hist_bins[bin_] += 1

    # Modify `src_min`
    nth = int(f_low * x.size)
    n_passed = 0
    bin_ = 0
    while (n_passed < nth) and (bin_ < N_HIST_BINS):
        n_passed += hist_bins[bin_]
        bin_ += 1
    src_min = float(bin_ * bin_size + src_min)

    # Modify `src_max`
    nth = int((1.0 - f_high) * nonzero)
    n_passed = 0
    bin_ = N_HIST_BINS - 1
    while (n_passed < nth) and bin_ > 0:
        n_passed += hist_bins[bin_]
        bin_ -= 1
    src_max = float(bin_ * bin_size + src_min)
    
    assert src_max > src_min
    scale = (dest_max - dest_min) / (src_max - src_min)
    xnew =  scale * (x - src_min) + dest_min

    # This is not part of this function, but somewhere, there must be a
    # clipping operation, such as
    # xnew.clip(dest_min, dest_max)
    
    return xnew

[satra]

@kaczmarj - this looks like a straight up np.histogram operation. with the given keyword args. the only thing that would change would be src_min and src_max given the parameters f_low and f_high. and that would require the clipping bit.

and it seems f_low and f_high are controlling for number of elements to be discarded from the histogram in bins to determine src_min and src_max.

so should be easily reimplementable with np.histogram.

[satra]

@kaczmarj - this looks like a straight up np.histogram operation. with the given keyword args. the only thing that would change would be src_min and src_max given the parameters f_low and f_high. and that would require the clipping bit.

and it seems f_low and f_high are controlling for number of elements to be discarded from the histogram in bins to determine src_min and src_max.

so should be easily reimplementable with np.histogram.