Resulting warp unexpectedly has a smaller field of view than input images

bcdarwin commented 4 years ago

Describe the problem

ANTS produces a warp field that is smaller than the input images, so the resampled moving image is not fully aligned to the fixed image.

To Reproduce

Make some random data using numpy and nibabel and run a trivial registration of the image to itself:

$ python3 -c "import nibabel as nib; import numpy as np; nib.save(nib.nifti1.Nifti1Image(np.random.random((500, 500, 500)), np.eye(4) * 2), '/tmp/img.nii')"

$ ANTS 3 --number-of-affine-iterations 0 -m 'CC[/tmp/img.nii,/tmp/img.nii,1.0,3]' -t SyN[0.2] -r Gauss[2,1] -i 2x2x1x1 -o /tmp/ants.nii
 Run Reg 
 values 1
  Fixed image file: /tmp/img.nii
  Moving image file: /tmp/img.nii
Metric 0:  Not a Point-set
  Fixed image file: /tmp/img.nii
  Moving image file: /tmp/img.nii
  similarity metric weight: 1
  Radius: [3, 3, 3]
  radius: [3, 3, 3]
Use identity affine transform as initial affine para.
aff_init.IsNull()==1
Use identity affine transform as initial fixed affine para.
fixed_aff_init.IsNull()==1
Continue affine registration from the input
affine_opt.use_rotation_header = 0
affine_opt.ignore_void_orgin = 0
transform_initial: IsNotNull():0
OptAffine: metric_type=AffineWithMutualInformation
MI_bins=32 MI_samples=32000
number_of_seeds=0 time_seed=1599255771
number_of_levels=1
number_of_iteration_list=[0]
graident_scales=[1,1,1,1,1,1,1,1,1,1,0.0001,0.0001,0.0001]
is_rigid = 0
mask null: 1
maximum_step_length=0.1
relaxation_factor=0.5
minimum_step_length=0.0001
translation_scales=0.0001
opt.transform_initial.IsNull(): 1
 opt.use_rotation_header: 0
 opt.ignore_void_orgin: 0
input affine center: [-499.008, -499.005, 498.986]
input affine para: [0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]
level 0, iter 0, size: fix[500, 500, 500]-mov[500, 500, 500], affine para: [0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]
    does not reach oscillation, current step: 0.1>0.0001
A=1 0 0
0 1 0
0 0 1

rotation R1 0 0
0 1 0
0 0 1

upper R1 0 0
0 1 0
0 0 1

s=0.25 u=0 v=0 w0 r=1
m_Rotation from vnl0 0 0 1
level 0, iter 0, size: fix[500, 500, 500]-mov[500, 500, 500], affine para: [0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]
    does not reach oscillation, current step: 0.1>0.0001
 v1 0 v2 0
final [0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]
outputput affine center: [-499.008, -499.005, 498.986]
output affine para: [0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]
initial measure value (MMI): rval = -2.39058
final measure value (MMI): rval = -2.38938
finish affine registeration...
 Requested Transformation Model:  SyN : Using 
SyN diffeomorphic model for transformation. 
  Grad Step 0.2 total-smoothing 1 gradient-smoothing 2
 setting N-TimeSteps = 1 trunc 256
 ScaleFactor 8 nlev 4 curl 0
 allocated def field -1 0 0
0 -1 0
0 0 1

 Its at this level 2
 Allocating 
 Allocating Done 
 iteration 1 energy 0 : 0
 iteration 2 energy 0 : 0
 tired convergence: reached max iterations 
 ScaleFactor 4 nlev 4 curl 1
 Its at this level 2
 iteration 1 energy 0 : 0
 iteration 2 energy 0 : 0
 tired convergence: reached max iterations 
 ScaleFactor 2 nlev 4 curl 2
 Its at this level 1
 iteration 1 energy 0 : 0
 tired convergence: reached max iterations 
 ScaleFactor 1 nlev 4 curl 3
 Its at this level 1
 iteration 1 energy 0 : 0
 tired convergence: reached max iterations 
 Registration Done 
 begin writing /tmp/ants.nii
 writing /tmp/ants affine 
 writing /tmp/ants def 
filename /tmp/antsWarp.nii

$ python3 -c "import nibabel as nib; print(nib.load('/tmp/antsWarp.nii').header)"
<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'
sizeof_hdr      : 348
data_type       : b''
db_name         : b''
extents         : 0
session_error   : 0
regular         : b'r'
dim_info        : 0
dim             : [  5 378 378 378   1   3   1   1]     # <- the array size has decreased from 500x500x500
intent_p1       : 0.0
intent_p2       : 0.0
intent_p3       : 0.0
intent_code     : vector
datatype        : float32
bitpix          : 32
slice_start     : 0
pixdim          : [1. 2. 2. 2. 0. 0. 0. 0.]     # <- while the voxel size is unchanged
vox_offset      : 0.0
scl_slope       : nan
scl_inter       : nan
slice_end       : 0
slice_code      : unknown
xyzt_units      : 2
cal_max         : 0.0
cal_min         : 0.0
slice_duration  : 0.0
toffset         : 0.0
glmax           : 0
glmin           : 0
descrip         : b''
aux_file        : b''
qform_code      : scanner
sform_code      : unknown
quatern_b       : 0.0
quatern_c       : 0.0
quatern_d       : 0.0
qoffset_x       : -0.0
qoffset_y       : -0.0
qoffset_z       : 0.0
srow_x          : [0. 0. 0. 0.]
srow_y          : [0. 0. 0. 0.]
srow_z          : [0. 0. 0. 0.]
intent_name     : b''
magic           : b'n+1'

This seems to be independent of IO (e.g., MINC vs NIFTI inputs/outputs produce the same result).

System information

OS: Ubuntu Linux
OS version: 16.04 (system gcc 5.4.0)
Type of system: desktop
OS: CentOS
OS version: 6 (gcc 7.2.0)
Type of system: cluster
OS: CentOS + Nixpkgs
OS version: 7 (gcc 8.3)
Type of system: cluster

ANTs version information

ANTs code version: [output of antsRegistration --version]
```
ANTs Version: 2.3.4.dev170-g11953
Compiled: Sep  3 2020 12:51:41
```
(also present with 2.3.1.dev90-g55fa5, but not in 2.2.0)
ANTs installation type: compiled from source using superbuild and nearly default options (except enabling shared libs)

Additional information

Running with a debug build produces no errors, but valgrind produces a couple suspicious things such as:

==28741== Conditional jump or move depends on uninitialised value(s)
==28741==    at 0x4D7FA66: SetEdgePaddingValue (itkWarpImageMultiTransformFilter.h:217)
==28741==    by 0x4D7FA66: itk::ANTSImageRegistrationOptimizer<3u, float>::WarpMultiTransform(itk::SmartPointer<itk::Image<float, 3u> >, itk::SmartPointer<itk::Image<float, 3u> >, itk::SmartPointer<itk::MatrixOffsetTransformBase<double, 3u, 3u> >, itk::SmartPointer<itk::Image<itk::Vector<float, 3u>, 3u> >, bool, itk::SmartPointer<itk::MatrixOffsetTransformBase<double, 3u, 3u> >) (itkANTSImageRegistrationOptimizer.h:644)
 ...

but I didn't investigate this further.

Initial discovery by @mcvaneede.

ntustison commented 4 years ago

That will happen if you set the number of iterations at the top levels to 0.

bcdarwin commented 4 years ago

Sorry, I disabled the finest resolution iterations just for the purposes of speed. The output shape is unchanged if you run some. (I've updated the code example and output above.)

(Also, even if true, wouldn't the output header above still imply the voxel spacing of the output had been set incorrectly?)

ntustison commented 4 years ago

500 is way too large for a debugging example. I can reproduce with an image of 100x100x100.

The pyramid scheme is not resampling properly for this particular number of levels and image size. Definitely a bug and something for which a fix and pull request would be nice. But it's probably not something anybody will be fixing on our end as ANTS is somewhat deprecated.

bcdarwin commented 4 years ago

Apparently more than somewhat: I didn't adjust any parameters above to trigger the bug, so this means that an arbitrary ANTS call is now likely broken.

stnava commented 4 years ago

As nick said, we don’t maintain that program , we leave it up to the community.

The old program is superseded by antsRegistration.

On Fri, Sep 4, 2020 at 6:41 PM Ben Darwin notifications@github.com wrote:

Apparently more than somewhat: I didn't adjust any parameters above to trigger the bug, so this means that an arbitrary ANTS call is now likely broken.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/ANTsX/ANTs/issues/1076#issuecomment-687440448, or unsubscribe https://github.com/notifications/unsubscribe-auth/AACPE7TYDJACKBMICS5M2PLSEFUJHANCNFSM4QZJNA6A .

--

brian

bcdarwin commented 4 years ago

Of course I don't mind that, I just have some code that I didn't switch to antsRegistration since I didn't know this was the case (it's not gated in the cmake config, doesn't print any warning, I didn't see a changelog, etc.)

Thanks for clarifying!

ntustison commented 4 years ago

One thought, though, is that this is an obvious bug---something that we would've caught before now. I'm wondering if something was checked in recently or if ITK made some changes which are now causing this issue. Perhaps I'll check if I can find some time.

gdevenyi commented 4 years ago

Git bisect would help here

cookpa commented 4 years ago

Bisecting now, will update when I find the problem.

cookpa commented 4 years ago

f0977f8c3c69fdb6e7b98231f85f9a491007041e is the first bad commit
commit f0977f8c3c69fdb6e7b98231f85f9a491007041e
Author: Hans Johnson <hans.j.johnson@gmail.com>
Date:   Sun Aug 18 14:18:04 2019 -0500

    COMP: Allow building with ITK_LEGACY_REMOVE turned on

    The VectorExpandImage filter specialization is no longer needed because
    the ExpandImage filter now works with both vector and scalar types
    uniformly.

 ImageRegistration/itkANTSImageRegistrationOptimizer.h | 6 +++---
 SuperBuild/External_ITKv5.cmake                       | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

https://github.com/ANTsX/ANTs/commit/f0977f8c3c69fdb6e7b98231f85f9a491007041e

ntustison commented 4 years ago

That's awesome, @cookpa and makes sense. Do you want me to fix it or do you want to give it a go?

cookpa commented 4 years ago

Would you be OK with putting the old vector expander class back in? That's all I can see that might explain the change.

ntustison commented 4 years ago

Yeah, that's the first thing I would try and, if it works, I would just leave the old one.

cookpa commented 4 years ago

@ntustison looks like that class has been removed from ITK

ntustison commented 4 years ago

Hey @cookpa , Sorry just got back from a run but I'm pretty sure I found the issue without even looking at the ANTs code. Look at the different signatures for the likely culprit function:

itkVectorExpandImageFilter.h

# Line 134
SetExpandFactors(const float factor);

vs.

itkExpandImageFilter.h

# Line 112
SetExpandFactors(const unsigned int factor);

Going to do some digging in the optimizer class to verify that this is the issue. Tagging @hjmjohnson since he made the original change and I don't know yet if this needs to be reported back to ITK.

cookpa commented 4 years ago

That function is overloaded, the other version takes a templated array. I'm going to try to get that to take a float

EDIT: it's a static attribute, so I don't think I can change it.

ntustison commented 4 years ago

Right. But am I not reading correctly that a related difference is the template type for that array (Line 106 vs. Lines 116-117)?

cookpa commented 4 years ago

No you're right. I was looking at the other method:

virtual void | SetExpandFactors (const unsigned int factor)

But you're right, it's declared as float in the vector class and unsigned int in the regular class. So unless there's some way to change that at compile time, it's not going to work for us.

ntustison commented 4 years ago

Okay, then I'm thinking we simply replace the ExpandImageFilter approach with appropriate calls to the ResampleImageFilter class. How does that sound?

ntustison commented 4 years ago

https://github.com/ANTsX/ANTs/pull/1078

ANTsX / ANTs

Resulting warp unexpectedly has a smaller field of view than input images #1076