mitk: mitkDemonsRegistration.cpp Source File

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003 Program:   Medical Imaging & Interaction Toolkit
00004 Language:  C++
00005 Date:      $Date$
00006 Version:   $Revision$
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 "mitkDemonsRegistration.h"
00023 
00024 namespace mitk {
00025 
00026   DemonsRegistration::DemonsRegistration():
00027     m_Iterations(50),
00028     m_StandardDeviation(1.0),
00029     m_FieldName("newField.mhd"),
00030     m_ResultName("deformedImage.mhd"),
00031     m_SaveField(true),
00032     m_SaveResult(true),
00033     m_DeformationField(NULL)
00034   {
00035 
00036   }
00037 
00038   DemonsRegistration::~DemonsRegistration()
00039   {
00040   }
00041 
00042   void DemonsRegistration::SetNumberOfIterations(int iterations)
00043   {
00044     m_Iterations = iterations;
00045   }
00046 
00047   void DemonsRegistration::SetStandardDeviation(float deviation)
00048   {
00049     m_StandardDeviation = deviation;
00050   }
00051 
00052   void DemonsRegistration::SetSaveDeformationField(bool saveField)
00053   {
00054     m_SaveField = saveField;
00055   }
00056 
00057   void DemonsRegistration::SetDeformationFieldFileName(const char* fieldName)
00058   {
00059     m_FieldName = fieldName;
00060   }
00061 
00062   void DemonsRegistration::SetSaveResult(bool saveResult)
00063   {
00064     m_SaveResult = saveResult;
00065   }
00066 
00067   void DemonsRegistration::SetResultFileName(const char* resultName)
00068   {
00069     m_ResultName = resultName;
00070   }
00071 
00072   itk::Image<itk::Vector<float, 3>,3>::Pointer DemonsRegistration::GetDeformationField()
00073   {
00074     return m_DeformationField;
00075   }
00076 
00077   template < typename TPixel, unsigned int VImageDimension >
00078     void DemonsRegistration::GenerateData2( itk::Image<TPixel, VImageDimension>* itkImage1)
00079   {
00080     typedef typename itk::Image< TPixel, VImageDimension >  FixedImageType;
00081     typedef typename itk::Image< TPixel, VImageDimension >  MovingImageType;
00082 
00083     typedef float InternalPixelType;
00084     typedef typename itk::Image< InternalPixelType, VImageDimension > InternalImageType;
00085     typedef typename itk::CastImageFilter< FixedImageType, 
00086                                   InternalImageType > FixedImageCasterType;
00087     typedef typename itk::CastImageFilter< MovingImageType, 
00088                                   InternalImageType > MovingImageCasterType;
00089     typedef typename itk::Image< InternalPixelType, VImageDimension > InternalImageType;
00090     typedef typename itk::Vector< float, VImageDimension >    VectorPixelType;
00091     typedef typename itk::Image<  VectorPixelType, VImageDimension > DeformationFieldType;
00092     typedef typename itk::DemonsRegistrationFilter<
00093                                   InternalImageType,
00094                                   InternalImageType,
00095                                   DeformationFieldType>   RegistrationFilterType;
00096     typedef typename itk::WarpImageFilter<
00097                             MovingImageType, 
00098                             MovingImageType,
00099                             DeformationFieldType  >     WarperType;
00100     typedef typename itk::LinearInterpolateImageFunction<
00101                                     MovingImageType,
00102                                     double          >  InterpolatorType;
00103 
00104     typedef  TPixel  OutputPixelType;
00105     typedef typename itk::Image< OutputPixelType, VImageDimension > OutputImageType;
00106     typedef typename itk::CastImageFilter< 
00107                           MovingImageType,
00108                           OutputImageType > CastFilterType;
00109     typedef typename itk::ImageFileWriter< OutputImageType >  WriterType;
00110     typedef typename itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
00111 
00112     typename FixedImageType::Pointer fixedImage = FixedImageType::New();
00113     mitk::CastToItkImage(m_ReferenceImage, fixedImage);
00114     typename MovingImageType::Pointer movingImage = itkImage1;
00115 
00116     if (fixedImage.IsNotNull() && movingImage.IsNotNull())
00117     {
00118       typename RegistrationFilterType::Pointer filter = RegistrationFilterType::New();
00119 
00120       this->AddStepsToDo(4);
00121       typename itk::ReceptorMemberCommand<DemonsRegistration>::Pointer command = itk::ReceptorMemberCommand<DemonsRegistration>::New();
00122       command->SetCallbackFunction(this, &DemonsRegistration::SetProgress);
00123       filter->AddObserver( itk::IterationEvent(), command );
00124 
00125       typename FixedImageCasterType::Pointer fixedImageCaster = FixedImageCasterType::New();
00126       fixedImageCaster->SetInput(fixedImage);
00127       filter->SetFixedImage( fixedImageCaster->GetOutput() );
00128       typename MovingImageCasterType::Pointer movingImageCaster = MovingImageCasterType::New();
00129       movingImageCaster->SetInput(movingImage);
00130       filter->SetMovingImage(movingImageCaster->GetOutput());
00131       filter->SetNumberOfIterations( m_Iterations );
00132       filter->SetStandardDeviations( m_StandardDeviation );
00133       filter->Update();
00134 
00135       typename WarperType::Pointer warper = WarperType::New();
00136       typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
00137 
00138       warper->SetInput( movingImage );
00139       warper->SetInterpolator( interpolator );
00140       warper->SetOutputSpacing( fixedImage->GetSpacing() );
00141       warper->SetOutputOrigin( fixedImage->GetOrigin() );
00142       warper->SetOutputDirection( fixedImage->GetDirection());
00143       warper->SetDeformationField( filter->GetOutput() );
00144       warper->Update();
00145       Image::Pointer outputImage = this->GetOutput();
00146       mitk::CastToMitkImage( warper->GetOutput(), outputImage );
00147 
00148 
00149       typename WriterType::Pointer      writer =  WriterType::New();
00150       typename CastFilterType::Pointer  caster =  CastFilterType::New(); 
00151 
00152       writer->SetFileName( m_ResultName );
00153   
00154       caster->SetInput( warper->GetOutput() );
00155       writer->SetInput( caster->GetOutput()   );
00156       if(m_SaveResult)
00157       {
00158         writer->Update();
00159       }
00160 
00161       if (VImageDimension == 2)
00162       {
00163         typedef DeformationFieldType  VectorImage2DType;
00164         typedef typename DeformationFieldType::PixelType Vector2DType;
00165 
00166         typename VectorImage2DType::ConstPointer vectorImage2D = filter->GetOutput();
00167 
00168         typename VectorImage2DType::RegionType  region2D = vectorImage2D->GetBufferedRegion();
00169         typename VectorImage2DType::IndexType   index2D  = region2D.GetIndex();
00170         typename VectorImage2DType::SizeType    size2D   = region2D.GetSize(); 
00171 
00172 
00173         typedef typename itk::Vector< float,       3 >  Vector3DType;
00174         typedef typename itk::Image< Vector3DType, 3 >  VectorImage3DType;
00175 
00176         typedef typename itk::ImageFileWriter< VectorImage3DType > WriterType;
00177 
00178         typename WriterType::Pointer writer3D = WriterType::New();
00179 
00180         typename VectorImage3DType::Pointer vectorImage3D = VectorImage3DType::New();
00181         
00182         typename VectorImage3DType::RegionType  region3D;
00183         typename VectorImage3DType::IndexType   index3D;
00184         typename VectorImage3DType::SizeType    size3D;
00185 
00186         index3D[0] = index2D[0];
00187         index3D[1] = index2D[1];
00188         index3D[2] = 0;
00189         
00190         size3D[0]  = size2D[0];
00191         size3D[1]  = size2D[1];
00192         size3D[2]  = 1;
00193 
00194         region3D.SetSize( size3D );
00195         region3D.SetIndex( index3D );
00196 
00197         typename VectorImage2DType::SpacingType spacing2D   = vectorImage2D->GetSpacing();
00198         typename VectorImage3DType::SpacingType spacing3D;
00199 
00200         spacing3D[0] = spacing2D[0];
00201         spacing3D[1] = spacing2D[1];
00202         spacing3D[2] = 1.0;
00203 
00204         vectorImage3D->SetSpacing( spacing3D );
00205 
00206         vectorImage3D->SetRegions( region3D );
00207         vectorImage3D->Allocate();
00208         
00209         typedef typename itk::ImageRegionConstIterator< VectorImage2DType > Iterator2DType;
00210 
00211         typedef typename itk::ImageRegionIterator< VectorImage3DType > Iterator3DType;
00212 
00213         Iterator2DType  it2( vectorImage2D, region2D );
00214         Iterator3DType  it3( vectorImage3D, region3D );
00215 
00216         it2.GoToBegin();
00217         it3.GoToBegin();
00218 
00219         Vector2DType vector2D;
00220         Vector3DType vector3D;
00221 
00222         vector3D[2] = 0; // set Z component to zero.
00223 
00224         while( !it2.IsAtEnd() )
00225         {
00226           vector2D = it2.Get();
00227           vector3D[0] = vector2D[0];  
00228           vector3D[1] = vector2D[1];  
00229           it3.Set( vector3D );
00230           ++it2;
00231           ++it3;
00232         }
00233 
00234         writer3D->SetInput( vectorImage3D );
00235         m_DeformationField = vectorImage3D;
00236 
00237         writer3D->SetFileName( m_FieldName );
00238 
00239         try
00240         {
00241           if(m_SaveField)
00242           {
00243             writer3D->Update();
00244           }
00245         }
00246         catch( itk::ExceptionObject & excp )
00247         {
00248           MITK_ERROR << excp << std::endl;
00249         }
00250       }
00251       else
00252       {
00253         typename FieldWriterType::Pointer      fieldwriter =  FieldWriterType::New();
00254         fieldwriter->SetFileName(m_FieldName);
00255         fieldwriter->SetInput( filter->GetOutput() );
00256         m_DeformationField = (itk::Image<itk::Vector<float, 3>,3> *)(filter->GetOutput());
00257         if(m_SaveField)
00258         {
00259           fieldwriter->Update();
00260         }
00261       }
00262       this->SetRemainingProgress(4);
00263     }
00264   }
00265 } // end namespace