f.r.e.e. - Flexible Registration and Evaluation Engine: freFEMRegistrationProcessor.txx Source File

00001 /*=========================================================================
00002  
00003   Program:   F.R.E.E. - flexible registration evaluation engine
00004   Version:   v.1.0.0
00005   Date:      $Date: 2006/09/01 12:00:00 $
00006   Module:    $RCSfile: freFEMRegistrationProcessor.h,v $
00007   Language:  C++
00008  
00009  
00010  
00011   Copyright (c) 2007 Ralf o Floca (Department of Medical Informatics,
00012   Institute for Medical Biometry and Informatics, University of Heidelberg,
00013   Germany). All rights reserved.
00014   See FREECopyright.txt or http://www.mi.med.uni-hd.de/free/copyright.htm
00015   for details.
00016  
00017      This software is distributed WITHOUT ANY WARRANTY; without even 
00018      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
00019      PURPOSE.  See the above copyright notices for more information.
00020  
00021 =========================================================================*/
00022 #ifndef __freFEMRegistrationProcessor_txx
00023 #define __freFEMRegistrationProcessor_txx
00024 
00025 #include "freFEMRegistrationProcessor.h"
00026 #include "freOptimizerControllerBase.h"
00027 #include "freTransformSetupAdaptor.h"
00028 
00029 #include "itkPasteImageFilter.h"
00030 #include "itkNormalizeImageFilter.h"
00031 #include "itkVectorCastImageFilter.h"
00032 
00033 namespace FREE
00034 {
00035 
00039 
00040 
00041 
00042 template <unsigned int VImageDimension>
00043 const int
00044 FEMRegistrationProcessor<VImageDimension>::
00045 GetResolutionLevelCount() const
00046 {
00047     int iLevelCount;
00048     if (this->GetComponentSetup())
00049         if (this->GetComponentSetup()->Parameters().GetParameterValue(cParam_FEMRegLevels,iLevelCount))
00050             return iLevelCount;
00051     return 1;
00052 };
00053 
00054 template <unsigned int VImageDimension>
00055 const long
00056 FEMRegistrationProcessor<VImageDimension>::
00057 GetMaxIterationCount()
00058 {
00059     if (this->GetComponentSetup())
00060     {
00061         ComponentSetup* pComponentSetup = this->GetComponentSetup()->Components().GetElement(cComp_MainOptimizer);
00062         if(pComponentSetup)
00063         {
00064             long lIterations;
00065             if (pComponentSetup->Parameters().GetParameterValue(cParam_Iterations,lIterations,0,this->m_iResolutionLevel))
00066                 return lIterations;
00067         }
00068     }
00069     return -1;
00070 };
00071 
00072 template <unsigned int VImageDimension>
00073 void
00074 FEMRegistrationProcessor<VImageDimension>::
00075 ComputeRegistration()
00076 {
00077     long optimizerObserverID = 0;
00078     try
00079     {
00080         Superclass::ComputeRegistration();
00081         typename HistoMatchFilterType::Pointer histoMatch = HistoMatchFilterType::New();
00082 
00083         //Set initials parameters
00084         if(this->fnOnProgress.IsNotNull())
00085         {
00086             this->fnOnProgress->Execute(RSTInitProcessor, "Setting initial setup values of components",this);
00087         }
00088 
00089         SessionBuilder::ActualizeComponent(this->m_pComponentCache, this->m_pSessionInfo);
00090 
00091         //Extract needed ImageParts;
00092         if(this->fnOnProgress.IsNotNull())
00093             this->fnOnProgress->Execute(RSTInitProcessor, "Prepare images",this);
00094         typename ImageType::Pointer smpMovingPart = this->m_smpMovingImage;
00095         typename ImageType::Pointer smpFixedPart = this->m_smpFixedImage;
00096 
00097         bool bHistoMatch;
00098         if(!this->GetComponentSetup()->Parameters().GetParameterValue(cParam_HistoMatch,bHistoMatch))
00099             throwExceptionMacro("Missing parameter: " << cParam_HistoMatch);
00100         long lHistoLevels;
00101         if(!this->GetComponentSetup()->Parameters().GetParameterValue(cParam_HistoLevels,lHistoLevels))
00102             throwExceptionMacro("Missing parameter: " << cParam_HistoLevels);
00103 
00104         long lHistoPoints;
00105         if(!this->GetComponentSetup()->Parameters().GetParameterValue(cParam_HistoMatchPoints,lHistoPoints))
00106             throwExceptionMacro("Missing parameter: " << cParam_HistoMatchPoints);
00107 
00108         bool bHistoThreshold;
00109         if(!this->GetComponentSetup()->Parameters().GetParameterValue(cParam_HistoThreshold,bHistoThreshold))
00110             throwExceptionMacro("Missing parameter: " << cParam_HistoThreshold);
00111 
00112         //Connect images
00113         if(bHistoMatch)
00114         {
00115             if(this->fnOnProgress.IsNotNull())
00116                 this->fnOnProgress->Execute(RSTInitProcessor, "Histogram match enabled",this);
00117 
00118             histoMatch->SetInput( smpMovingPart );
00119             histoMatch->SetReferenceImage( smpFixedPart );
00120 
00121             histoMatch->SetNumberOfHistogramLevels( lHistoLevels );
00122             histoMatch->SetNumberOfMatchPoints( lHistoPoints );
00123             if ( bHistoThreshold )
00124                 histoMatch->ThresholdAtMeanIntensityOn();
00125 
00126             //Calculate histogram matching
00127             histoMatch->Update();
00128             //Connect it to registration
00129             this->smpRegistration->SetFixedImage(smpFixedPart);
00130             this->smpRegistration->SetMovingImage(histoMatch->GetOutput());
00131         }
00132         else
00133         {
00134             this->smpRegistration->SetFixedImage(smpFixedPart);
00135             this->smpRegistration->SetMovingImage(smpMovingPart);
00136         };
00137 
00138         SetRegistrationParameters();
00139 
00140         SetFEMElements();
00141 
00142         //Register observer
00143         m_IterationObserver->fnOnNotify = m_IterationNotificationEvent;
00144         optimizerObserverID = this->smpRegistration->AddObserver(itk::IterationEvent(), this->m_IterationObserver);
00145     }
00146     catchAllNPassMacro("Unknown Error while preparing the registration.");
00147 
00148     try
00149     {
00150         if(this->fnOnProgress.IsNotNull())
00151             this->fnOnProgress->Execute(RSTProcessing, "Solving FEM",this);
00152         this->m_lCurIteration = 0;
00153 
00154         //Call the OnNextLevel one time before starting the registration.
00155         //even a singel resolution Registration has at least one level to notify
00156         if(this->fnOnNextLevel.IsNotNull())
00157             this->fnOnNextLevel->Execute(0,this);
00158 
00159         this->smpRegistration->RunRegistration();
00160 
00161         this->smpRegistration->RemoveObserver(optimizerObserverID);
00162     }
00163     catchAllNPassMacro("Unknown Error while itk processing the registration.");
00164 };
00165 
00166 template <unsigned int VImageDimension>
00167 FEMRegistrationProcessor<VImageDimension>::
00168 FEMRegistrationProcessor()
00169 {
00170     smpRegistration = RegistrationType::New();
00171     smpMetric = NULL;
00172     m_IterationObserver = IterationObserver::New();
00173     m_IterationNotificationEvent = NotificationEvent<FEMRegistrationProcessor>::New(this, &Self::OnNewIteration);
00174 };
00175 
00176 template <unsigned int VImageDimension>
00177 FEMRegistrationProcessor<VImageDimension>::
00178 ~FEMRegistrationProcessor()
00179 {}
00180 ;
00181 
00182 template <unsigned int VImageDimension>
00183 void
00184 FEMRegistrationProcessor<VImageDimension>::
00185 SetRegistrationParameters()
00186 {
00187     std::string sParam;
00188 
00189     ComponentSetup* pSetup = this->GetComponentSetup();
00190 
00191     std::string sSetupFile;
00192     if(!pSetup->Parameters().GetParameterValue(cParam_FEMSetup,sSetupFile))
00193         throwExceptionMacro("Missing parameter: " <<cParam_FEMSetup);
00194 
00195     int iRegLevels;
00196     if(!pSetup->Parameters().GetParameterValue(cParam_FEMRegLevels,iRegLevels))
00197         throwExceptionMacro("Missing parameter: " <<cParam_FEMRegLevels);
00198 
00199     double dAlpha;
00200     if(!pSetup->Parameters().GetParameterValue(cParam_FEMAlpha,dAlpha))
00201         throwExceptionMacro("Missing parameter: " << cParam_FEMAlpha);
00202 
00203     bool bDescentDirectionMinimize;
00204     if(!pSetup->Parameters().GetParameterValue(cParam_FEMDescentDirection,bDescentDirectionMinimize))
00205         throwExceptionMacro("Missing parameter: " << cParam_FEMDescentDirection);
00206 
00207     int iDoLineSearch;
00208     if(!pSetup->Parameters().GetParameterValue(cParam_FEMDoLineSearch,iDoLineSearch))
00209         throwExceptionMacro("Missing parameter: " << cParam_FEMDoLineSearch);
00210 
00211     double dTimeStep;
00212     if(!pSetup->Parameters().GetParameterValue(cParam_FEMTimeStep,dTimeStep))
00213         throwExceptionMacro("Missing parameter: " << cParam_FEMTimeStep);
00214 
00215     double dEnergyReduction;
00216     if(!pSetup->Parameters().GetParameterValue(cParam_FEMEnergyReduction,dEnergyReduction))
00217         throwExceptionMacro("Missing parameter: " << cParam_FEMEnergyReduction);
00218 
00219     std::string sLandmarkFile;
00220     if(!pSetup->Parameters().GetParameterValue(cParam_FEMLandmarkFile,sLandmarkFile))
00221         throwExceptionMacro("Missing parameter: " << cParam_FEMLandmarkFile);
00222 
00223     int iLandmarkType;
00224     if(!pSetup->Parameters().GetParameterValue(cParam_FEMLandmarkType,iLandmarkType))
00225         throwExceptionMacro("Missing parameter: " << cParam_FEMLandmarkType);
00226 
00227     int iMetricType;
00228     if(!pSetup->Parameters().GetParameterValue(cParam_FEMMetric,iMetricType))
00229         throwExceptionMacro("Missing parameter: " << cParam_FEMMetric);
00230 
00231     int iScaleGradient;
00232     if(!pSetup->Parameters().GetParameterValue(cParam_FEMScaleGradient,iScaleGradient))
00233         throwExceptionMacro("Missing parameter: " << cParam_FEMScaleGradient);
00234 
00235     int iRegridding;
00236     if(!pSetup->Parameters().GetParameterValue(cParam_FEMRegridding,iRegridding))
00237         throwExceptionMacro("Missing parameter: " << cParam_FEMRegridding);
00238 
00239     if (sSetupFile!="")
00240     {
00241         smpRegistration->SetConfigFileName(sSetupFile.c_str());
00242         if ( !smpRegistration->ReadConfigFile( (smpRegistration->GetConfigFileName()).c_str() ) )
00243             throwExceptionMacro("Error while reading fem setup file: " << smpRegistration->GetConfigFileName());
00244     }
00245     else
00246     {
00247         //Use parameters of free-file
00248         //Multi Res is always used. With one Level, it is a single resolution registration
00249         smpRegistration->DoMultiRes(true);
00250         smpRegistration->SetNumLevels(iRegLevels);
00251         smpRegistration->SetMaxLevel(iRegLevels);
00252 
00253         for (int iLevel=0; iLevel<iRegLevels; iLevel++)
00254         {
00255             int iNrIterations;
00256             if(!pSetup->Parameters().GetParameterValue(cParam_Iterations,iNrIterations,0,iLevel))
00257                 throwExceptionMacro("Missing parameter: " << cParam_Iterations);
00258 
00259             int iPPE;
00260             if(!pSetup->Parameters().GetParameterValue(cParam_FEMPPE,iPPE,0,iLevel))
00261                 throwExceptionMacro("Missing parameter: " << cParam_FEMPPE);
00262 
00263             double dElasticity;
00264             if(!pSetup->Parameters().GetParameterValue(cParam_FEMElasticity,dElasticity,0,iLevel))
00265                 throwExceptionMacro("Missing parameter: " << cParam_FEMElasticity);
00266 
00267             double dDensity;
00268             if(!pSetup->Parameters().GetParameterValue(cParam_FEMDensity,dDensity,0,iLevel))
00269                 throwExceptionMacro("Missing parameter: " << cParam_FEMDensity);
00270 
00271             double dEnergyScale;
00272             if(!pSetup->Parameters().GetParameterValue(cParam_FEMEnergyScale,dEnergyScale,0,iLevel))
00273                 throwExceptionMacro("Missing parameter: " << cParam_FEMEnergyScale);
00274 
00275             int iNrIntPoints;
00276             if(!pSetup->Parameters().GetParameterValue(cParam_FEMNrIntPoints,iNrIntPoints,0,iLevel))
00277                 throwExceptionMacro("Missing parameter: " << cParam_FEMNrIntPoints);
00278 
00279             int iWidthOfMetric;
00280             if(!pSetup->Parameters().GetParameterValue(cParam_FEMWidthOfMetric,iWidthOfMetric,0,iLevel))
00281                 throwExceptionMacro("Missing parameter: " << cParam_FEMWidthOfMetric);
00282 
00283             smpRegistration->SetMaximumIterations(iNrIterations,iLevel);
00284             smpRegistration->SetMeshPixelsPerElementAtEachResolution(iPPE,iLevel);
00285             smpRegistration->SetElasticity(dElasticity,iLevel);
00286             smpRegistration->SetRho(dDensity,iLevel);
00287             smpRegistration->SetGamma(dEnergyScale,iLevel);
00288             smpRegistration->SetNumberOfIntegrationPoints(iNrIntPoints,iLevel);
00289             smpRegistration->SetWidthOfMetricRegion(iWidthOfMetric,iLevel);
00290         };
00291 
00292         smpRegistration->SetAlpha(dAlpha);
00293         smpRegistration->SetTimeStep(dTimeStep);
00294         smpRegistration->SetEnergyReductionFactor(dEnergyReduction);
00295         smpRegistration->EmployRegridding(iRegridding);
00296 
00297         if (iLandmarkType==0)
00298             smpRegistration->UseLandmarks(false);
00299         else
00300         {
00301             if (iLandmarkType==2)
00302             {
00303                 sLandmarkFile = GetGeneralFREEPath()+"fem.lms";
00304                 SaveLandmarksToFile(sLandmarkFile);
00305             }
00306 
00307             smpRegistration->UseLandmarks(true);
00308             smpRegistration->SetLandmarkFile(sLandmarkFile.c_str());
00309         };
00310 
00311         //Connect metric
00312         smpRegistration->SetTemp(iScaleGradient);
00313         smpRegistration->ChooseMetric(iMetricType);
00314         if (bDescentDirectionMinimize)
00315             smpRegistration->SetDescentDirectionMinimize();
00316         else
00317             smpRegistration->SetDescentDirectionMaximize();
00318 
00319         smpRegistration->DoLineSearch(iDoLineSearch);
00320     }
00321 };
00322 
00323 template <unsigned int VImageDimension>
00324 void
00325 FEMRegistrationProcessor<VImageDimension>::
00326 SaveLandmarksToFile(const std::string& sFileName)
00327 {
00328     std::ios_base::openmode  iOpenFlag = std::ios_base::out | std::ios_base::trunc;
00329     ;
00330 
00331     std::ofstream landmarkFile;
00332     bool result = false;
00333 
00334     landmarkFile.open(sFileName.c_str(), iOpenFlag );
00335 
00336     if (!landmarkFile.is_open())
00337         throwExceptionMacro("Unable to open landmark file: " << sFileName);
00338 
00339     landmarkFile<<"% Landmarkfile for F.R.E.E. / itk" << std::endl;
00340 
00341     for (int iLevel=0; iLevel< this->GetComponentSetup()->Parameters().ParameterLayerCount(cParam_FEMLandmarks); iLevel++)
00342     {
00343         std::string sLine = FREE::Convert::ToStr(VImageDimension);
00344         landmarkFile<<"<LoadLandmark>" << std::endl;
00345         landmarkFile<<"-1   % GN of element (undefined)" << std::endl;
00346         for (int iIndex=0; iIndex< VImageDimension; iIndex++)
00347         {
00348             double dValue;
00349             this->GetComponentSetup()->Parameters().GetParameterValue(cParam_FEMLandmarks,dValue,iIndex,iLevel);
00350             sLine = sLine + "    " + FREE::Convert::ToStr((float)dValue);
00351         }
00352         landmarkFile<< sLine << "   % Undeformed landmark\n" << std::endl;
00353 
00354         sLine = FREE::Convert::ToStr(VImageDimension);
00355         for (int iIndex=VImageDimension; iIndex< this->GetComponentSetup()->Parameters().ParameterSize(cParam_FEMLandmarks); iIndex++)
00356         {
00357             double dValue;
00358             this->GetComponentSetup()->Parameters().GetParameterValue(cParam_FEMLandmarks,dValue,iIndex,iLevel);
00359             sLine = sLine + "    " + FREE::Convert::ToStr((float)dValue);
00360         }
00361         landmarkFile<< sLine << "   % Deformed landmark" << std::endl;
00362 
00363         landmarkFile<< "     1.000000" <<std::endl;
00364     }
00365 
00366     landmarkFile<< "<END>" << std::endl;
00367 
00368     landmarkFile.close();
00369 }
00370 
00371 template <unsigned int VImageDimension>
00372 void
00373 FEMRegistrationProcessor<VImageDimension>::
00374 SetFEMElements()
00375 {
00376     itk::fem::MaterialLinearElasticity::Pointer m;
00377     m = itk::fem::MaterialLinearElasticity::New();
00378     m->GN = 0;                                                            // Global number of the material
00379     m->E = smpRegistration->GetElasticity();  // Young's modulus -- used in the membrane
00380     m->A = 1.0;                                                           // Cross-sectional area
00381     m->h = 1.0;                                                           // Thickness
00382     m->I = 1.0;                                                           // Moment of inertia
00383     m->nu = 0.;                                                           // Poisson's ratio -- DONT CHOOSE 1.0!!
00384     m->RhoC = 1.0;                                                        // Density
00385 
00386     if (VImageDimension==2)
00387     {
00388         //2D Version
00389         Element2DType::LoadImplementationFunctionPointer fp =
00390             &itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
00391         Dispatcher2DType::RegisterVisitor((ImageLoadType*)0,fp);
00392 
00393         Element2DType2::LoadImplementationFunctionPointer fp2 =
00394             &itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
00395         Dispatcher2DType2::RegisterVisitor((ImageLoadType*)0,fp2);
00396 
00397         // Create the element type
00398         Element2DType::Pointer e1=Element2DType::New();
00399         e1->m_mat=dynamic_cast<itk::fem::MaterialLinearElasticity*>( m );
00400         smpRegistration->SetElement(e1);
00401         smpRegistration->SetMaterial(m);
00402     }
00403     else
00404     {
00405         //3D Version
00406         Element3DType::LoadImplementationFunctionPointer fp =
00407             &itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
00408         Dispatcher3DType::RegisterVisitor((ImageLoadType*)0,fp);
00409 
00410         Element3DType2::LoadImplementationFunctionPointer fp2 =
00411             &itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
00412         Dispatcher3DType2::RegisterVisitor((ImageLoadType*)0,fp2);
00413 
00414         // Create the element type
00415         Element3DType::Pointer e1=Element3DType::New();
00416         e1->m_mat=dynamic_cast<itk::fem::MaterialLinearElasticity*>( m );
00417         smpRegistration->SetElement(e1);
00418         smpRegistration->SetMaterial(m);
00419     }
00420 };
00421 
00422 template <unsigned int VImageDimension>
00423 typename FEMRegistrationProcessor<VImageDimension>::TransformFieldPointer
00424 FEMRegistrationProcessor<VImageDimension>::
00425 ComputeTransformationField(const PointType& origin, const RegionType& region, const SpacingType& spacing)
00426 {
00427     long optimizerObserverID = 0;
00428     if(this->fnOnProgress.IsNotNull())
00429         this->fnOnProgress->Execute(RSTFinalizing, "Computing transformation field",this);
00430 
00431     TransformFieldPointer smpFEMField;
00432     try
00433     {
00434         typedef itk::VectorCastImageFilter< itk::Image<itk::Vector<float,VImageDimension>,VImageDimension>, TransformFieldType > VectorCastFilterType;
00435         typename VectorCastFilterType::Pointer  deformCaster =  VectorCastFilterType::New();
00436 
00437         deformCaster->SetInput( this->smpRegistration->GetDeformationField());
00438         deformCaster->Update();
00439 
00440         smpFEMField = deformCaster->GetOutput();
00441 
00442         //the deform field computed be itk has always a spacing of 1 in all dimensions
00443         smpFEMField->SetSpacing(this->m_smpFixedImage->GetSpacing());
00444     }
00445     catchAllNPassMacro("Unknown error, while casting transformation field.");
00446 
00447     try
00448     {
00449     if (!( (region==smpFEMField->GetLargestPossibleRegion())&
00450            (spacing==smpFEMField->GetSpacing())&
00451            (origin==smpFEMField->GetOrigin()) ) )
00452       { // have to resize the field.
00453         //Create the field in correct size
00454         typename TransformFieldType::Pointer smpField = TransformFieldType::New();
00455         smpField->SetOrigin(origin);
00456         smpField->SetSpacing(spacing);
00457         smpField->SetRegions(region);
00458         smpField->Allocate();
00459         //Reset the transformation field to no transformation (all vectors components are 0)
00460         typename TransformFieldType::PixelType nullVector;
00461         for (int iIndex=0; iIndex < VImageDimension; iIndex++)
00462             nullVector[iIndex]=0;
00463 
00464         // iterator for the output image
00465         itk::ImageRegionIteratorWithIndex<TransformFieldType> outputIt(
00466           smpField, smpField->GetLargestPossibleRegion() );
00467 
00468         // interpolator
00469         typedef itk::LinearInterpolateImageFunction<TransformFieldType,ScalarType> InterpolatorType;
00470         typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
00471         interpolator->SetInputImage(smpFEMField);
00472 
00473         typename ImageType::IndexType index;
00474         typename ImageType::PointType point;
00475         typename TransformFieldType::PixelType value;
00476 
00477         while( !outputIt.IsAtEnd() )
00478         {
00479           // get the output image index
00480           index = outputIt.GetIndex();
00481           smpField->TransformIndexToPhysicalPoint( index, point );
00482 
00483           if( interpolator->IsInsideBuffer( point ) )
00484           {
00485             value = interpolator->Evaluate( point );
00486             outputIt.Set( value );
00487           }
00488           else
00489           {
00490             outputIt.Set( nullVector );
00491           }  
00492 
00493           ++outputIt;
00494         }
00495 
00496         return smpField;
00497       }
00498     }
00499     catchAllNPassMacro("Unknown error, while pasting transformation field.");
00500 
00501     return smpFEMField;
00502 }
00503 ;
00504 
00505 template <unsigned int VImageDimension>
00506 void
00507 FEMRegistrationProcessor<VImageDimension>::
00508 OnNewIteration(void* pSender, long threadID)
00509 {
00510     if(this->fnOnNextIteration.IsNotNull())
00511         this->fnOnNextIteration->Execute(this->m_lCurIteration, NULL, this, 0);
00512     this->m_lCurIteration++;
00513 };
00514 
00515 }//End of Namespace free
00516 
00517 #endif