mitk: mitkBSplineRegistration.cpp Source File

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003 Program:   Medical Imaging & Interaction Toolkit
00004 Language:  C++
00005 Date:      $Date: 2009-05-05 19:02:32 +0200 (Tue, 05 May 2009) $
00006 Version:   $Revision: 17105 $
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 "itkImageFileWriter.h"
00019 #include "itkWarpImageFilter.h"
00020 #include "itkImageRegionIterator.h"
00021 
00022 #include "mitkBSplineRegistration.h"
00023 #include "itkBSplineDeformableTransform.h"
00024 //#include "itkLBFGSOptimizer.h"
00025 #include "itkMeanSquaresImageToImageMetric.h"
00026 #include "itkMattesMutualInformationImageToImageMetric.h"
00027 #include "itkResampleImageFilter.h"
00028 #include "itkImageRegistrationMethod.h"
00029 
00030 #include "itkBSplineDeformableTransformInitializer.h"
00031 
00032 #include "mitkOptimizerFactory.h"
00033 #include "mitkMetricFactory.h"
00034 
00035 #include <itkRescaleIntensityImageFilter.h>
00036 #include <itkHistogramMatchingImageFilter.h>
00037 
00038 
00039 namespace mitk {
00040 
00041   
00042   BSplineRegistration::BSplineRegistration():
00043     m_Iterations(50),    
00044     m_ResultName("deformedImage.mhd"),   
00045     m_SaveResult(true),
00046     m_SaveDeformationField(false),
00047     m_UpdateInputImage(false),
00048     m_MatchHistograms(true),    
00049     m_Metric(0)
00050   {
00051     m_Observer = mitk::RigidRegistrationObserver::New();
00052   }
00053 
00054   BSplineRegistration::~BSplineRegistration()
00055   {
00056   }
00057 
00058   void BSplineRegistration::SetNumberOfIterations(int iterations)
00059   {
00060     m_Iterations = iterations;
00061   }  
00062 
00063   void BSplineRegistration::SetSaveResult(bool saveResult)
00064   {
00065     m_SaveResult = saveResult;
00066   }
00067 
00068   void BSplineRegistration::SetResultFileName(const char* resultName)
00069   {
00070     m_ResultName = resultName;
00071   }
00072 
00073 
00074   
00075 
00076   template < typename TPixel, unsigned int VImageDimension >
00077     void BSplineRegistration::GenerateData2( itk::Image<TPixel, VImageDimension>* itkImage1)
00078   {
00079     std::cout << "start bspline registration" << std::endl;
00080     
00081     // Typedefs
00082     typedef typename itk::Image< TPixel, VImageDimension >  InternalImageType;
00083        
00084     typedef typename itk::Vector< float, VImageDimension >    VectorPixelType;
00085     typedef typename itk::Image<  VectorPixelType, VImageDimension > DeformationFieldType;
00086 
00087     typedef itk::BSplineDeformableTransform<
00088                                 double,
00089                                 VImageDimension,
00090                                 3 >                         TransformType;
00091 
00092     typedef typename TransformType::ParametersType          ParametersType;
00093 
00094 
00095 
00096     
00097     //typedef itk::LBFGSOptimizer                             OptimizerType;
00098     typedef itk::SingleValuedNonLinearOptimizer             OptimizerType;
00099     //typedef itk::SingleValuedCostFunction                   MetricType;
00100 
00101     typedef itk::MattesMutualInformationImageToImageMetric<
00102                                 InternalImageType,
00103                                 InternalImageType >           MetricType;
00104 
00105     typedef itk::MeanSquaresImageToImageMetric<
00106                                 InternalImageType,
00107                                 InternalImageType >           MetricTypeMS;
00108 
00109     typedef itk::LinearInterpolateImageFunction<
00110                                 InternalImageType,
00111                                 double >                    InterpolatorType;
00112 
00113     typedef itk::ImageRegistrationMethod<
00114                                 InternalImageType,
00115                                 InternalImageType >           RegistrationType;
00116 
00117     typedef typename itk::WarpImageFilter<
00118                             InternalImageType, 
00119                             InternalImageType,
00120                             DeformationFieldType  >         WarperType;
00121 
00122     typedef typename TransformType::SpacingType                      SpacingType;
00123 
00124     typedef typename TransformType::OriginType                       OriginType;
00125 
00126     typedef itk::ResampleImageFilter< 
00127                                 InternalImageType, 
00128                                 InternalImageType >            ResampleFilterType;
00129 
00130     typedef itk::Image< TPixel, VImageDimension >           OutputImageType;
00131   
00132 
00133     // Sample new image with the same image type as the fixed image
00134     typedef itk::CastImageFilter< 
00135                                 InternalImageType,
00136                                 InternalImageType >            CastFilterType;
00137                     
00138     
00139     typedef itk::Vector< float, VImageDimension >           VectorType;
00140     typedef itk::Image< VectorType, VImageDimension >       DeformationFieldType;
00141 
00142 
00143     typedef itk::BSplineDeformableTransformInitializer <
00144                                 TransformType,
00145                                 InternalImageType >            InitializerType;
00146     
00147 
00148     typename InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
00149     typename RegistrationType::Pointer   registration  = RegistrationType::New();   
00150     typename InitializerType::Pointer    initializer   = InitializerType::New();
00151     typename TransformType::Pointer      transform     = TransformType::New();
00152     
00153     
00154     if(m_Metric==0 || m_Metric==1)
00155     {
00156       typename MetricType::Pointer metric = MetricType::New();
00157       metric->SetNumberOfHistogramBins( 32);
00158       metric->SetNumberOfSpatialSamples(90000);
00159       registration->SetMetric(          metric       );
00160     }
00161     else{
00162       typename MetricTypeMS::Pointer metric = MetricTypeMS::New();  
00163       registration->SetMetric(          metric       );
00164     }
00165      
00166     typename OptimizerFactory::Pointer optFac = OptimizerFactory::New();
00167     optFac->SetOptimizerParameters(m_OptimizerParameters);
00168     optFac->SetNumberOfTransformParameters(transform->GetNumberOfParameters());
00169     OptimizerType::Pointer optimizer = optFac->GetOptimizer();    
00170 
00171     optimizer->AddObserver(itk::AnyEvent(), m_Observer);
00172     
00173     //typedef mitk::MetricFactory <TPixel, VImageDimension> MetricFactoryType;
00174     //typename MetricFactoryType::Pointer metricFac = MetricFactoryType::New();
00175     //metricFac->SetMetricParameters(m_MetricParameters);
00177 
00178 
00179     typename InternalImageType::Pointer fixedImage = InternalImageType::New();
00180     mitk::CastToItkImage(m_ReferenceImage, fixedImage);    
00181     typename InternalImageType::Pointer movingImage = itkImage1;
00182     typename InternalImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
00183     typename InternalImageType::RegionType movingRegion = movingImage->GetBufferedRegion();
00184 
00185     
00186     if(m_MatchHistograms)
00187     {
00188       typedef itk::RescaleIntensityImageFilter<InternalImageType,InternalImageType> FilterType;   
00189       typedef itk::HistogramMatchingImageFilter<InternalImageType,InternalImageType> HEFilterType;
00190 
00191       typename FilterType::Pointer inputRescaleFilter = FilterType::New();  
00192       typename FilterType::Pointer referenceRescaleFilter = FilterType::New();  
00193 
00194       referenceRescaleFilter->SetInput(fixedImage);
00195       inputRescaleFilter->SetInput(movingImage);
00196 
00197       TPixel desiredMinimum =  0;
00198       TPixel desiredMaximum =  255;
00199       
00200       referenceRescaleFilter->SetOutputMinimum( desiredMinimum );
00201       referenceRescaleFilter->SetOutputMaximum( desiredMaximum );
00202       referenceRescaleFilter->UpdateLargestPossibleRegion();  
00203       inputRescaleFilter->SetOutputMinimum( desiredMinimum );
00204       inputRescaleFilter->SetOutputMaximum( desiredMaximum );
00205       inputRescaleFilter->UpdateLargestPossibleRegion();
00206 
00207       // Histogram match the images
00208       typename HEFilterType::Pointer intensityEqualizeFilter = HEFilterType::New();
00209 
00210       intensityEqualizeFilter->SetReferenceImage( inputRescaleFilter->GetOutput() );
00211       intensityEqualizeFilter->SetInput( referenceRescaleFilter->GetOutput() );
00212       intensityEqualizeFilter->SetNumberOfHistogramLevels( 64 );
00213       intensityEqualizeFilter->SetNumberOfMatchPoints( 12 );
00214       intensityEqualizeFilter->ThresholdAtMeanIntensityOn();
00215       intensityEqualizeFilter->Update();
00216 
00217       //fixedImage = referenceRescaleFilter->GetOutput();
00218       //movingImage = IntensityEqualizeFilter->GetOutput();
00219 
00220       fixedImage = intensityEqualizeFilter->GetOutput();
00221       movingImage = inputRescaleFilter->GetOutput();
00222     }
00223 
00224 
00225     //
00226     registration->SetOptimizer(       optimizer     );
00227     registration->SetInterpolator(    interpolator  );  
00228     registration->SetFixedImage(      fixedImage    );
00229     registration->SetMovingImage(     movingImage   );    
00230     registration->SetFixedImageRegion(fixedRegion   );
00231 
00232     initializer->SetTransform(transform);
00233     initializer->SetImage(fixedImage);
00234     initializer->SetNumberOfGridNodesInsideTheImage( m_NumberOfGridPoints );
00235     initializer->InitializeTransform();
00236 
00237     registration->SetTransform( transform );    
00238 
00239     const unsigned int numberOfParameters = transform->GetNumberOfParameters();
00240     
00241     typename itk::BSplineDeformableTransform<
00242                                 double,
00243                                 VImageDimension,
00244                                 3 >::ParametersType  parameters;
00245 
00246     parameters.set_size(numberOfParameters);
00247     parameters.Fill( 0.0 );
00248     transform->SetParameters( parameters );
00249 
00250     // We now pass the parameters of the current transform as the initial
00251     // parameters to be used when the registration process starts.
00252     registration->SetInitialTransformParameters( transform->GetParameters() );
00253     
00254     std::cout << "Intial Parameters = " << std::endl;
00255     std::cout << transform->GetParameters() << std::endl;
00256  
00257 
00258     std::cout << std::endl << "Starting Registration" << std::endl;
00259 
00260     try 
00261     { 
00262       double tstart(clock());     
00263       registration->StartRegistration();    
00264       double time = clock() - tstart;
00265       time = time / CLOCKS_PER_SEC;
00266       MITK_INFO << "Registration time: " << time;
00267     } 
00268     catch( itk::ExceptionObject & err ) 
00269     { 
00270       std::cerr << "ExceptionObject caught !" << std::endl; 
00271       std::cerr << err << std::endl;       
00272     } 
00273     
00274     typename OptimizerType::ParametersType finalParameters = 
00275                       registration->GetLastTransformParameters();
00276 
00277     std::cout << "Last Transform Parameters" << std::endl;
00278     std::cout << finalParameters << std::endl;
00279 
00280     transform->SetParameters( finalParameters );
00281 
00282 /*
00283     ResampleFilterType::Pointer       resampler = ResampleFilterType::New();
00284     resampler->SetTransform(          transform );
00285     resampler->SetInput(              movingImage );
00286     resampler->SetSize(               fixedImage->GetLargestPossibleRegion().GetSize() );
00287     resampler->SetOutputOrigin(       fixedImage->GetOrigin() );
00288     resampler->SetOutputSpacing(      fixedImage->GetSpacing() );
00289     resampler->SetOutputDirection(    fixedImage->GetDirection() );
00290     resampler->SetDefaultPixelValue(  100 );
00291     resampler->SetInterpolator(       interpolator);
00292     resampler->Update();*/
00293 
00294 
00295 
00296     // Generate deformation field
00297     typename DeformationFieldType::Pointer field = DeformationFieldType::New();    
00298     field->SetRegions( movingRegion );
00299     field->SetOrigin( movingImage->GetOrigin() );
00300     field->SetSpacing( movingImage->GetSpacing() );
00301     field->SetDirection( movingImage->GetDirection() );   
00302     field->Allocate();
00303 
00304 
00305     typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;
00306     FieldIterator fi( field, movingRegion );
00307     fi.GoToBegin();
00308 
00309     typename TransformType::InputPointType  fixedPoint;
00310     typename TransformType::OutputPointType movingPoint;
00311     typename DeformationFieldType::IndexType index;
00312 
00313     VectorType displacement;
00314 
00315     while( ! fi.IsAtEnd() )
00316     {
00317       index = fi.GetIndex();
00318       field->TransformIndexToPhysicalPoint( index, fixedPoint );
00319       movingPoint = transform->TransformPoint( fixedPoint );
00320       displacement = movingPoint - fixedPoint;
00321       fi.Set( displacement );
00322       ++fi;
00323     }
00324 
00325 
00326     // Use the deformation field to warp the moving image
00327     typename WarperType::Pointer warper = WarperType::New();    
00328     warper->SetInput( movingImage );
00329     warper->SetInterpolator( interpolator );
00330     warper->SetOutputSpacing( movingImage->GetSpacing() );
00331     warper->SetOutputOrigin( movingImage->GetOrigin() );
00332     warper->SetOutputDirection( movingImage->GetDirection() );
00333     warper->SetDeformationField( field );
00334     warper->Update();
00335 
00336     typename InternalImageType::Pointer result = warper->GetOutput();    
00337 
00338     if(m_UpdateInputImage)
00339     {   
00340       Image::Pointer outputImage = this->GetOutput();
00341       mitk::CastToMitkImage( result, outputImage );
00342     }
00343 
00344 
00345     // Save the deformationfield resulting from the registration    
00346     if(m_SaveDeformationField)
00347     {      
00348       typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
00349       typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
00350 
00351       fieldWriter->SetInput( field );
00352       
00353       fieldWriter->SetFileName( m_DeformationFileName );
00354       try
00355       {
00356         fieldWriter->Update();
00357       }
00358       catch( itk::ExceptionObject & excp )
00359       {
00360         std::cerr << "Exception thrown " << std::endl;
00361         std::cerr << excp << std::endl;
00362         //return EXIT_FAILURE;
00363       }
00364     }
00365 
00366 
00367 
00368   }
00369 } // end namespace