I was producing inconsistent results in my medical image processing softwar= e 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 giv= en my application great speed increase, and no longer segmentation faults ra= ndomly. 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 a= re 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 mu= ch, and the method converges after a bunch of iterations, and occasionally y= ou are unlucky and the the initialization is too crazy big/small and the fun= ction being minmized evaluates to nan or inf and all heck breaks loose.

The least invasive way to avoid this bug is to initialize the variable R= 20;bx=3D0” 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 t= he deprecated vnl_brent class member functions into the vnl_brent_minimizer = class. This, however, is not the best solution in my opinion. &n= bsp;Currently the the powell optimizer uses vn_brent.h which states in it= 217;s header

=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D

//: 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.

=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D

I’ve made a minimally modified version of vnl_powell.cxx that uses vn= l_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 d= eprecated class.

The bug fixed files can be found at:

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

## Valgrind is our friend!

=3D=3D15048=3D=3D Conditional jump or move depend= s on uninitialised value(s)

=3D=3D15048=3D=3D at 0x890BE9: vnl_bracket_minimum(= vnl_cost_function&, double&, double&, double&, double&, = double&, double&) (vnl_bracket_minimum.cxx:46)

=3D=3D15048=3D=3D by 0x890938: vnl_brent::bracket_min= imum(double*, double*, double*, double*, double*, double*) (vnl_brent.cxx:58= )

=3D=3D15048=3D=3D by 0x8908E0: vnl_brent::bracket_min= imum(double*, double*, double*) (vnl_brent.cxx:52)

=3D=3D15048=3D=3D by 0x88F876: vnl_powell::minimize(v= nl_vector<double>&) (vnl_powell.cxx:85)

=3D=3D15048=3D=3D by 0x676E4D: PowellOptimizeMSP(itk:= :SmartPointer<itk::OrientedImage<short, 3> >&, itk::Point<= ;double, 3> const&, itk::Point<double, 3> const&, itk::Poin= t<double, 3> const&, PlaneObject const&, itk::Point<double,= 3>&, Plan

eObject&, double) (acpcdetect_itk.cxx:361)

=3D=3D15048=3D=3D by 0x678079: ComputeMSP(itk::SmartP= ointer<itk::OrientedImage<short, 3> >, PlaneObject&, itk::Po= int<double, 3>&, itk::Matrix<double, 3, 3>&, unsigned) (= acpcdetect_itk.cxx:679)

=3D=3D15048=3D=3D by 0x678571: computeTmsp(itk::Smart= Pointer<itk::OrientedImage<short, 3> >&, PlaneObject&, i= tk::Point<double, 3>&) (acpcdetect_itk.cxx:739)

=3D=3D15048=3D=3D by 0x56D4C8: main (acpcResampleMSP.= cxx:74)

vnl_powell.cxx

=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D

lines 82-85 && lines 122-125

double ax =3D 0.0;

double xx =3D initial_step_;

double bx; &nb= sp; &= nbsp; /**NOTE bx is not initialized!*/

brent.bracket_minimum(&ax, = &xx, &bx);

vnl_brent.cxx

=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D

lines 49-59

void vnl_brent::bracket_minimum(double *ax, double *bx, double *cx) //**NOTE cx is not initial= ized

{

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 initia= lized

&= nbsp; = ; double *fa, double *fb, dou= ble *fc)

{

vnl_bracket_minimum( *f_, *cx, *bx, *ax, *fc, *fb, *fa );

&nb= sp; <= FONT COLOR=3D"#FE0000"> /*^----------- NOTE: cx is not initialized= */

}

vnl_bracket_minimum.cxx

=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D

lines 38-48

void vnl_bracket_minimum(vnl_cost_function& fn,

&nb= sp; d= ouble& a, double& b, double& c, //*= * Note a is not initialized

&nb= sp; d= ouble& fa, double& fb, double& fc)

{

// Set up object to evaluate function

// Note that fn takes a vector input - f converts a scala= r to a vector

vnl_bm_func f(fn);

if (b=3D=3Da) b=3Da+1.0;

/*^-------= ------- Note comparison of value where a is not initialized. In my cas= e, a=3D324e+124 sometimes

fa =3D f(a);

/*^--= -----A not initialized here, and should not be used (this is bx in the origi= nal calling funciton.*/

fb =3D f(b);

--B_3308241462_61765521--