[Insight-users] Transforming a Volume

[suresh]

Hi All,

Tahnsks Luis.Thanks Sayan.. With your help i could complete the 
segmentation using NeighborhoodConnectedImageFilter.

Now i'm facing aproblem in tarnsforming. My problem is ..

I have an MRI and SPECT scan of different sizes. I'm registering 
SPECT to MRI using the code listed at the end.
Then i'm usingthe output matrix to tarnsform the SPECT . What i'm 
excpecting is the SPECT to GROW to the MRI volume size.
I tried with ResampleImageFilter but i could not get it to work.

Please help me with the following issues..
1. Is my Registration code correct.?Or do i need change any 
parameters.
2. How do i apply a transformation matrix to a vaolume.?

____________________________Registration Code 
___________________________________________________
//******This code is written based on itk MIRegistartion 
example.***/

 	   typedef itk::QuaternionRigidTransform<double> 
TransformType;
 	   typedef itk::QuaternionRigidTransformGradientDescentOptimizer 
OptimizerType;
 	   typedef itk::LinearInterpolateImageFunction <UCharImage, 
double> InterpolatorType;
 	   typedef itk::ImageRegistrationMethod <UCharImage , UCharImage 
 > RegistrationType1;

 	   typedef itk::MutualInformationImageToImageMetric<UCharImage , 
UCharImage> MetricType;
 	   ProgressUpdate(50, "Performing MutualInformation 
Registration", MRIVolumeName);
 	   MetricType::Pointer         metric        = 
MetricType::New();
 	   TransformType::Pointer      transform     = 
TransformType::New();
 	   OptimizerType::Pointer      optimizer     = 
OptimizerType::New();
 	   InterpolatorType::Pointer   interpolator  = 
InterpolatorType::New();
 	   RegistrationType1::Pointer  registration  = 
RegistrationType1::New();

 	   ConverterType converter; // Image source
 	   UCharImage::Pointer fixedImage = converter.GetImage(MRI);// 
converter is my Image Source
 	   UCharImage::Pointer movingImage = 
converter.GetImage(SPECT);

 	   RegistrationType1::ParametersType 
guess(transform->GetNumberOfParameters() );

//******What is this guess for????//
 	   guess[0] = 0.0;	guess[1] = 0.0;	 guess[2] = 0.0; guess[3] = 
1.0;
 	   guess[4] = 20.0; guess[5] = 40.0; guess[6] = 0.0;
 	   try{
 		   registration->SetInitialTransformParameters (guess);
 		   metric->SetMovingImageStandardDeviation( 2 );
 		   metric->SetFixedImageStandardDeviation( 2 );
 		   metric->SetNumberOfSpatialSamples( 1000 );
 		   typedef OptimizerType::ScalesType ScaleType;
 		   ScaleType scales(transform->GetNumberOfParameters());
 		   scales.Fill( 1.0 );
 		   for( unsigned j = 4; j < 7; j++ )
 		   {
 			   scales[j] = 1.0 / vnl_math_sqr(300.0);
 		   }
 		   optimizer->SetNumberOfIterations( 7 );

//******What is this learning rate??
 		   optimizer->SetLearningRate( 0.0000001 );

 		   registration->SetMetric(metric);
 		   registration->SetOptimizer(optimizer);
 		   registration->SetTransform(transform);
 		   registration->SetInterpolator(interpolator);
 		   registration->SetFixedImage(fixedImage);
 		   registration->SetMovingImage(movingImage);
 		   optimizer->SetScales(scales);
 		   registration->StartRegistration();
 		   ProgressUpdate(60, "Performing Normalized Registration15", 
MRIVolumeName);
 		   itk::Array<double>  *arr=new itk::Array<double>;
 		   RegistrationType1::ParametersType solution 		 
=registration->GetLastTransformParameters();
 		   vnl_quaternion<double> 
quat(solution[0],solution[1],solution[2],solution[3]);
 		   vnl_matrix_fixed<double,3,3> vnlmat = 
quat.rotation_matrix();

 		   result[0] = vnlmat[0][0]; result[1] = vnlmat[0][1];
 		   result[2] = vnlmat[0][2]; result[3] = solution[4];
 		   result[4] = vnlmat[1][0]; result[5] = vnlmat[1][1];
 		   result[6] = vnlmat[1][2]; result[7] = solution[5];
 		   result[8] = vnlmat[2][0]; result[9] = vnlmat[2][1];
 		   result[10] = vnlmat[2][2]; result[11] = solution[6];
 		   result[12] = 0;	result[13] = 0;	result[14] = 0; result[15] = 
1;
_______________________________________________________________________________________________
***************Transform Code*********************
 	typedef itk::ResampleImageFilter<UCharImage, UCharImage>  
ResampleFilter;
 	typedef ResampleFilter::TransformType  TransformType  ;
 	typedef UCharImage::SizeType  ImageSizeType;
 	typedef  TransformType::MatrixType MatrixType;
 	typedef  TransformType::OffsetType  VectorType;
 	ConverterType converter ;
 	ConverterType ::ImagePointer image = converter.GetImage(pVol);
 	ImageSizeType size;
 	size[0]=pTemp->width, size[1]= pTemp->height, size[2] = 
pTemp->depth;

 	ResampleFilter::Pointer resampleFilter = 
ResampleFilter::New();

 	TransformType::Pointer transform = 
resampleFilter->GetTransform();
 	VectorType vector2 = transform->GetOffset();
 	MatrixType matrix2 = transform->GetMatrix();

// data is my Registration output matrix of size 4X4
 	matrix2[0][0] = data[0][0];  matrix2[0][1] = matrix2[0][1] ;  
matrix2[0][2] = matrix2[0][2] ;
 	matrix2[1][0] = data[1][0];  matrix2[1][1] = matrix2[1][1] ;  
matrix2[1][2] = matrix2[1][2] ;
 	matrix2[2][0] = data[2][0];  matrix2[2][1] = matrix2[2][1] ;  
matrix2[2][2] = matrix2[2][2] ;
 	transform->SetOffset(vector2);
 	transform->SetMatrix(matrix2);
 	UCharImage::Pointer outImage  = NULL;
 	resampleFilter ->SetTransform(transform);
     resampleFilter->SetInput(image);
 	try{
 		resampleFilter->Update();
 	}catch(itk::ExceptionObject &Eo){
 		Eo.GetDescription();
 	}
 	outImage   = resampleFilter->GetOutput();
_______________________________________________________________________________________________

thanks,

suresh