mitk: mitkPyramidalRegistrationMethod.txx Source File

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003 Program:   Medical Imaging & Interaction Toolkit
00004 Language:  C++
00005 Date:      $Date: 2007-09-05 13:45:48 +0200 (Mi, 05 Sep 2007) $
00006 Version:   $Revision: 9502 $
00007 
00008 Copyright (c) German Cancer Research Center, Division of Medical and
00009 Biological Informatics. All rights reserved.
00010 See MITKCopyright.txt or https://www.mitk.org/copyright.html for details.
00011 
00012 This software is distributed WITHOUT ANY WARRANTY; without even
00013 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00014 PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 
00018 #include "mitkPyramidalRegistrationMethod.h"
00019 
00020 #include "mitkMetricFactory.h"
00021 #include "mitkTransformFactory.h"
00022 #include "mitkOptimizerFactory.h"
00023 #include "mitkTransformParameters.h"
00024 #include "mitkMetricParameters.h"
00025 #include <mitkITKImageImport.h>
00026 #include "mitkRigidRegistrationObserver.h"
00027 
00028 #include "itkResampleImageFilter.h"
00029 #include "itkImageRegistrationMethod.h"
00030 #include "itkMultiResolutionImageRegistrationMethod.h"
00031 #include "itkLinearInterpolateImageFunction.h"
00032 #include "itkNearestNeighborInterpolateImageFunction.h"
00033 #include "itkImage.h"
00034 
00035 #include "itkImageToImageMetric.h"
00036 #include "itkSingleValuedNonLinearOptimizer.h"
00037 
00038 
00039 #include "itkDiscreteGaussianImageFilter.h"
00040 
00041 #include "mitkRegistrationInterfaceCommand.h"
00042 
00043 #include <itkRescaleIntensityImageFilter.h>
00044 #include <itkHistogramMatchingImageFilter.h>
00045 #include <itkCastImageFilter.h>
00046 
00047 #include <itkStatisticsImageFilter.h>
00048 #include <itkIntensityWindowingImageFilter.h>
00049 
00050 namespace mitk {
00051 
00052  
00053   template<typename TPixel, unsigned int VImageDimension>
00054   void PyramidalRegistrationMethod::GenerateData2(itk::Image<TPixel, VImageDimension>* itkImage1)
00055   {    
00056 
00057     // TYPEDEFS        
00058     typedef itk::Image< TPixel, VImageDimension >             ImageType;
00059                     
00060                 typedef float                                             InternalPixelType;
00061                 typedef itk::Image< InternalPixelType, VImageDimension >  InternalImageType;
00062 
00063     typedef typename itk::Transform
00064       < double, VImageDimension, VImageDimension >            TransformType;
00065     typedef mitk::TransformFactory<InternalPixelType, 
00066                                    VImageDimension>           TransformFactoryType;
00067     typedef itk::LinearInterpolateImageFunction
00068       < InternalImageType, double >                           InterpolatorType;
00069     typedef mitk::MetricFactory
00070       <InternalPixelType, VImageDimension>                    MetricFactoryType;
00071     typedef itk::RecursiveMultiResolutionPyramidImageFilter<
00072                                                                                                                                                         InternalImageType,
00073                                                                                                                                                         InternalImageType >         ImagePyramidType;
00074     typedef itk::DiscreteGaussianImageFilter
00075       <ImageType, InternalImageType>                          GaussianFilterType;    
00076    
00077 
00078 
00079     typedef itk::MultiResolutionImageRegistrationMethod< 
00080                                                                                                                                                         InternalImageType, 
00081                                                                                                                                                         InternalImageType >     RegistrationType;       
00082     typedef RegistrationInterfaceCommand
00083       <RegistrationType, TPixel>                              CommandType;
00084 
00085     typedef itk::CastImageFilter<ImageType, 
00086                                   InternalImageType>          CastImageFilterType;
00087 
00088     itk::Array<double> initialParameters;
00089     if(m_TransformParameters->GetInitialParameters().size())
00090     {
00091       initialParameters = m_TransformParameters->GetInitialParameters();
00092     }
00093    
00094     // LOAD PARAMETERS
00095     itk::Array<double> transformValues = m_Preset->getTransformValues(m_Presets[0]);
00096     itk::Array<double> metricValues = m_Preset->getMetricValues(m_Presets[0]);
00097     itk::Array<double> optimizerValues = m_Preset->getOptimizerValues(m_Presets[0]);
00098     m_TransformParameters = ParseTransformParameters(transformValues);
00099     m_MetricParameters = ParseMetricParameters(metricValues);
00100     m_OptimizerParameters = ParseOptimizerParameters(optimizerValues);                  
00101 
00102 
00103     // The fixed and the moving image
00104                 typename InternalImageType::Pointer fixedImage = InternalImageType::New();
00105     typename InternalImageType::Pointer movingImage = InternalImageType::New();
00106 
00107   
00108 
00109 
00110                 mitk::CastToItkImage(m_ReferenceImage, fixedImage); 
00111 
00112      // Blur the moving image        
00113     if(m_BlurMovingImage)
00114     {
00115       typename GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New();
00116       gaussianFilter->SetInput(itkImage1);
00117       gaussianFilter->SetVariance(6.0);
00118       gaussianFilter->SetMaximumError( 0.1 );
00119       //gaussianFilter->SetMaximumKernelWidth ( 3 );
00120       gaussianFilter->Update();
00121       movingImage = gaussianFilter->GetOutput();
00122     } else{
00123       typename CastImageFilterType::Pointer castImageFilter = CastImageFilterType::New();
00124       castImageFilter->SetInput(itkImage1);
00125       castImageFilter->Update();
00126       movingImage = castImageFilter->GetOutput();
00127     }
00128     
00129     if(m_MatchHistograms)
00130     {
00131       typedef itk::RescaleIntensityImageFilter<InternalImageType,InternalImageType> FilterType;   
00132       typedef itk::HistogramMatchingImageFilter<InternalImageType,InternalImageType> HEFilterType;
00133 
00134       typename FilterType::Pointer inputRescaleFilter = FilterType::New();  
00135       typename FilterType::Pointer referenceRescaleFilter = FilterType::New();  
00136 
00137       referenceRescaleFilter->SetInput(fixedImage);
00138       inputRescaleFilter->SetInput(movingImage);
00139 
00140       const float desiredMinimum =  0.0;
00141       const float desiredMaximum =  255.0;
00142       
00143       referenceRescaleFilter->SetOutputMinimum( desiredMinimum );
00144       referenceRescaleFilter->SetOutputMaximum( desiredMaximum );
00145       referenceRescaleFilter->UpdateLargestPossibleRegion();  
00146       inputRescaleFilter->SetOutputMinimum( desiredMinimum );
00147       inputRescaleFilter->SetOutputMaximum( desiredMaximum );
00148       inputRescaleFilter->UpdateLargestPossibleRegion();
00149 
00150       // Histogram match the images
00151       typename HEFilterType::Pointer intensityEqualizeFilter = HEFilterType::New();
00152 
00153       intensityEqualizeFilter->SetReferenceImage( inputRescaleFilter->GetOutput() );
00154       intensityEqualizeFilter->SetInput( referenceRescaleFilter->GetOutput() );
00155       intensityEqualizeFilter->SetNumberOfHistogramLevels( 64 );
00156       intensityEqualizeFilter->SetNumberOfMatchPoints( 12 );
00157       intensityEqualizeFilter->ThresholdAtMeanIntensityOn();
00158       intensityEqualizeFilter->Update();
00159 
00160       //fixedImage = referenceRescaleFilter->GetOutput();
00161       //movingImage = IntensityEqualizeFilter->GetOutput();
00162 
00163       fixedImage = intensityEqualizeFilter->GetOutput();
00164       movingImage = inputRescaleFilter->GetOutput();
00165      
00166     }
00167   
00168     
00169     
00170     typename TransformFactoryType::Pointer transFac = TransformFactoryType::New();
00171     transFac->SetTransformParameters(m_TransformParameters);
00172     transFac->SetFixedImage(fixedImage);
00173     transFac->SetMovingImage(movingImage);
00174     typename TransformType::Pointer transform = transFac->GetTransform();
00175     
00176                 
00177     typename InterpolatorType::Pointer   interpolator  = InterpolatorType::New();         
00178     typename MetricFactoryType::Pointer metFac = MetricFactoryType::New();
00179     metFac->SetMetricParameters(m_MetricParameters);
00180 
00181     
00182     typename OptimizerFactory::Pointer optFac = OptimizerFactory::New();
00183     optFac->SetOptimizerParameters(m_OptimizerParameters);
00184     optFac->SetNumberOfTransformParameters(transform->GetNumberOfParameters());
00185     typename OptimizerType::Pointer optimizer = optFac->GetOptimizer(); 
00186 
00187     // optimizer scales
00188     if (m_TransformParameters->GetUseOptimizerScales())
00189     {
00190       itk::Array<double> optimizerScales = m_TransformParameters->GetScales();
00191       typename OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
00192       for (unsigned int i = 0; i < scales.Size(); i++)
00193       {
00194         scales[i] = optimizerScales[i];
00195       }
00196       optimizer->SetScales( scales );
00197     }
00198         
00199                 typename ImagePyramidType::Pointer fixedImagePyramid =  ImagePyramidType::New();
00200                 typename ImagePyramidType::Pointer movingImagePyramid = ImagePyramidType::New();  
00201     
00202     if(m_FixedSchedule.size() > 0 && m_MovingSchedule.size() > 0)
00203     {
00204       fixedImagePyramid->SetSchedule(m_FixedSchedule);
00205       movingImagePyramid->SetSchedule(m_MovingSchedule);
00206       // Otherwise just use the default schedule
00207     }
00208 
00209 
00210     
00211 
00212                 typename RegistrationType::Pointer registration = RegistrationType::New();
00213                 registration->SetOptimizer(           optimizer     );
00214                 registration->SetTransform(           transform     );
00215                 registration->SetInterpolator(        interpolator  );
00216     registration->SetMetric(              metFac->GetMetric()  );        
00217     registration->SetFixedImagePyramid(   fixedImagePyramid );
00218                 registration->SetMovingImagePyramid(  movingImagePyramid );
00219           registration->SetFixedImage(          fixedImage );     
00220                 registration->SetMovingImage(         movingImage );    
00221                 registration->SetFixedImageRegion(    fixedImage->GetBufferedRegion() );     
00222     
00223     if(transFac->GetTransformParameters()->GetInitialParameters().size())
00224     {
00225       registration->SetInitialTransformParameters( transFac->GetTransformParameters()->GetInitialParameters() );
00226     }
00227     else
00228     {
00229       itk::Array<double> zeroInitial;
00230       zeroInitial.set_size(transform->GetNumberOfParameters());
00231       zeroInitial.fill(0.0);
00232       zeroInitial[0] = 1.0;
00233       zeroInitial[4] = 1.0;
00234       zeroInitial[8] = 1.0;
00235       registration->SetInitialTransformParameters( zeroInitial );
00236     }
00237 
00238    
00239     if(m_UseMask)
00240     {
00241       itk::ImageMaskSpatialObject< VImageDimension>* mask = 
00242         dynamic_cast< itk::ImageMaskSpatialObject< VImageDimension>* > ( m_BrainMask.GetPointer() );
00243       registration->GetMetric()->SetFixedImageMask(mask);
00244     }
00245     
00246     // registering command observer with the optimizer        
00247                 if (m_Observer.IsNotNull())
00248                 {
00249                   m_Observer->AddStepsToDo(20);
00250                   optimizer->AddObserver(itk::AnyEvent(), m_Observer);
00251                   registration->AddObserver(itk::AnyEvent(), m_Observer);
00252                   transform->AddObserver(itk::AnyEvent(), m_Observer);
00253                 }                       
00254 
00255     typename CommandType::Pointer command = CommandType::New();
00256     command->observer = m_Observer;
00257     command->m_Presets = m_Presets;
00258     command->m_UseMask = m_UseMask;
00259     command->m_BrainMask = m_BrainMask;
00260 
00261                 registration->AddObserver( itk::IterationEvent(), command );        
00262     registration->SetSchedules(m_FixedSchedule, m_MovingSchedule);        
00263     
00264 
00265     // Start the registration process
00266                 try 
00267                 { 
00268                   registration->StartRegistration(); 
00269                 } 
00270                 catch( itk::ExceptionObject & err ) 
00271                 { 
00272                   std::cout << "ExceptionObject caught !" << std::endl; 
00273                   std::cout << err << std::endl;                                
00274                 }
00275     if (m_Observer.IsNotNull())
00276     {
00277       optimizer->RemoveAllObservers();
00278       registration->RemoveAllObservers();
00279       transform->RemoveAllObservers();
00280       m_Observer->SetRemainingProgress(15);
00281     }
00282     if (m_Observer.IsNotNull())
00283     {
00284       m_Observer->SetRemainingProgress(5);
00285     }
00286 
00287   } 
00288 
00289  
00290 
00291 } // end namespace