[Insight-users] registration using Mutual Information
Massinissa Bandou
Massinissa.Bandou at USherbrooke.ca
Fri Dec 13 14:59:01 EST 2013
Thank you for replying Mr Luis!!!!
I'm using win7 x64 and I'm compiling on debug mode.
The *source* image is from a series of dicom
format: float
image dimension: 120 x 120 x 137
voxel dimension: 0.5 x 0.5 x 0.596756
The *target* image is a raw file.
format: short
image dimension: 256 x 256 x 256
voxel dimension: 0.308 x 0.308 x 0.308
I decided to register 2 same images and use the same code as
*Registration/ImageRegistration8.cxx*
but I got the following error:
<http://itk-users.7.n7.nabble.com/file/n32985/error.png>
and my code is:
int main(int argc, char *argv[])
{
vtkSmartPointer<vtkImageReader> imagetarget =
vtkSmartPointer<vtkImageReader>::New();
imagetarget->SetFileName("C:/Users/Massi/Desktop/Engineering/IRM_CT_FIDUCIALS/CT_(130802_m_FDGF18_004ct)/image/CT_volume.raw");
imagetarget->SetFileDimensionality(3);
imagetarget->SetDataExtent(0,511,0,511,0,511);
imagetarget->SetDataSpacing(0.154,0.154,0.154);
imagetarget->SetDataByteOrderToLittleEndian();
imagetarget->SetDataScalarTypeToShort();
imagetarget->Update();
vtkSmartPointer<vtkImageReslice> resliceMyTargetImage =
vtkSmartPointer<vtkImageReslice>::New();
resliceMyTargetImage ->SetInputConnection(imagetarget->GetOutputPort());
resliceMyTargetImage->SetOutputExtent(0,255, 0,255, 0,255);
double p[3];
for(unsigned int i=0;i<3;i++){
p[i] = (512 * 0.154) / 256;
}
resliceMyTargetImage->SetOutputSpacing(p[0],p[1],p[2]);
resliceMyTargetImage->SetOutputDimensionality(3);
resliceMyTargetImage->SetInterpolationModeToCubic();
resliceMyTargetImage->AutoCropOutputOn();
resliceMyTargetImage->SetResliceAxesOrigin(0,0,0);
resliceMyTargetImage->Update();
vtkSmartPointer<vtkDICOMImageReader> imagesource =
vtkSmartPointer<vtkDICOMImageReader>::New();
imagesource->SetDirectoryName("C:/Users/Massi/Downloads/M016_Labpet8_20130425/T_6");
imagesource->Update();
typedef itk::Image<float, 3> FixedImageType;
typedef itk::Image<short, 3> MovingImageType;
/*transform vtk image to itk image*/
typedef itk::VTKImageToImageFilter<FixedImageType> FixedImageReaderType;
typedef itk::VTKImageToImageFilter<MovingImageType> MovingImageReaderType;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
fixedImageReader->SetInput(imagesource->GetOutput()); // source
movingImageReader->SetInput(resliceMyTargetImage->GetOutput()); //target
typedef itk::VersorRigid3DTransform< double > TransformType;
typedef itk::VersorRigid3DTransformOptimizer
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 );
fixedImageReader->Update();
registration->SetFixedImageRegion(fixedImageReader->GetOutput()->GetBufferedRegion());
typedef
itk::CenteredTransformInitializer<TransformType,FixedImageType,MovingImageType>
TransformInitializerType;
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
initializer->SetTransform(transform);
initializer->SetFixedImage(fixedImageReader->GetOutput());
initializer->SetMovingImage(movingImageReader->GetOutput());
initializer->MomentsOn();
initializer->InitializeTransform();
typedef TransformType::VersorType VersorType;
typedef VersorType::VectorType VectorType;
VersorType rotation;
VectorType axis;
axis[0] = 0.0;
axis[1] = 0.0;
axis[2] = 1.0;
const double angle = 0;
rotation.Set( axis, angle );
transform->SetRotation( rotation );
registration->SetInitialTransformParameters( transform->GetParameters() );
typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
const double translationScale = 1.0 / 1000.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = translationScale;
optimizerScales[4] = translationScale;
optimizerScales[5] = translationScale;
optimizer->SetScales( optimizerScales );
optimizer->SetMaximumStepLength( 0.2000 );
optimizer->SetMinimumStepLength( 0.0001 );
optimizer->SetNumberOfIterations( 200 );
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
try{
registration->Update();
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;
}
OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
const double versorX = finalParameters[0];
const double versorY = finalParameters[1];
const double versorZ = finalParameters[2];
const double finalTranslationX = finalParameters[3];
const double finalTranslationY = finalParameters[4];
const double finalTranslationZ = finalParameters[5];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
std::cout << std::endl << std::endl;
std::cout << " Result = " << std::endl;
std::cout << " versor X = " << versorX << std::endl;
std::cout << " versor Y = " << versorY << std::endl;
std::cout << " versor Z = " << versorZ << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Translation Z = " << finalTranslationZ << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
transform->SetParameters( finalParameters );
TransformType::MatrixType matrix = transform->GetMatrix();
TransformType::OffsetType offset = transform->GetOffset();
std::cout << "Matrix = " << std::endl << matrix << std::endl;
std::cout << "Offset = " << std::endl << offset << std::endl;
typedef itk::ResampleImageFilter< MovingImageType,FixedImageType >
ResampleFilterType;
TransformType::Pointer finalTransform = TransformType::New();
finalTransform->SetCenter( transform->GetCenter() );
finalTransform->SetParameters( finalParameters );
finalTransform->SetFixedParameters( transform->GetFixedParameters() );
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetTransform( finalTransform );
resampler->SetInput( movingImageReader->GetOutput() );
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resampler->SetOutputOrigin( fixedImage->GetOrigin() );
resampler->SetOutputSpacing( fixedImage->GetSpacing() );
resampler->SetOutputDirection( fixedImage->GetDirection() );
resampler->SetDefaultPixelValue( 100 );
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, 3 > OutputImageType;
typedef itk::CastImageFilter< FixedImageType, OutputImageType >
CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
//WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
//writer->SetFileName( "result.raw" );
caster->SetInput( resampler->GetOutput() );
//writer->SetInput( caster->GetOutput() );
//writer->Update();
typedef
itk::SubtractImageFilter<FixedImageType,FixedImageType,FixedImageType >
DifferenceFilterType;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
typedef itk::RescaleIntensityImageFilter<FixedImageType,OutputImageType >
RescalerType;
RescalerType::Pointer intensityRescaler = RescalerType::New();
intensityRescaler->SetInput( difference->GetOutput() );
intensityRescaler->SetOutputMinimum( 0 );
intensityRescaler->SetOutputMaximum( 255 );
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resampler->GetOutput() );
resampler->SetDefaultPixelValue( 1 );
/*connect to VTK*/
typedef itk::ImageToVTKImageFilter<OutputImageType> ConnectorType;
ConnectorType::Pointer connector = ConnectorType::New();
connector->SetInput(caster->GetOutput());
connector->Update();
vtkSmartPointer<vtkGPUVolumeRayCastMapper> volumeMapper =
vtkSmartPointer<vtkGPUVolumeRayCastMapper>::New();
volumeMapper->SetInput(connector->GetOutput());
volumeMapper->AutoAdjustSampleDistancesOff();
volumeMapper->SetBlendModeToComposite();
/********************set color & opacity*******************/
vtkSmartPointer<vtkVolumeProperty> volumeProperty =
vtkSmartPointer<vtkVolumeProperty>::New();
vtkSmartPointer<vtkPiecewiseFunction> compositeOpacity =
vtkSmartPointer<vtkPiecewiseFunction>::New();
vtkSmartPointer<vtkColorTransferFunction> color =
vtkSmartPointer<vtkColorTransferFunction>::New();
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[0],0.0);
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*3/255,0.0);
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*15/255,0.1);
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*70/255,0.2);
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*75/255,0.3);
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*0/255,0.0,0.0,0.0);
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*3/255,1.0,0.3,0.3);//
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*20/255,0.0,0.0,1.0);
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*50/255,0.9,0.37,0.27);
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*80.0/255,0.15,0.15,0.5);
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*120.0/255,1.0,1.0,1.0);
volumeProperty->SetAmbient(0.5);
volumeProperty->SetDiffuse(0.8);
volumeProperty->SetSpecular(8.0);
volumeProperty->SetInterpolationTypeToLinear();
volumeProperty->SetColor(color);
volumeProperty->SetScalarOpacity(compositeOpacity);
volumeProperty->ShadeOn();
vtkSmartPointer<vtkVolume> volume = vtkSmartPointer<vtkVolume>::New();
volume->SetMapper(volumeMapper);
volume->SetProperty(volumeProperty);
vtkSmartPointer<vtkRenderWindow> renderWindow =
vtkSmartPointer<vtkRenderWindow>::New();
vtkSmartPointer<vtkRenderer> renderer =
vtkSmartPointer<vtkRenderer>::New();
renderWindow->AddRenderer(renderer);
renderer->AddVolume(volume);
renderer->SetBackground(0.9,0.9,0.9);
renderer->GradientBackgroundOn();
vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
vtkSmartPointer<vtkRenderWindowInteractor>::New();
renderWindowInteractor->SetRenderWindow(renderWindow);
/*Set camera position*/
renderer->GetActiveCamera()->Azimuth(30);
renderer->GetActiveCamera()->Elevation(30);
renderer->ResetCamera();
renderWindow->Render();
renderWindowInteractor->Start();
return EXIT_SUCCESS;
}
--
View this message in context: http://itk-users.7.n7.nabble.com/registration-using-Mutual-Information-tp32972p32985.html
Sent from the ITK - Users mailing list archive at Nabble.com.
More information about the Insight-users
mailing list