ITK DICOM series регистрация

Question

Я пытаюсь читать в двух разных наборах серии изображений DICOM в двух разных папках и выполнять регистрацию. Кажется, что строка читается правильно, все идет гладко до тех пор, пока не вызывается регистрация -> Update(). Он подавляется напрямую и, вероятно, функция abort вызывается внутри Update(). Программа отлично работала с 2D-изображениями. Вот код

#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 

Я боролся с этим в течение многих недель и до сих пор не мог понять, в чем проблема. Я попытался найти в руководстве пользователя и примерах, но ни один из них не читает серии изображений DICOM.

Было бы также полезно, если бы кто-нибудь мог предоставить мне пример регистрации ITK в серии изображений.

Заранее спасибо.

Answer 1

Если вы не уверены, что серии DICOM читаются правильно, возможно, вам стоит попробовать сначала проверить эту часть кода. Я использовал подобный пример без проблем. Тем не менее, я читал в обеих сериях DICOM и писал их в 3D-набор данных, а затем выполнил регистрацию. Это может помочь.