itkEMLocalOptimization.h

Go to the documentation of this file.
00001 /*=auto=========================================================================
00002 
00003 (c) Copyright 2001 Massachusetts Institute of Technology 
00004 
00005 Permission is hereby granted, without payment, to copy, modify, display 
00006 and distribute this software and its documentation, if any, for any purpose, 
00007 provided that the above copyright notice and the following three paragraphs 
00008 appear on all copies of this software.  Use of this software constitutes 
00009 acceptance of these terms and conditions.
00010 
00011 IN NO EVENT SHALL MIT BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, 
00012 INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF THIS SOFTWARE 
00013 AND ITS DOCUMENTATION, EVEN IF MIT HAS BEEN ADVISED OF THE POSSIBILITY OF 
00014 SUCH DAMAGE.
00015 
00016 MIT SPECIFICALLY DISCLAIMS ANY EXPRESS OR IMPLIED WARRANTIES INCLUDING, 
00017 BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR 
00018 A PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
00019 
00020 THE SOFTWARE IS PROVIDED "AS IS."  MIT HAS NO OBLIGATION TO PROVIDE 
00021 MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
00022 
00023 =========================================================================auto=*/
00024 
00025 #ifndef _ITKEMLOCALREOPTIMIZATION_H_INCLUDED
00026 #define _ITKEMLOCALREOPTIMIZATION_H_INCLUDED 1
00027 
00028 #include "vtkEMSegment.h"
00029 #include "itkSingleValuedCostFunction.h"
00030 #include "EMLocalRegistrationCostFunction.h"
00031 #include "EMLocalShapeCostFunction.h"
00032 #include "itkPowellOptimizer.h"
00033 
00034 namespace itk
00035 {
00050 class VTK_EMSEGMENT_EXPORT  EMLocalCostFunctionWrapper :  public SingleValuedCostFunction
00051 {
00052 public:
00053 
00055   typedef EMLocalCostFunctionWrapper    Self;
00056   typedef SingleValuedCostFunction  Superclass;
00057 
00058   typedef SmartPointer<Self>         Pointer;
00059   typedef SmartPointer<const Self>   ConstPointer;
00060 
00062   itkNewMacro(Self);
00063  
00065   itkTypeMacro(EMLocalCostFunctionWrapper, SingleValuedCostFunction);
00066 
00069   typedef Superclass::ParametersType      ParametersType;
00070   typedef Superclass::DerivativeType      DerivativeType;
00071   typedef Superclass::MeasureType         MeasureType ;
00072 
00073 
00075   void GetDerivative( const ParametersType & vtkNotUsed(parameters),
00076                       DerivativeType & vtkNotUsed(derivative) ) const
00077   {
00078      itkExceptionMacro( "GetDerivative not supported!" );
00079   }
00080 
00082   MeasureType GetValue( const ParametersType & parameters ) const
00083   {
00084     itkDebugMacro("GetValue( " << parameters << " ) ");
00085     const double* para_double = parameters.data_block();
00086 
00087     if (this->m_Registration) return this->m_Registration->ComputeCostFunction(para_double);
00088     if (this->m_Shape) return this->m_Shape->ComputeCostFunction(para_double);
00089     itkExceptionMacro( "Neither registration nor shape cost function is set!" );
00090   }
00091 
00093   void GetValueAndDerivative( const ParametersType & vtkNotUsed(parameters),
00094                               MeasureType& vtkNotUsed(Value),
00095                               DerivativeType& vtkNotUsed(Derivative) ) const
00096   {
00097     itkExceptionMacro( "GetValueAndDerivative not supported!" );
00098   }
00099 
00100   unsigned int GetNumberOfParameters(void) const {
00101     if (m_Registration) return  this->m_Registration->GetNumberOfTotalParameters(); 
00102     if (m_Shape) return  this->m_Shape->GetPCATotalNumOfShapeParameters(); 
00103     itkExceptionMacro( "Neither registration nor shape cost function is set!" );
00104   }
00105   
00106   void SetRegistrationCostFunction(EMLocalRegistrationCostFunction* init) {assert(!this->m_Shape); this->m_Registration = init;} 
00107   void SetShapeCostFunction(EMLocalShapeCostFunction* init) {assert(!this->m_Registration); this->m_Shape = init;} 
00108 
00109 protected:
00110   EMLocalCostFunctionWrapper() { m_Registration = NULL; m_Shape = NULL;}; 
00111   virtual ~EMLocalCostFunctionWrapper() { };
00112   EMLocalRegistrationCostFunction* m_Registration;
00113   EMLocalShapeCostFunction* m_Shape;
00114 
00115 private:
00116   EMLocalCostFunctionWrapper(const Self&); //purposely not implemented
00117   void operator=(const Self&); //purposely not implemented
00118 
00119 };
00120 
00121 } // end namespace itk
00122 
00123 /* Run optimzation */
00124 // Seperate from itk function bc it is only a wrapper around 
00125 void  itkEMLocalOptimization_Registration_Start(EMLocalRegistrationCostFunction* RegCostFunction,
00126                                                 double* Parameters,
00127                                                 float &vtkNotUsed(Cost))
00128 {
00129   std::cerr << "==================== Start Registration =========================== " << endl;
00130 
00131   RegCostFunction->InitializeCostFunction();
00132   int NumOfFunctionEvaluations;
00133 
00134   {
00135      // EMLocalRegistration::Pointer em = EMLocalRegistration::New();
00136      itk::EMLocalCostFunctionWrapper::Pointer itkRegCostFunction = itk::EMLocalCostFunctionWrapper::New(); 
00137      itkRegCostFunction->SetRegistrationCostFunction(RegCostFunction);
00138      int NumPara     =  itkRegCostFunction->GetNumberOfParameters(); 
00139 
00140      itk::PowellOptimizer::Pointer optimizer = itk::PowellOptimizer::New();
00141      optimizer->SetCostFunction(itkRegCostFunction);
00142      optimizer->SetMaximize( false );
00143 
00144      // From NR 
00145      optimizer->SetStepLength( 1.0 );
00146      optimizer->SetStepTolerance(float(2.0e-4) );
00147      optimizer->SetValueTolerance( 0.01 );
00148      optimizer->SetMaximumIteration( 200 );
00149      optimizer->SetMaximumLineIteration( 100 );
00150 
00151      typedef itk::EMLocalCostFunctionWrapper::ParametersType ParametersType;
00152 
00153      ParametersType InitialPara(NumPara); 
00154      memcpy(InitialPara.data_block(),Parameters,sizeof(double) *  NumPara);
00155     
00156      optimizer->SetInitialPosition(InitialPara);
00157 
00158      try 
00159      {
00160        optimizer->StartOptimization();
00161      }
00162      catch( itk::ExceptionObject & e )
00163      {
00164        std::cerr << "Exception thrown ! " << std::endl;
00165        std::cerr << "An error ocurred during Optimization" << std::endl;
00166        std::cerr << "Location    = " << e.GetLocation()    << std::endl;
00167        std::cerr << "Description = " << e.GetDescription() << std::endl;
00168        return;
00169      }
00170 
00171      NumOfFunctionEvaluations = optimizer->GetCurrentIteration();
00172 
00173      ParametersType finalPosition = optimizer->GetCurrentPosition();
00174      memcpy(Parameters,finalPosition.data_block(), sizeof(double) *  NumPara);
00175   }
00176   RegCostFunction->FinalizeCostFunction(Parameters, NumOfFunctionEvaluations);
00177   std::cerr << "==================== End Registration =========================== " << endl;
00178 }
00179 
00180 
00181 void  itkEMLocalOptimization_Shape_Start(EMLocalShapeCostFunction* ShapeCostFunction, float **PCAShapeParameters, int PCAMaxX, int PCAMinX, 
00182                        int PCAMaxY, int PCAMinY, int PCAMaxZ, int PCAMinZ, int BoundaryMinX, int BoundaryMinY, int BoundaryMinZ, 
00183                        int Boundary_LengthX, int Boundary_LengthY, float** w_m, unsigned char* PCA_ROI, void  **initProbDataPtr, 
00184                        float** initPCAMeanShapePtr, int* initPCAMeanShapeIncY, int *initPCAMeanShapeIncZ, float*** initPCAEigenVectorsPtr, 
00185                        int **initPCAEigenVectorsIncY,  int** initPCAEigenVectorsIncZ, float& Cost) {
00186   std::cerr << "==================== Start Shape Deformation  =========================== " << endl;
00187   std::cerr << "Implementation: ITK" << endl;
00188 
00189   ShapeCostFunction->InitializeCostFunction(PCAMaxX, PCAMinX, PCAMaxY, PCAMinY, PCAMaxZ, PCAMinZ, BoundaryMinX, BoundaryMinY, BoundaryMinZ,
00190                       Boundary_LengthX, Boundary_LengthY, w_m, PCA_ROI, initProbDataPtr, initPCAMeanShapePtr, initPCAMeanShapeIncY, 
00191                       initPCAMeanShapeIncZ, initPCAEigenVectorsPtr, initPCAEigenVectorsIncY, initPCAEigenVectorsIncZ);
00192 
00193   itk::EMLocalCostFunctionWrapper::Pointer itkShapeCostFunction = itk::EMLocalCostFunctionWrapper::New(); 
00194   itkShapeCostFunction->SetShapeCostFunction(ShapeCostFunction);
00195   int NumPara     =  itkShapeCostFunction->GetNumberOfParameters(); 
00196     
00197   itk::PowellOptimizer::Pointer optimizer = itk::PowellOptimizer::New();
00198   optimizer->SetCostFunction(itkShapeCostFunction);
00199   optimizer->SetMaximize( false );
00200 
00201   // From NR 
00202   optimizer->SetStepLength( 1.0 );
00203   optimizer->SetStepTolerance(float(2.0e-4) );
00204   optimizer->SetValueTolerance( 0.01 );
00205   optimizer->SetMaximumIteration( 200 );
00206   optimizer->SetMaximumLineIteration( 100 );
00207 
00208   typedef itk::EMLocalCostFunctionWrapper::ParametersType ParametersType;
00209   
00210   // Transfere Variables
00211   float* PCAParameters  = new float[NumPara]; 
00212   ShapeCostFunction->TransferePCAShapeParametersIntoArray(PCAShapeParameters,PCAParameters); 
00213 
00214   ParametersType InitialPara(NumPara); 
00215   double *InitialPara_db = InitialPara.data_block();
00216   for (int i = 0 ; i < NumPara; i++) InitialPara_db[i] = (double) PCAParameters[i]; 
00217   optimizer->SetInitialPosition(InitialPara);
00218 
00219   // run Algorithm,
00220   try 
00221     {
00222       optimizer->StartOptimization();
00223     }
00224   catch( itk::ExceptionObject & e )
00225     {
00226       std::cerr << "Exception thrown ! " << std::endl;
00227       std::cerr << "An error ocurred during Optimization" << std::endl;
00228       std::cerr << "Location    = " << e.GetLocation()    << std::endl;
00229       std::cerr << "Description = " << e.GetDescription() << std::endl;
00230       return;
00231     }
00232 
00233   Cost = float(optimizer->GetCurrentCost()); 
00234   std::cerr << "Number of Evaluations :" << optimizer->GetCurrentIteration() << endl;
00235 
00236   
00237   const double* finalPosition_db = optimizer->GetCurrentPosition().data_block();
00238   for (int i = 0 ; i < NumPara; i++) PCAParameters[i] = float(finalPosition_db[i]);
00239 
00240   ShapeCostFunction->TransfereArrayIntoPCAShapeParameters(PCAParameters,PCAShapeParameters); 
00241 
00242   delete[] PCAParameters;
00243 
00244   std::cerr << "==================== End Shape Deformation =========================== " << endl;
00245 }
00246 
00247 
00248 #endif

Generated on 6 Apr 2011 for Slicer3 by  doxygen 1.6.1