From a7b0a734cf79ac82e70a61321b41ccfee4d5feec Mon Sep 17 00:00:00 2001 From: Stephen Aylward Date: Wed, 23 Apr 2008 08:49:11 -0400 Subject: [PATCH] ENH: Eliminates variable construct/destruct - important for large parameters spaces, e.g., during bspline transform optimization --- Code/Numerics/itkFRPROptimizer.cxx | 38 ++++++++------ Code/Numerics/itkFRPROptimizer.h | 3 ++ Code/Numerics/itkOptimizer.h | 6 ++- Code/Numerics/itkPowellOptimizer.cxx | 77 +++++++++++++++++++--------- Code/Numerics/itkPowellOptimizer.h | 20 +++++++- 5 files changed, 103 insertions(+), 41 deletions(-) diff --git a/Code/Numerics/itkFRPROptimizer.cxx b/Code/Numerics/itkFRPROptimizer.cxx index 13c6c2b310a..fa6b0e7ccb5 100644 --- a/Code/Numerics/itkFRPROptimizer.cxx +++ b/Code/Numerics/itkFRPROptimizer.cxx @@ -69,6 +69,15 @@ ::GetValueAndDerivative(ParametersType & p, double * val, void FRPROptimizer ::LineOptimize(ParametersType * p, ParametersType & xi, double * val) +{ + ParametersType tempCoord( this->GetSpaceDimension() ); + this->LineOptimize( p, xi, val, tempCoord ); +} + +void +FRPROptimizer +::LineOptimize(ParametersType * p, ParametersType & xi, double * val, + ParametersType & tempCoord ) { this->SetLine(*p, xi); @@ -79,15 +88,14 @@ ::LineOptimize(ParametersType * p, ParametersType & xi, double * val) double bx; double fb; - ParametersType pp = (*p); - - this->LineBracket(&ax, &xx, &bx, &fa, &fx, &fb); + this->LineBracket(&ax, &xx, &bx, &fa, &fx, &fb, tempCoord); this->SetCurrentLinePoint(xx, fx); double extX = 0; double extVal = 0; - this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &extX, &extVal); + this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &extX, &extVal, + tempCoord); this->SetCurrentLinePoint(extX, extVal); (*p) = this->GetCurrentPosition(); @@ -110,21 +118,21 @@ ::StartOptimization() this->SetSpaceDimension(m_CostFunction->GetNumberOfParameters()); - const unsigned int SpaceDimension = this->GetSpaceDimension(); + FRPROptimizer::ParametersType tempCoord( this->GetSpaceDimension() ); double gg, gam, dgg; - FRPROptimizer::ParametersType g( SpaceDimension ); - FRPROptimizer::ParametersType h( SpaceDimension ); - FRPROptimizer::ParametersType xi( SpaceDimension ); + FRPROptimizer::ParametersType g( this->GetSpaceDimension() ); + FRPROptimizer::ParametersType h( this->GetSpaceDimension() ); + FRPROptimizer::ParametersType xi( this->GetSpaceDimension() ); - FRPROptimizer::ParametersType p( SpaceDimension ); + FRPROptimizer::ParametersType p( this->GetSpaceDimension() ); p = this->GetInitialPosition(); this->SetCurrentPosition(p); double fp; this->GetValueAndDerivative(p, &fp, &xi); - for( i = 0; i < SpaceDimension; i++ ) + for( i = 0; i < this->GetSpaceDimension(); i++ ) { g[i] = -xi[i]; xi[i] = g[i]; @@ -141,12 +149,12 @@ ::StartOptimization() double fret; fret = fp; - this->LineOptimize(&p, xi, &fret); + this->LineOptimize(&p, xi, &fret, tempCoord); if ( 2.0 * vcl_abs(fret - fp) <= this->GetValueTolerance() * (vcl_abs(fret)+ vcl_abs(fp) + FRPR_TINY) ) { - if( limitCount < SpaceDimension ) + if( limitCount < this->GetSpaceDimension() ) { this->GetValueAndDerivative(p, &fp, &xi); xi[limitCount] = 1; @@ -170,7 +178,7 @@ ::StartOptimization() if( m_OptimizationType == PolakRibiere ) { - for( i=0; i< SpaceDimension; i++ ) + for( i=0; i< this->GetSpaceDimension(); i++ ) { gg += g[i] * g[i]; dgg += (xi[i] + g[i]) * xi[i]; @@ -178,7 +186,7 @@ ::StartOptimization() } if( m_OptimizationType == FletchReeves ) { - for( i=0; i< SpaceDimension; i++ ) + for( i=0; i< this->GetSpaceDimension(); i++ ) { gg += g[i] * g[i]; dgg += xi[i] * xi[i]; @@ -193,7 +201,7 @@ ::StartOptimization() } gam = dgg/gg; - for( i = 0; i < SpaceDimension; i++) + for( i = 0; i < this->GetSpaceDimension(); i++) { g[i] = -xi[i]; xi[i] = g[i] + gam * h[i]; diff --git a/Code/Numerics/itkFRPROptimizer.h b/Code/Numerics/itkFRPROptimizer.h index 8dba4268973..8f95b730767 100644 --- a/Code/Numerics/itkFRPROptimizer.h +++ b/Code/Numerics/itkFRPROptimizer.h @@ -96,6 +96,9 @@ class ITK_EXPORT FRPROptimizer: virtual void LineOptimize(ParametersType * p, ParametersType & xi, double * val ); + virtual void LineOptimize(ParametersType * p, ParametersType & xi, + double * val, + ParametersType & tempCoord ); private: diff --git a/Code/Numerics/itkOptimizer.h b/Code/Numerics/itkOptimizer.h index d7b0835bdd0..0d721c90cd3 100644 --- a/Code/Numerics/itkOptimizer.h +++ b/Code/Numerics/itkOptimizer.h @@ -89,12 +89,16 @@ class ITK_EXPORT Optimizer : public Object bool m_ScalesInitialized; + // Keep m_CurrentPosition as a protected var so that subclasses can + // have fast access. This is important when optimizing high-dimensional + // spaces, e.g. bspline transforms. + ParametersType m_CurrentPosition; + private: Optimizer(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented ParametersType m_InitialPosition; - ParametersType m_CurrentPosition; ScalesType m_Scales; }; diff --git a/Code/Numerics/itkPowellOptimizer.cxx b/Code/Numerics/itkPowellOptimizer.cxx index 6fb38de0de7..1b6b5ce3a58 100644 --- a/Code/Numerics/itkPowellOptimizer.cxx +++ b/Code/Numerics/itkPowellOptimizer.cxx @@ -55,11 +55,10 @@ PowellOptimizer ::SetLine(const PowellOptimizer::ParametersType & origin, const vnl_vector & direction) { - m_LineOrigin = origin; - m_LineDirection = direction; for(unsigned int i=0; iGetScales()[i]; + m_LineOrigin[i] = origin[i]; + m_LineDirection[i] = direction[i] / this->GetScales()[i]; } } @@ -67,18 +66,25 @@ double PowellOptimizer ::GetLineValue(double x) const { - PowellOptimizer::ParametersType xCoord( m_SpaceDimension ); + PowellOptimizer::ParametersType tempCoord( m_SpaceDimension ); + return this->GetLineValue(x, tempCoord); +} + +double +PowellOptimizer +::GetLineValue(double x, ParametersType & tempCoord) const +{ for(unsigned int i=0; im_LineOrigin[i] + x * this->m_LineDirection[i]; + tempCoord[i] = this->m_LineOrigin[i] + x * this->m_LineDirection[i]; } if(m_Maximize) { - return -(this->m_CostFunction->GetValue(xCoord)); + return -(this->m_CostFunction->GetValue(tempCoord)); } else { - return this->m_CostFunction->GetValue(xCoord); + return this->m_CostFunction->GetValue(tempCoord); } } @@ -86,12 +92,10 @@ void PowellOptimizer ::SetCurrentLinePoint(double x, double fx) { - PowellOptimizer::ParametersType xCoord( m_SpaceDimension ); for(unsigned int i=0; im_LineOrigin[i] + x * this->m_LineDirection[i]; + this->m_CurrentPosition[i] = this->m_LineOrigin[i] + x * this->m_LineDirection[i]; } - this->SetCurrentPosition(xCoord); if(m_Maximize) { this->SetCurrentCost(-fx); @@ -100,6 +104,7 @@ ::SetCurrentLinePoint(double x, double fx) { this->SetCurrentCost(fx); } + this->Modified(); } void @@ -143,6 +148,16 @@ void PowellOptimizer ::LineBracket(double * x1, double * x2, double * x3, double * f1, double * f2, double * f3) +{ + PowellOptimizer::ParametersType tempCoord( m_SpaceDimension ); + this->LineBracket( x1, x2, x3, f1, f2, f3, tempCoord); +} + +void +PowellOptimizer +::LineBracket(double * x1, double * x2, double * x3, + double * f1, double * f2, double * f3, + ParametersType & tempCoord) { // // Compute the golden ratio as a constant to be @@ -153,7 +168,7 @@ ::LineBracket(double * x1, double * x2, double * x3, // // Get the value of the function for point x2 // - *f2 = this->GetLineValue( *x2 ); + *f2 = this->GetLineValue( *x2, tempCoord ); // // Compute extrapolated point using the golden ratio @@ -162,7 +177,7 @@ ::LineBracket(double * x1, double * x2, double * x3, { // compute x3 on the side of x2 *x3 = *x1 + goldenRatio * ( *x2 - *x1 ); - *f3 = this->GetLineValue( *x3 ); + *f3 = this->GetLineValue( *x3, tempCoord ); // If the new point is a minimum // then continue extrapolating in @@ -174,14 +189,14 @@ ::LineBracket(double * x1, double * x2, double * x3, *x2 = *x3; *f2 = *f3; *x3 = *x1 + goldenRatio * ( *x2 - *x1 ); - *f3 = this->GetLineValue( *x3 ); + *f3 = this->GetLineValue( *x3, tempCoord ); } } else { // compute x3 on the side of x1 *x3 = *x2 + goldenRatio * ( *x1 - *x2 ); - *f3 = this->GetLineValue( *x3 ); + *f3 = this->GetLineValue( *x3, tempCoord ); // If the new point is a minimum // then continue extrapolating in @@ -193,7 +208,7 @@ ::LineBracket(double * x1, double * x2, double * x3, *x1 = *x3; *f1 = *f3; *x3 = *x2 + goldenRatio * ( *x1 - *x2 ); - *f3 = this->GetLineValue( *x3 ); + *f3 = this->GetLineValue( *x3, tempCoord ); } } @@ -226,7 +241,7 @@ ::LineBracket(double * x1, double * x2, double * x3, // is no need for using it as a factor. // *x4 = *x1 - *x2 + *x3; - *f4 = this->GetLineValue( *x4 ); + *f4 = this->GetLineValue( *x4, tempCoord ); // // If the new value f4 at x4 is larger than f2 @@ -280,8 +295,19 @@ ::LineBracket(double * x1, double * x2, double * x3, void PowellOptimizer ::BracketedLineOptimize(double ax, double bx, double cx, - double itkNotUsed(fa), double functionValueOfb, double itkNotUsed(fc), + double fa, double functionValueOfb, double fc, double * extX, double * extVal) +{ + PowellOptimizer::ParametersType tempCoord( m_SpaceDimension ); + this->BracketedLineOptimize( ax, bx, cx, fa, functionValueOfb, fc, extX, extVal, tempCoord); +} + +void +PowellOptimizer +::BracketedLineOptimize(double ax, double bx, double cx, + double itkNotUsed(fa), double functionValueOfb, double itkNotUsed(fc), + double * extX, double * extVal, + ParametersType & tempCoord) { double x; double v = 0.0; @@ -389,7 +415,7 @@ ::BracketedLineOptimize(double ax, double bx, double cx, double functionValueOft; - functionValueOft = this->GetLineValue(t); + functionValueOft = this->GetLineValue(t, tempCoord); if( functionValueOft <= functionValueOfX ) { @@ -459,6 +485,7 @@ ::StartOptimization() m_SpaceDimension = m_CostFunction->GetNumberOfParameters(); m_LineOrigin.set_size(m_SpaceDimension); m_LineDirection.set_size(m_SpaceDimension); + this->m_CurrentPosition.set_size(m_SpaceDimension); vnl_matrix xi(m_SpaceDimension, m_SpaceDimension); vnl_vector xit(m_SpaceDimension); @@ -466,6 +493,8 @@ ::StartOptimization() xit.fill(0); xit[0] = 1; + PowellOptimizer::ParametersType tempCoord(m_SpaceDimension); + PowellOptimizer::ParametersType p(m_SpaceDimension); PowellOptimizer::ParametersType pt(m_SpaceDimension); PowellOptimizer::ParametersType ptt(m_SpaceDimension); @@ -479,7 +508,7 @@ ::StartOptimization() xx = 0; this->SetLine(p, xit); - fx = this->GetLineValue(0); + fx = this->GetLineValue(0, tempCoord); for (m_CurrentIteration = 0; m_CurrentIteration <= m_MaximumIteration; @@ -502,8 +531,8 @@ ::StartOptimization() ax = 0.0; fa = fx; xx = m_StepLength; - this->LineBracket( &ax, &xx, &bx, &fa, &fx, &fb); - this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &xx, &fx); + this->LineBracket( &ax, &xx, &bx, &fa, &fx, &fb, tempCoord); + this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &xx, &fx, tempCoord); this->SetCurrentLinePoint(xx, fx); p = this->GetCurrentPosition(); @@ -529,7 +558,7 @@ ::StartOptimization() } this->SetLine(ptt, xit); - fptt = this->GetLineValue(0); + fptt = this->GetLineValue(0, tempCoord); if (fptt < fp) { double t = 2.0 * (fp - 2.0*fx + fptt) @@ -542,8 +571,8 @@ ::StartOptimization() ax = 0.0; fa = fx; xx = 1; - this->LineBracket( &ax, &xx, &bx, &fa, &fx, &fb); - this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &xx, &fx); + this->LineBracket( &ax, &xx, &bx, &fa, &fx, &fb, tempCoord); + this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &xx, &fx, tempCoord); this->SetCurrentLinePoint(xx, fx); p = this->GetCurrentPosition(); diff --git a/Code/Numerics/itkPowellOptimizer.h b/Code/Numerics/itkPowellOptimizer.h index bc6c8b1ee7a..f5c3aa3abec 100644 --- a/Code/Numerics/itkPowellOptimizer.h +++ b/Code/Numerics/itkPowellOptimizer.h @@ -145,6 +145,8 @@ class ITK_EXPORT PowellOptimizer: * Line origin and distances are set via SetLine */ double GetLineValue(double x) const; + double GetLineValue(double x, ParametersType & tempCoord) const; + /** Set the given scalar step distance (x) and function value (fx) as the * "best-so-far" optimizer values. */ void SetCurrentLinePoint(double x, double fx); @@ -169,6 +171,10 @@ class ITK_EXPORT PowellOptimizer: virtual void LineBracket(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc); + virtual void LineBracket(double *ax, double *bx, double *cx, + double *fa, double *fb, double *fc, + ParametersType & tempCoord); + /** Given a bracketing triple of points and their function values, returns * a bounded extreme. These values are in parameter space, along the * current line and wrt the current origin set via SetLine. Optimization @@ -178,8 +184,20 @@ class ITK_EXPORT PowellOptimizer: double fa, double fb, double fc, double * extX, double * extVal); + virtual void BracketedLineOptimize(double ax, double bx, double cx, + double fa, double fb, double fc, + double * extX, double * extVal, + ParametersType & tempCoord); + itkGetMacro(SpaceDimension, unsigned int); - itkSetMacro(SpaceDimension, unsigned int); + void SetSpaceDimension( unsigned int dim ) + { + this->m_SpaceDimension = dim; + this->m_LineDirection.set_size( dim ); + this->m_LineOrigin.set_size( dim ); + this->m_CurrentPosition.set_size( dim ); + this->Modified(); + } itkSetMacro(CurrentIteration, unsigned int);