0
2つの異なるフォルダにある2つの異なるDICOMイメージシリーズを読み込み、登録を実行しようとしています。このシリーズは正しく読み込まれているようですが、登録 - > Update（）が呼び出されるまでスムーズに進みます。それは直接破壊され、中止関数がUpdate（）の内部で呼び出される可能性があります。プログラムは2D画像でうまく機能していました。ここにコードがありますITK DICOMシリーズ登録
#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登録の例を教えてもらえると便利です。
ありがとうございます。
32ビットウィンドウで実行しているか、32ビットでコンパイルしていますか？通常、32ビットプログラムは2X 3D DICOM +フロートイメージの場合に素早く使い果たされる可能性がある（RAMまたはスワップの量に関係なく）2GBの断片化されたアドレス空間に制限されるためです。 – drescherjm