[ITK Community] [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.
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-users



More information about the Community mailing list