问ITK DICOM系列注册

EN

   #include "itkImageRegistrationMethod.h"
   #include "itkMeanSquaresImageToImageMetric.h"

   #include "itkTimeProbesCollectorBase.h"
   #include "itkMemoryProbesCollectorBase.h"

   #include "itkBSplineTransform.h"
   #include "itkLBFGSBOptimizer.h"

   #include "itkImageFileReader.h"
   #include "itkImageFileWriter.h"

   #include "itkResampleImageFilter.h"
   #include "itkCastImageFilter.h"
   #include "itkSquaredDifferenceImageFilter.h"

   #include "itkGDCMImageIO.h"
   #include "itkGDCMSeriesFileNames.h"
   #include "itkNumericSeriesFileNames.h"
   #include "gdcmUIDGenerator.h"

   #include "itkImageSeriesReader.h"
   #include "itkImageSeriesWriter.h"

   #define SRI

   #include "itkCommand.h"
   class CommandIterationUpdate : public itk::Command
   {
   public:
     typedef  CommandIterationUpdate   Self;
     typedef  itk::Command             Superclass;
     typedef itk::SmartPointer<Self>   Pointer;
     itkNewMacro( Self );
   protected:
     CommandIterationUpdate() {};
   public:
     typedef itk::LBFGSBOptimizer       OptimizerType;
     typedef   const OptimizerType *    OptimizerPointer;

     void Execute(itk::Object *caller, const itk::EventObject & event)
       {
       Execute( (const itk::Object *)caller, event);
       }

     void Execute(const itk::Object * object, const itk::EventObject & event)
       {
       OptimizerPointer optimizer =
         dynamic_cast< OptimizerPointer >( object );
       if( !(itk::IterationEvent().CheckEvent( &event )) )
         {
         return;
         }
       std::cout << optimizer->GetCurrentIteration() << "   ";
       std::cout << optimizer->GetValue() << "   ";
       std::cout << optimizer->GetInfinityNormOfProjectedGradient() << std::endl;
       }
   };


   int main( int argc, char *argv[] )
   {
     if( argc < 4 )
       {
       std::cerr << "Missing Parameters " << std::endl;
       std::cerr << "Usage: " << argv[0];
       std::cerr << " fixedImageFile  movingImageFile outputImagefile  ";
       std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";
       std::cerr << " [deformationField] ";
       return EXIT_FAILURE;
       }

     const    unsigned int    ImageDimension = 3;
     typedef  float           PixelType;

     typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
     typedef itk::Image< PixelType, ImageDimension >  MovingImageType;


     const unsigned int SpaceDimension = ImageDimension;
     const unsigned int SplineOrder = 3;
     typedef double CoordinateRepType;

     typedef itk::BSplineTransform<
                               CoordinateRepType,
                               SpaceDimension,
                               SplineOrder >     TransformType;

     typedef itk::LBFGSBOptimizer       OptimizerType;

     typedef itk::MeanSquaresImageToImageMetric<
                                       FixedImageType,
                                       MovingImageType >    MetricType;

     typedef itk:: LinearInterpolateImageFunction<
                                       MovingImageType,
                                       double          >    InterpolatorType;

     typedef itk::ImageRegistrationMethod<
                                       FixedImageType,
                                       MovingImageType >    RegistrationType;

     MetricType::Pointer         metric        = MetricType::New();
     OptimizerType::Pointer      optimizer     = OptimizerType::New();
     InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
     RegistrationType::Pointer   registration  = RegistrationType::New();

     registration->SetMetric(        metric        );
     registration->SetOptimizer(     optimizer     );
     registration->SetInterpolator(  interpolator  );

     TransformType::Pointer  transform = TransformType::New();
     registration->SetTransform( transform );

   #ifdef SRI
     typedef itk::ImageSeriesReader< FixedImageType > FixedImageReaderType;
     typedef itk::ImageSeriesReader< MovingImageType > MovingImageReaderType;

     typedef itk::GDCMImageIO ImageIOType;
     typedef itk::GDCMSeriesFileNames InputNamesGeneratorType;
     typedef itk::NumericSeriesFileNames OutputNameGeneratorType;
     typedef itk::ImageSeriesWriter< MovingImageType, MovingImageType >             SeriesWriterType;

     ImageIOType::Pointer gdcmIO = ImageIOType::New();

     InputNamesGeneratorType::Pointer fixedImageNames = InputNamesGeneratorType::New();
     InputNamesGeneratorType::Pointer movingImageNames =        InputNamesGeneratorType::New();

     fixedImageNames->SetInputDirectory( argv[1] );
     movingImageNames->SetInputDirectory( argv[2] );

     const FixedImageReaderType::FileNamesContainer & fixedNames = fixedImageNames-       >GetInputFileNames();
     const MovingImageReaderType::FileNamesContainer & movingNames = movingImageNames->GetInputFileNames();

   #else
     typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
     typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
   #endif

     FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
     MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();

   #ifdef SRI
     fixedImageReader->SetImageIO( gdcmIO );
     movingImageReader->SetImageIO( gdcmIO );

     fixedImageReader->SetFileNames( fixedNames );
     movingImageReader->SetFileNames( movingNames );

   #else
     fixedImageReader->SetFileName(  argv[1] );
     movingImageReader->SetFileName( argv[2] );
   #endif

     FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();

     registration->SetFixedImage(  fixedImage   );
     registration->SetMovingImage(   movingImageReader->GetOutput()   );

     fixedImageReader->Update();

     FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();

    registration->SetFixedImageRegion( fixedRegion );

     typedef TransformType::RegionType RegionType;

     unsigned int numberOfGridNodes = 8;

     TransformType::PhysicalDimensionsType   fixedPhysicalDimensions;
     TransformType::MeshSizeType             meshSize;
     TransformType::OriginType               fixedOrigin;

     for( unsigned int i=0; i< SpaceDimension; i++ )
       {
       fixedOrigin = fixedImage->GetOrigin()[i];
       fixedPhysicalDimensions[i] = fixedImage->GetSpacing()[i] *
         static_cast<double>(
         fixedImage->GetLargestPossibleRegion().GetSize()[i] - 1 );
       }
     meshSize.Fill( numberOfGridNodes - SplineOrder );

     transform->SetTransformDomainOrigin( fixedOrigin );
     transform->SetTransformDomainPhysicalDimensions(
       fixedPhysicalDimensions );
     transform->SetTransformDomainMeshSize( meshSize );
     transform->SetTransformDomainDirection( fixedImage->GetDirection() );

     typedef TransformType::ParametersType     ParametersType;

     const unsigned int numberOfParameters =
                  transform->GetNumberOfParameters();

     ParametersType parameters( numberOfParameters );

     parameters.Fill( 0.0 );

     transform->SetParameters( parameters );

     registration->SetInitialTransformParameters( transform->GetParameters() );

     std::cout << "Intial Parameters = " << std::endl;
     std::cout << transform->GetParameters() << std::endl;

     OptimizerType::BoundSelectionType boundSelect( transform->GetNumberOfParameters()        );
     OptimizerType::BoundValueType upperBound( transform->GetNumberOfParameters() );
     OptimizerType::BoundValueType lowerBound( transform->GetNumberOfParameters() );

     boundSelect.Fill( 0 );
     upperBound.Fill( 0.0 );
     lowerBound.Fill( 0.0 );

     optimizer->SetBoundSelection( boundSelect );
     optimizer->SetUpperBound( upperBound );
     optimizer->SetLowerBound( lowerBound );

     optimizer->SetCostFunctionConvergenceFactor( 1e+12 );
     optimizer->SetProjectedGradientTolerance( 1.0 );
     optimizer->SetMaximumNumberOfIterations( 500 );
     optimizer->SetMaximumNumberOfEvaluations( 500 );
     optimizer->SetMaximumNumberOfCorrections( 5 );

     CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
     optimizer->AddObserver( itk::IterationEvent(), observer );

     itk::TimeProbesCollectorBase chronometer;
     itk::MemoryProbesCollectorBase memorymeter;

     std::cout << std::endl << "Starting Registration" << std::endl;

     try
       {
       memorymeter.Start( "Registration" );
       chronometer.Start( "Registration" );

       registration->Update();

       chronometer.Stop( "Registration" );
       memorymeter.Stop( "Registration" );

       std::cout << "Optimizer stop condition = "
                 << registration->GetOptimizer()->GetStopConditionDescription()
                 << std::endl;
       }
     catch( itk::ExceptionObject & err )
       {
       std::cerr << "ExceptionObject caught !" << std::endl;
       std::cerr << err << std::endl;
       return EXIT_FAILURE;
       }

     OptimizerType::ParametersType finalParameters =
                       registration->GetLastTransformParameters();

     std::cout << "Last Transform Parameters" << std::endl;
     std::cout << finalParameters << std::endl;


     // Report the time taken by the registration
     chronometer.Report( std::cout );
     memorymeter.Report( std::cout );

     transform->SetParameters( finalParameters );

     // resample and output