[Insight-users] MRI-SPCET Registration & Transform

[suresh]

Hi Luis,

Thanks for the support you have been extending for me.
I have made all changes suggested by you.(code added at the end of 
this mail)
But still for some reason i could not get proper results.When i 
try extract the pixels in the output image it's failing at 
runtime.

After registration the GetlastTransformParameters returned the 
following values ,
when fixed=MRI, and moving = SPECT

1.000000  0.000000  0.000000
0.000000  1.000000  0.000000
0.000000  0.000000  1.000000
0.000000  0.000000  0.000000

when fixed=SPECT, and moving = MRI, this is result

1.000447  0.002335  0.000873
-0.002983 0.990681  -0.001873
-0.016508  -0.002983 0.996628
0.134443  -0.462841  -0.002983


/************************
REGISTRATION CODE
************************/
 	   double result[16];

 	   typedef BufferToImageConversion<unsigned char> 
ConverterType;
 	   typedef ConverterType::ImageType UCharImage;

    	   ConverterType converter;
 	   UCharImage::Pointer fixedImage = converter.GetImage(MRI);
 	   UCharImage::Pointer movingImage = 
converter.GetImage(SPECT);

 	   typedef itk::AffineTransform<double> TransformType;
 	   typedef itk::RegularStepGradientDescentOptimizer 
OptimizerType;
 	   typedef itk::LinearInterpolateImageFunction <UCharImage, 
double> InterpolatorType;
 	   typedef itk::ImageRegistrationMethod <UCharImage , UCharImage 
 > RegistrationType;
 	   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();
 	   RegistrationType::Pointer  registration  = 
RegistrationType::New();

 	   
/******************************************************************
 	   * Set up the optimizer.
 	   
******************************************************************/

 	   // set the translation scale
 	   typedef OptimizerType::ScalesType ScalesType;
 	   ScalesType parametersScales( 
transform->GetNumberOfParameters() );

 	   parametersScales.Fill( 1.0 );

 	   for (int  j = 9; j < 12; j++ )
 	   {
        	     parametersScales[j] = 0.0001;
 	   }

 	   optimizer->SetScales( parametersScales );

 	   // need to maximize for mutual information
 	   optimizer->MaximizeOn();
 	   
/******************************************************************
 	   * Set up the metric.
            
******************************************************************/
 	   metric->SetMovingImageStandardDeviation( 5.0 );
 	   metric->SetFixedImageStandardDeviation( 5.0 );
 	   metric->SetNumberOfSpatialSamples( 50);
 	   metric->SetFixedImageRegion( fixedImage->GetBufferedRegion() 
);


     
/******************************************************************
      * Set up the registrator.
      
******************************************************************/

 	 // connect up the components
 	  registration->SetMetric( metric );
 	  registration->SetOptimizer( optimizer );
 	  registration->SetTransform( transform );
 	  registration->SetFixedImage(fixedImage );
 	  registration->SetMovingImage( movingImage  );
 	  registration->SetInterpolator( interpolator );

 	  // set initial parameters to identity
 	  RegistrationType::ParametersType initialParameters(
 					transform->GetNumberOfParameters() );

 	  initialParameters.Fill( 0.0 );
 	  initialParameters[0] = 1.0;
 	  initialParameters[4] = 1.0;
 	  initialParameters[8] = 1.0;

       registration->SetInitialTransformParameters( 
initialParameters );
 	  SHOWLINE
       optimizer->SetNumberOfIterations(400);

 	  AfxMessageBox("Start Regsitartion");
 	  try
 	  {
 		  registration->StartRegistration();

 	  }catch(itk::ExceptionObject& Eo)
 	  {
 	     AfxMessageBox(Eo.GetDescription());
 		 return NULL;
 	  }

 	  AfxMessageBox("End Regsitartion");
 	  SHOWLINE
 	  ProgressUpdate(60, "Completed Normalized Registration15", 
MRIVolumeName);

 	  TransformType::ParametersType lastParams = 
registration->GetLastTransformParameters();		  // Get the 
results
 	  for (int pi=0;pi<12;pi++)
 		  result[pi] = lastParams[pi];

 	 /* This Result matrix i'm using to construct the tarnsform in 
Resample code*/

/******************************************
 		TRANSFORM  CODE
******************************************/
 	typedef BufferToImageConversion <unsigned char> ConverterType;
 	typedef ConverterType::ImageType UCharImage;

 	typedef itk::AffineTransform<double, 3> TransformType;

 	typedef itk::ResampleImageFilter<UCharImage, UCharImage>  
ResampleFilter;

 	typedef UCharImage::SizeType  ImageSizeType;
 	typedef TransformType::MatrixType MatrixType;
 	typedef TransformType::OffsetType  VectorType;

 	/*******************
 		initializations
 	********************/
 	ConverterType converter ;
 	ConverterType ::ImagePointer movingImage = 
converter.GetImage(pMoving);
 	ConverterType ::ImagePointer mriImage = 
converter.GetImage(pFixed);
 	ImageSizeType size;
 	size[0]=pTemp->width, size[1]= pTemp->height, size[2] = 
pTemp->depth;

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

 	TransformType::Pointer transform = 
TransformType::New();;//resampleFilter->GetTransform();
 	TransformType::ParametersType 
params(transform->GetNumberOfParameters()) ;

     for (int pi=0;pi<12;pi++)
 		  params[pi] = result[pi];

 	   transform->SetParameters( params);

 	UCharImage::Pointer outImage  = NULL;
 	resampleFilter ->SetTransform(transform);
 	resampleFilter->SetInput(movingImage);

     resampleFilter->SetSize( 
mriImage->GetLargestPossibleRegion().GetSize() );
     resampleFilter->SetOutputOrigin( mriImage->GetOrigin() );
     resampleFilter->SetOutputSpacing( mriImage->GetSpacing() );
 	try
 	{
 		resampleFilter->Update();

 	}catch(itk::ExceptionObject &Eo){
 		AfxMessageBox(Eo.GetDescription());
 	}
 	outImage   = resampleFilter->GetOutput();

/********************************
 	Another important issue iwould like you to observe is my Image 
source code. I have got a buffer
 	from i need to create a itk::Image .This is the code i'm 
using.
 	please Correct if i'm wrong.

********************************/
 		ImagePointer GetImage(Volume* pVol)
 		{
 			PixelType* Buffer = (PixelType*)pVol->Mem;
 			ImagePointer image = ImageType::New();
 			ImageType::SizeType size = {{pVol->width, pVol->height, 
pVol->depth}};
 			ImageType::IndexType index = 
ImageType::IndexType::ZeroIndex;
 			ImageType::RegionType region;

 			region.SetSize(size);
 			region.SetIndex(index);

 			image ->SetLargestPossibleRegion( region );
 			image ->SetBufferedRegion( region );
 			image ->SetRequestedRegion( region );
 			image ->Allocate();

 			itk::ImageRegionIteratorWithIndex <ImageType> it(image , 
region);
 			int k=0;
 			while( !it.IsAtEnd()) {
 				it.Set(Buffer[k]);
 				k++;
 				++it;
 			}
 			return image;
 		}

I could not find out what was wroing in this.What i guess is may 
be i need to modify the Optimizer scales and transform parameters 
values.

Please help me in solving this ..

I Thank youagain  for the help

-suresh