[Insight-developers] Problems with the vnl library

[Paul Hughett]

Here is the current status of the problems I claimed in vnl, following
some offline discussions including one of the vnl maintainers.  Note
particularly number 2, in which I recommend a change to the way we
are writing the Insight code.


> 1. vnl_matrix_inverse does not check for singular matrices.

This is indeed a problem with vnl.  They don't appear to have any
current plans to fix it, though I have the impression they would
be more than happy to accept someone else's fix.

Quoting from one of the VXL developers, after I commented that an
exception or error flag would be better than printing out an error
message for singular matrices:

"In my opinion VXL should start thinking about exceptions soon but at the
moment it is one of the "new" C++ features that we are ignoring because
some users still use old compilers (removing support for old compilers can
be hard to push though). If more users say they want exceptions it will
happen sooner, though."


> 2. The eigenvalues computed by vnl_symmetric_eigensystem are not
very accurate.

I misdiagnosed this one.  It's not a vnl problem but an Insight problem.
The principal moments (which are eigenvalues of the central moments
matrix) computed by itkImageMomentsTest were much less accurate than
they should have been, considering that I was doing all the computations
in double.  I finally traced the problem to the fact that the Origin
and Scaling attributes of the Image class are stored as float rather
than double; this imprecision was contaminating the moments matrix
when it was converted to physical coordinates, and preventing the
eigenvalues from being computed to the expected accuracy.

I would recommend that the Origin and Spacing of the Image class be
converted from float to double.  The space penalty is negligible
compared to the space taken up by the image data, and the time
penalty is negligible, if not negative, on modern workstations.
Can anyone see a compelling reason not to convert these to double?

More generally, I would recommend that we make double rather than
float the standard within Insight, except where there is a compelling
reason for float.  Even there, I would make double the default and
float an option available to the user.  Examples:  Voxel data should
be either float or double, at the user's option.  Point coordinates
(e.g. the Point class) should be either double only, or double/float
at the user's option.  I would argue double only, but there is a case
to be made for allowing the option.  For things like origin/spacing,
moments, etc, I can't see any advantage to using anything other than
double; there's not enough data there that the space savings are worth
anything.


> 3. The diagonalizing transformation computed by the function
vnl_symmetric_eigensystem is improper.

This is also a problem of vnl.  The FORTRAN subroutine that they are
using internally to do the actual eigensystem computation does not
make the necessary information available.  As with the first problem,
they would be more than happy to accept contributions...


Paul Hughett