N4BiasFieldCorrection fails on images with negative values

[Shotgunosine]

Describe the bug
I'm running N4BiasFieldCorrection as part of fmriprep and it produced strange results on one of my images. The input image looks pretty much fine, but the output image has negative values pretty much everywhere inside the skull except for a ring and a small spot in the cerebellum with values in the 1000s. At the suggestion of @effigies I looked for negative values and found that there were a good number of voxels with values less than -100 with a minimum of -321. I added 500 to the input data and it worked.

To Reproduce
Steps to reproduce the behavior:

1. Download and unzip the data needed to reproduce from here: https://drive.google.com/file/d/1zaZ2Qug3Tt7_I6zQcPRCT5a-1mlTng-P/view?usp=sharing
2. Run N4BiasFieldCorrection --bspline-fitting [ 200 ] -d 3 --input-image ref_bold_corrected_trans.nii.gz --mask-image tpl-MNI152NLin2009cAsym_res-02_desc-brain_mask_trans_dil_hdr.nii.gz --output ref_bold_corrected_trans_corrected_debug.nii.gz to see the bad normalization result.
3. Run N4BiasFieldCorrection --bspline-fitting [ 200 ] -d 3 --input-image ref_bold_corrected_trans_embiggen.nii.gz --mask-image tpl-MNI152NLin2009cAsym_res-02_desc-brain_mask_trans_dil_hdr.nii.gz --output ref_bold_corrected_trans_corrected_embiggen_debug.nii.gz to produce the bettter results.

Expected behavior
It would be good if this corner case could either be caught and avoided. Alternatively, it could be caught and a warning could be provided.

Screenshots
Here's the input image:

And the mask overlaid on the input image in red:

The results of the bias correction. note the spot in the cerebellum:

Finally, the results after adding 500 to the input image:

Desktop (please complete the following information):
This is running inside a singularity container based on the fmriprep v 1.4.1 docker image. Here's the ANTs version information:
ANTs Version: 2.2.0.dev815-g0740f
Compiled: Jun 27 2017 17:39:25

[gdevenyi]

N4 includes a log transform so this is expected

[ntustison]

@gdevenyi , this is correct. Although this log transform is mentioned in the paper and documentation, should N4 check for negative values and issue a warning? And/or should we include an explicit mention of this issue in the help menu? If the former, is this something that should be added to the core N4 class in ITK? In addition to @gdevenyi , @blowekamp and/or @hjmjohnson -- do you have an opinion?

[gdevenyi]

In general, I need to wrap my usage of a number of ANTs tools to check for them producing negative values and chopping those to zero (interpolation artifacts in antsApplyTransforms for example). I'd be all for making an attempt to check for negatives and handle them in a consistent manner.

At least, I fully support adding some details to the help regarding log transform and negative values.

Side note: is N4(image + 500) == N4(image) + 500? Should we consider additional internal intensity transformations to standardize values before correction?

[Shotgunosine]

N4(image + 500) != N4(image) + 500

N4(image) gives uniform values of -162.4263 over the entire black region in those screenshots above. Adding 500 to that just makes the black region uniformly 337.5737.

The fmriprep pipeline is running N4 correction on the output of ANTS registrations. Those registrations appear to be the source of the negative values in this case.

[cookpa]

The negative values are likely the result of using BSpline interpolation to compute the deformed image.

I think it would be nice to have in the N4BiasFieldCorrection application a means to deal with negative values as well as winsorize outliers like antsRegistration does.

I think a lot of these issues got handled implicitly in the original N3 nu_correct program, because it converted things to MINC format and in doing so (usually) imposed a sensible range for the histogram.

[Shotgunosine]

Here's an example of a command that produces some rings of negative values around the brain:

antsRegistration --collapse-output-transforms 1 --dimensionality 3 --float 1 --initialize-transforms-per-stage 0 --interpolation LanczosWindowedSinc --output [ matching_transform, matching_transform_Warped.nii.gz ] --transform Translation[ 0.05 ] --metric Mattes[ bold_ref_uni_xform_masked.nii.gz, matching_uni_xform_masked.nii.gz, 1, 64, Random, 0.5 ] --convergence [ 500, 1e-07, 200 ] --smoothing-sigmas 8.0mm --shrink-factors 2 --use-estimate-learning-rate-once 1 --use-histogram-matching 1 --transform Rigid[ 0.01 ] --metric Mattes[ bold_ref_uni_xform_masked.nii.gz, matching_uni_xform_masked.nii.gz, 1, 64, Random, 0.5 ] --convergence [ 200, 1e-08, 100 ] --smoothing-sigmas 2.0mm --shrink-factors 1 --use-estimate-learning-rate-once 1 --use-histogram-matching 1 --winsorize-image-intensities [ 0.005, 0.998 ] --write-composite-transform 1

I can share the source files for that if needed.

[gdevenyi]

Thanks @Shotgunosine I'm referring to the more general case when there are already no negative values.

@cookpa indeed any interpolation other than nearest or linear often results in negative values which I have to deal with before passing forward in processing.

[cookpa]

Thanks @Shotgunosine @gdevenyi , you're correct, I should have said non-linear interpolation above.

I would hope that N4(image) and N4(image + constant) would be similar, but with a large enough constant I suppose it could shift the histogram in the log space sufficiently to cause problems.

[Shotgunosine]

Is there any chance that the registration and transformation could be updated with options to remove the negative values they introduce? Though maybe the most straight forward fix would be to have winsorize options in N4BiasFieldCorrection.

[gdevenyi]

fMRIprep should be handling this internally, or perhaps using Linear interpolation. But, as I understand, BOLD is not typically corrected for inhomogeneities, and you're supposed to correct/use the native space BOLD, not the template space images wherever possible.

I wonder whether BSpline has been validated for transforming BOLD without creating interesting effects.

[ntustison]

Just to answer a couple questions:

is N4(image + 500) == N4(image) + 500?

As @Shotgunosine pointed out, this is not true. But doing something like

0) get range of input image = [a,b]
1) rescale input image to be in a positive, e.g., [100,1000]
2) perform N4 on 1)
3) rescale image from 2) to be in [a,b]

would be reasonable although since the MRI scale is not absolute, step 3) would probably not be necessary. One thing about the core of this algorithm which goes back to N3 is that the range of the output could be very different from the range of the input. That's why we ultimately put in the -r option.

Is there any chance that the registration and transformation could be updated with options to remove the negative values they introduce?

No, this is feature creep which we try to avoid.

My current thought is to make the nonnegativity requirement explicit in the help menu and then throw a warning if there are negative values with the suggestion above included with that warning.

[cookpa]

Agree on not changing the registration, which is not just for MRI, or even just for medical images. You could register CT images where negative values are perfectly correct.

I don't understand why fmirprep wants to run N4 after resampling the image to template space, but doing it beforehand would solve the interpolation problem.

ImageMath's TruncateImageIntensity command might solve the problem as well as improve robustness to outlier intensities skewing the histogram.

I also notices the default shrink factor is being used, which is 4. With a low-res BOLD image you might be sampling way fewer voxels than you would for a T1. This might also increase susceptibility to outliers.

[gdevenyi]

I don't understand why fmirprep wants to run N4 after resampling the image to template space, but doing it beforehand would solve the interpolation problem.

I don't think fMRIprep does, I think @Shotgunosine does. fMRIprep did do a resample, but I don't think it does (or recommends) N4 on BOLD.

[Shotgunosine]

I'm not doing anything extra or special, this is all from within fMRIprep's workflows. It's happening because the original BOLD image is transformed and then skull stripping is repeated on the transformed image. N4BiasFieldCorrection is part of the BOLD skull stripping workflow. @oesteban may be able to comment further on how it's used.

[ntustison]

N4 can be pretty "aggressive" (relative to N3 or N3.5--N3BiasFieldCorrection) so often, if I remember correctly, @stnava will employ N3.5 for what might be considered non-traditional uses (e.g., bold or dwi).

[ntustison]

#824

[oesteban]

Thanks @Shotgunosine for the detailed description of the problem, and thank you all for the great suggestions.

I totally agree all the issues and recommendations given above pertain to fMRIPrep. I also think that a big warning when N4 finds negative values would be more than useful (and could've preemted this issue).