WIP: Estimate inverse warp field when invert is set in MINC Grid XFM

[gdevenyi]

Currently, the handling for the TransformMINC Grid nonlinear transforms fail to handle files with Invert_Flag = True; in them due to the following failure:

This line:
https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/IO/TransformMINC/include/itkMINCTransformIO.hxx#L163
Results in this failure in antsApplyTransforms:

terminate called after throwing an instance of 'itk::ExceptionObject'
  what():  /opt/quarantine/ANTs/git/src/ANTs-build/ITKv5-install/include/ITK-5.0/itkDisplacementFieldTransform.hxx:77:
itk::ERROR: DisplacementFieldTransform(0x3782ed0): No displacement field is specified.
Aborted

In this patch I add a itkIterativeInverseDisplacementFieldImageFilter stage to estimate an inverted displacement field and provide it to the displacement field input.

I originally tried this with InverseDisplacementFieldImageFilter however the filter update step caused a segfault with no exception thrown so I'm currently unable to debug any further. Another alternative implementation would use FixedPointInverseDisplacementFieldImageFilter however it is currently an external module so I don't want to depend on it.

This was tested by:

• generating a forward-inverse Warp field pair in XFM format from antsRegistration
• applying the forward transform to the input image
• text-editing the XFM inverse file to add Invert_Flag = True;
• applying the modified inverse XFM to the input file
• minccmp shows xcorr: 0.9999058331

Paging @vfonov

[gdevenyi]

This patch addresses half of the problem cited at #147 (comment)

The other half will be added in that patch (flipping vector field directions)

[gdevenyi]

Hi @vfonov this fix is independent of the RAS<->LPS work, it estimates the inverse transform via a method similar to grid_inverse_transform_point_with_input_steps from grid_transform.c

[vfonov]

1. How do you determine the domain and sampling where inverse transformation is calculated?
2. how long does it take now to open .xfm file that has .inverse flag set ?

[gdevenyi]

Hi @vfonov,

The default for the filter is to create a domain and sampling identical to the input deformation field.

The time added for a 1mm iso nlin registration invert for the default 5 iterations is ~7ms.

[vfonov]

well in my experiment, running InvertDisplacementFieldImageFilter with default parameters on a a 2x2x2 grid file produced by MINC longitudinal pipeline takes about 1 min. Is IterativeInverseDisplacementFieldImageFilter that much faster and still produces appropriate result?
Running minc-based inverse on a 2x2x2 grid file takes 20sec.

[gdevenyi]

Sorry, follow-up here, a bit of a mixup on the timings, I was using the internal timing from the Print() functionality of the IterativeInverse function, which seems to be the mean of some loop within the code, rather than the true value. I annotated the code with a itk::TimeProbe start(); and stop(); for the filter->Update();

Here are the results for a 1x1x1 mm deformation field generated from antsRegistration between the GAD and MNI model mentioned in the RAS/LPS discussion, applied using antsApplyTransforms, this is specifically the time spent in the inversion filter.

Settings for the filters,
InterativeInverse (settings roughly copied from Invert filter):

inverter->SetNumberOfIterations(50);
inverter->SetStopValue(0.1);

Invert (settings based on the ITK example):

inverter->SetMaximumNumberOfIterations( 50 );
inverter->SetMeanErrorToleranceThreshold( 0.001 );
inverter->SetMaxErrorToleranceThreshold( 0.1 );
inverter->SetEnforceBoundaryCondition( false );

InverativeInverse: Time elapsed: 8.78565
InvertDisplacement: Time elapsed: 3.48849

minccmp comparing the transformed images made with the "true" inverse field (generated by antsRegistration) and the inverted fields:
Iterative:

file[0]:      true_inverse.mnc
mask file:    (null)
file[1]:      Iterative.mnc
ssq:          597264.9631
rmse:         0.2646113505
xcorr:        0.9999792877
zscore:       0.00321001647

Invert:

file[0]:      true_inverse.mnc
mask file:    (null)
file[1]:      Invert.mnc
ssq:          32444.84004
rmse:         0.06167338255
xcorr:        0.9999988759
zscore:       0.0008179832775

Happy to switch to the Invert filter based on it being faster and more accurate.

[vfonov]

So, it looks like it works fast on a "good" deformation field from ANTs , how long would it take on a not so smooth field from minctracc?

[vfonov]

Here,

this is output of CIVET pipeline on SAMIR-99 data sample. Let's see what the proposed code will do.
001_nlfit_It_inv.tar.gz

[gdevenyi]

Thanks @vfonov

I generated nlpfit outputs (2x2x2) using the same GAD-MNI pairs and got the following results:
Invert: Time elapsed: 0.107751
IterativeInverse; Time elapsed: 0.836424

And generating inverse field files from the SAMIR-99 data you sent with antsApplyTransforms -o [output.mnc,1]
Invert: Time elapsed: 0.967698
Iterative: Time elapsed: 1.48106

[gdevenyi]

Updated PR to use Invert filter instead of IterativeInverse filter.

[vfonov]

is the result similar to mincresample ?

[gdevenyi]

I don't have/couldn't find a copy of SAMIR-99 so I couldn't directly compare resamples, but for the GAD-MNI pair from nlpfit:

file[0]:      itk_Invert.mnc
mask file:    (null)
file[1]:      forward-mincresample.mnc
ssq:          1359805581
rmse:         12.62592949
xcorr:        0.9974990007
zscore:       0.01719338444

Its not quite as nice a comparison as the antsRegistration one since minctracc doesn't produce a forward/reverse pair that I can use to do a more-head-to-head comparison, but everything looks fine.

[vfonov]

technically it is pretty trivial to see if the inverse is good - i.e X * inv(X) = identity
but here is the scan

001_t1_tal.tar.gz

[gdevenyi]

Thanks @vfonov, an excellent suggestion re: identity.

Here, I concat (via added two transforms to the text of the .xfm) the forward and reverse transforms and apply those transforms with mincresample, and antsApplyTransforms with the itkInvertDisplacementField.

I then compare these outputs using minccmp against the original image since this operation should be the identity operation:
mincresample:

file[0]:      identity.mnc
mask file:    (null)
file[1]:      001_t1_tal_1.0.mnc
ssq:          1415726556
rmse:         14.11176214
xcorr:        0.9969889137
zscore:       0.01226894951

antsApplyTransforms:

file[0]:      001_t1_tal_1.0.mnc
mask file:    (null)
file[1]:      itk_invert_identity.mnc
ssq:          13572824.13
rmse:         1.381741078
xcorr:        0.9999711985
zscore:       0.002255303473

Look like the InvertDisplacement filter does a good job, and does so in 0.973168 seconds according to my instrumentation of the code.

[vfonov]

ok, let's merge it into the mainline ITK.
Would be good to include identity test into the testing code.

[gdevenyi]

+1 on a test. I'll look at how to add one and whether some of the existing example files could be used.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[vfonov]

We still need a test code to check if concatenation of inverted transform with itself is identity.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[gdevenyi]

This code is valid, I still need to add tests