VXL Developers,
Sorry for the repeated mail message. The previous message was sent as HTML
and was not readable. I hope that this plain text message is easier to deal
with.
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).
I just found a nasty bug in the vxl_powell.cxx optimizer. I¹ve done
extensive testing on it, and figured out how fix it. This has given my
application great speed increase, and no longer segmentation faults
randomly. The problem was that a variable was not explicitly being set, and
from time to time it would have a non-zero initial value. If you are really
lucky then the default initialization is zero, and all is well; if you are
somewhat lucky, the random value does not affect the outcome too much, and
the method converges after a bunch of iterations, and occasionally you are
unlucky and the the initialization is too crazy big/small and the function
being minimized evaluates to nan or inf and all heck breaks loose.
The least invasive way to avoid this bug is to initialize the variable
³bx=0² on lines 87 and 124 of vnl_powell.cxx. This does work, but a little
more work also removes two façade functions that adapt the deprecated
vnl_brent class member functions into the vnl_brent_minimizer class. 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 2 levels of indirect function calls)
that makes it much easier to read and removes the errors and dependence on
the deprecated class.
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@...
## 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);
--
Hans J. Johnson, Ph.D.
Hans-johnson@...
278 GH
The University of Iowa
Iowa City, IA 52241
(319) 353 8587
Notice: This UI Health Care e-mail (including attachments) is covered by the Electronic Communications Privacy Act, 18 U.S.C. 2510-2521, is confidential and may be legally privileged. If you are not the intended recipient, you are hereby notified that any retention, dissemination, distribution, or copying of this communication is strictly prohibited. Please reply to the sender that you have received the message in error, then delete it. Thank you.
|