我正在尝试在两个不同的文件夹中读取两组不同的 DICOM 图像系列并执行注册。该系列似乎被正确阅读，一切顺利，直到注册->更新（）被调用。它直接崩溃，很可能在 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 注册示例，那也会很有帮助。
提前致谢。