[Insight-developers] Bug in vnl_powell.cxx

[Hans Johnson]

Hello VXL Developers,

I was producing inconsistent results in my medical image processing software
when using powell optimization in my code that prompted me to run valgrind
(See report at end of this message).

The least invasive way to avoid this bug is to initialize the variable
³bx=0² on lines 87 and 124.   This, however, is not the best solution in my
opinion.  Currently the the powell optimizer uses vn_brent.h which states in
it¹s header 
==============================
//: Brent 1D minimizer (deprecated)
//
// Please use vnl_brent_minimizer instead.
//
// This routine used to contain copyrighted code, and is deprecated.
// It is now simply a wrapper around vnl_brent_minimizer.
==============================

I¹ve made a minimally modified version of vnl_powell.cxx that uses
vnl_brent_minimizer directly (saving 3 levels of indirect function calls)
that makes it much easier to read and removes the errors.

The bug fixed files can be found at:
http://www.nitrc.org/plugins/scmsvn/viewcvs.php/trunk/lib/vnl_powell_fixed.c
xx?revision=172&root=art&sortdir=down&view=markup
Or as part of the ITK bug report at
http://public.kitware.com/Bug/view.php?id=7903.

In my program, this fix reliably runs to completion in 24 seconds where it
had been occasionally segmentation faulting or running somewhere between 120
seconds to 300 seconds depending on what value was being filled in.

Please consider implementing this bug fix in the next release of vxl.

Thanks,
Hans J. Johnson
hans-johnson at uiowa.edu

## Valgrind is our friend!

  ==15048== Conditional jump or move depends on uninitialised value(s)
  ==15048==    at 0x890BE9: vnl_bracket_minimum(vnl_cost_function&, double&,
double&, double&, double&, double&, double&) (vnl_bracket_minimum.cxx:46)
  ==15048==    by 0x890938: vnl_brent::bracket_minimum(double*, double*,
double*, double*, double*, double*) (vnl_brent.cxx:58)
  ==15048==    by 0x8908E0: vnl_brent::bracket_minimum(double*, double*,
double*) (vnl_brent.cxx:52)
  ==15048==    by 0x88F876: vnl_powell::minimize(vnl_vector<double>&)
(vnl_powell.cxx:85)
  ==15048==    by 0x676E4D:
PowellOptimizeMSP(itk::SmartPointer<itk::OrientedImage<short, 3> >&,
itk::Point<double, 3> const&, itk::Point<double, 3> const&,
itk::Point<double, 3> const&, PlaneObject const&, itk::Point<double, 3>&,
Plan
eObject&, double) (acpcdetect_itk.cxx:361)
  ==15048==    by 0x678079:
ComputeMSP(itk::SmartPointer<itk::OrientedImage<short, 3> >, PlaneObject&,
itk::Point<double, 3>&, itk::Matrix<double, 3, 3>&, unsigned)
(acpcdetect_itk.cxx:679)
  ==15048==    by 0x678571:
computeTmsp(itk::SmartPointer<itk::OrientedImage<short, 3> >&, PlaneObject&,
itk::Point<double, 3>&) (acpcdetect_itk.cxx:739)
  ==15048==    by 0x56D4C8: main (acpcResampleMSP.cxx:74)



vnl_powell.cxx
====================
lines 82-85 && lines 122-125
       double ax = 0.0;
       double xx = initial_step_;
       double bx;                             /**NOTE bx is not
initialized!*/
       brent.bracket_minimum(&ax, &xx, &bx);

vnl_brent.cxx
====================
lines   49-59
 void vnl_brent::bracket_minimum(double *ax, double *bx, double *cx)
//**NOTE cx is not initialized
 {
   double fa, fb, fc;
   bracket_minimum(ax,bx,cx,&fa,&fb,&fc);
 }

 void vnl_brent::bracket_minimum(double *ax, double *bx, double *cx,
//**NOTE cx is not initialized
                                 double *fa, double *fb, double *fc)
 {
   vnl_bracket_minimum( *f_, *cx, *bx, *ax, *fc, *fb, *fa );
                            /*^----------- NOTE: cx is not initialized */
 }

vnl_bracket_minimum.cxx
====================
lines 38-48

 void vnl_bracket_minimum(vnl_cost_function& fn,
                         double& a, double& b, double& c, //** Note a is not
initialized
                         double& fa, double& fb, double& fc)
 {
   // Set up object to evaluate function
   // Note that fn takes a vector input - f converts a scalar to a vector
   vnl_bm_func f(fn);

   if (b==a) b=a+1.0;
        /*^-------------- Note comparison of value where a is not
initialized.  In my case, a=324e+124 sometimes
   fa = f(a);
        /*^-------A not initialized here, and should not be used (this is bx
in the original calling funciton.*/
   fb = f(b);



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20081030/49609765/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: vnl_powell_fixed.cxx
Type: application/octet-stream
Size: 4378 bytes
Desc: not available
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20081030/49609765/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: vnl_powell_fixed.h
Type: application/octet-stream
Size: 1458 bytes
Desc: not available
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20081030/49609765/attachment-0001.obj>