[Insight-users] Debugging bspline registration

[Andriy Fedorov]

Hi,

I have a pair of 3D MRI images that I am registering using a custom
ITK registration code I developed. On one of the platforms I use, I
observe a very strange change in behavior of this code as a result of
a minor difference in the header of the input image. I've been trying
to understand this behavior for most of the day today ...

Specifically, this minor difference is the following (from .nrrd header):

Image1: "space origin:
(-79.285728454589872,-72.241912841796974,-38.198635101318366)"

vs

Image2: "space origin:
(-79.285728454589872,-72.241912841796989,-38.198635101318366)"

(the last two decimal digits of the Y coordinate of the origin)

If I choose Image1 as the reference image, I have a meaningful
registration result. If I choose Image2, I have the following painful
error message:

File: /home/prostate/fedorov/Slicer-3.6-IGT/Slicer3-lib/Insight/Code/Review/itkOptMattesMutualInformationImageToImageMetric.txx
Line: 1077
Description: itk::ERROR:
MattesMutualInformationImageToImageMetric(0x7f1b0c3dbcc0): Too many
samples map outside moving image buffer: 634 / 100000

All of the parameters of registration are otherwise identical (I run
both executions from the same command line just changing the reference
image). Some more technical details:

 * The random number generator is initialized with the same seed
 * same number of threads (16) is used
 * my code uses multi-resolution registration approach going rigid ->
rigid+scaling -> rigid+anisotropic scaling -> affine -> bspline
transformation
 * I confirmed that the affine transformation that I set as bulk for
bspline gives good initial alignment.
 * I use masks for both images, and these masks overlap very well
after affine registration.
 * The bspline grid is coarse, 4x4x4 knots.
 * the metric is MattesMutualInformation
 * both images have non-identity direction cosines
 * code compiled in Release mode, gcc 4.4.1, Ubuntu 9.04

When I look at the metric value during bspline registration, I see it
starts from the same values in both cases, but then diverges starting
at the third iteration:

Image1:

0   -0.036445
0   -0.0366935
0   -0.0366935
1   -0.0372417
1   -0.0372417
2   -0.0368946
2   -0.0373551
2   -0.0373551
3   -0.030338
3   -0.0373644
3   -0.0373644
4   -0.0319827
4   -0.037446
4   -0.037446
...
10   -0.0377343
10   -0.0383924
10   -0.0384768
10   -0.0384857
10   -0.0384966
10   -0.0384927
10   -0.0384966
10   -0.0384966
10   -0.0384966
10   -0.0384966

Image2:

0   -0.036445
0   -0.0366935
0   -0.0366935
1   -0.0372417
1   -0.0372417
2   -0.0368945
2   -0.0373552
2   -0.0373552
3   -0.030338
3   -0.0373643
3   -0.0373643
4   -0.0319829
4   -0.0374433
4   -0.0374433
...
12   -0.0363308
12   -0.0373945
12   -0.037511
12   -0.0375601
12   -0.0375916
12   -0.0375924
12   -0.0375819
12   -0.0375822
12   -0.0375924
12   -0.0375924
12   -0.0375926
12   -0.0375926
12   -0.0375926

When I compile the same code in Debug mode, bspline optimizer reaches
convergence without error independently of the image choice. When I
run the same example with 1 thread, I also have convergence and no
error.

Interestingly, the number of iterations with 1 thread is different for
each of the transformation levels. BSpline converges in 6 iterations.
Also, execution time is very similar for 1 thread as it is for 16
threads (just ~1 minute!), with the registration results very similar
visually.

The reason I came up with this degenerate case is that Image2 was read
and saved by a separate image visualization tool (3D Slicer) that uses
my code as a plugin. That execution initially failed with the error
message above. After a lengthy investigation I finally narrowed the
difference down to this discrepancy in the image header.

How should I debug this problem? Is there a known open bug in ITK that
would explain this and/or inconsistency of the optimizers/metric
behavior depending on the number of threads?

I tried 3.18 and 3.20 branches of ITK. I am running valgrind now to
find any possible leaks or data corruption, but this will take a
while.

I would appreciate any hints.

--
Andriy Fedorov, Ph.D.

Research Fellow
Brigham and Women's Hospital
Harvard Medical School
75 Francis Street
Boston, MA 02115 USA
fedorov at bwh.harvard.edu
(617) 525-6258 (office)