[MINC-users] failures when running N3withsmallsmoothingdistances

[Claude LEPAGE]

Andrew and Vitali,

> I took a bit of a look at this the other day and yes found that there
> is something strange going on with smaller step sizes.  Unbeknown to
> me Claude also was looking at N3 and not only found this bug too but
> fixed it.

I believe this is a different bug unrelated to the one I found. As I
explained to Oliver yesterday, let me explain it here.

N3 is solving a linear system for the coefficients of the splines. At
a small spline distance (relative to head size), there are more knots,
thus more degrees of freedom to solve. More precisely, for a cubic 
volume, let's say there are M*M*M knots (in 3-D). The matrix size is 
thus M*M*M by M*M*M for M^6 coefficients. Every time you reduce your 
distance by a factor of 2, you have 2^6 = 64 times more coefficients!
At a point, the system chokes - no more memory, swapping to disk and
very slow - or possibly crash. There is a fix, but it involves quite
a lot of work. It consists of using sparse matrix technology and 
iterative solvers to save space and time. However, this would be a
challenging programming task.

Another thing I noted that could explain the crash: valgrind (a 
memory checking tool) reports memory errors:

==22162== Invalid read of size 4
==22162==    at 0x818A3F7: TBSplineVolume::addDataPoint(int, int, int, double) (TBSpline.cc:1123)
==22162==    by 0x804BFC0: fitSplinesToVolumeLookup(TBSplineVolume*, volume_struct*, volume_struct*, Mat<double> const&, int) (splineSmooth.cc:477)
==22162==    by 0x804D028: main (splineSmooth.cc:154)
==22162==  Address 0x44D3104 is 4 bytes before a block of size 256 alloc'd
==22162==    at 0x401D14A: operator new[](unsigned) (vg_replace_malloc.c:195)
==22162==    by 0x818F9A3: TBSpline::TBSpline(Mat<double> const&, double, double, int) (TBSpline.cc:164)
==22162==    by 0x818FF2B: TBSplineVolume::TBSplineVolume(Mat<double> const&, double const*, double const*, int const*, double, double, int) (TBSpline.cc:887)
==22162==    by 0x804CCAB: main (splineSmooth.cc:144)

This is all very cryptic, but line 1123 of TBSpline.cc is a write
to the matrix. It looks like we are writing outside of the allocated
memory for the matrix when calculating the coefficients of the
matrix. I'm looking at this now.

As for the *other* bug I'm in the process of fixing, it will extend
the range of use of N3. In its current form, N3 is fine on human
brains at 1mm^3 sampling in stereotaxic MNI space. Beyond that, it
has difficulties in native space at different voxel sampling. For
example, try running N3 on a 1mm^3 sampling and repeat the experiment
after resampling your volume at 0.5mm^3. You'll see that N3 is not
invariant to voxel sampling. However, as I said, it is correct at
1mm^3 voxels. So no panic, please.

The purpose of my message is not to scare anyone, but to provide 
information on the problem. I hope you now better understand the
limitations of N3 when the spline distance is too small.

Yours,

Claude