[Insight-users] 3D bspline deformable registration works for example data but gives exception for breast mri data

Sachin Jambawalikar sachinjam at gmail.com
Mon Oct 4 11:41:19 EDT 2004


Hi Luis ,
The only thing needed to modify for deformable registration using
Bspline for  3D  breast MRI data   was image dimension.

this code works weel for Brain web simulated data  but  gives me
following error for my data sets (
http://www.ic.sunysb.edu/Stu/sjambawa/b_images/volumes.tar.gz )

  I   NFN    FUNC        GNORM       STEPLENGTH
 IFLAG= -1. LINE SEARCH FAILED. SEE DOCUMENTATION OF ROUTINE MCSRCH ERROR RETURN
  OF LINE SEARCH: INFO=3 POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE  INCORRECT OR
 INCORRECT TOLERANCESvnl_lbfgs: ** EEEK **
 ExceptionObject caught !

 itk::ExceptionObject (00F2FB80)
 Location: "Unknown"
 File: g:\itk_vtk\itk\insighttoolkit-1.6.0\code\common\itkBSplineDeformableTransf
 orm.txx
 Line: 232
 Description: itk::ERROR: BSplineDeformableTransform(003337D0):
Mismatched betwee
 n parameters size and region size

Breast MRI  data=512 x 512 x 60 (0.68x0.68x3)
BrainWeb Data   =181x217x180(1x1x1)

=======================
Code  for deformable registration
========================

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 1;
		}
		*/
		const    unsigned int    ImageDimension = 3;
		typedef  short           PixelType;

		typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
		typedef itk::Image< PixelType, ImageDimension >  MovingImageType;


		//  Software Guide : BeginLatex
		//
		//  We instantiate now the type of the \code{BSplineDeformableTransform} using
		//  as template parameters the type for coordinates representation, the
		//  dimension of the space, and the order of the BSpline. 
		// 
		//  \index{BSplineDeformableTransform|New}
		//  \index{BSplineDeformableTransform|Instantiation}
		//
		//  Software Guide : EndLatex 

		// Software Guide : BeginCodeSnippet
		const unsigned int SpaceDimension = ImageDimension;
		const unsigned int SplineOrder = 3;
		typedef double CoordinateRepType;

		typedef itk::BSplineDeformableTransform<
			CoordinateRepType,
			ImageDimension,
			SplineOrder >     TransformType;
		// Software Guide : EndCodeSnippet


		typedef itk::LBFGSOptimizer       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  );


		//  Software Guide : BeginLatex
		//
		//  The transform object is constructed below and passed to the registration
		//  method.
		//  \index{itk::RegistrationMethod!SetTransform()}
		//
		//  Software Guide : EndLatex 

		// Software Guide : BeginCodeSnippet
		TransformType::Pointer  transform = TransformType::New();
		registration->SetTransform( transform );
		// Software Guide : EndCodeSnippet

		typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
		typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;

		FixedImageReaderType::Pointer  fixedImageReader  =
			FixedImageReaderType::New();
		MovingImageReaderType::Pointer movingImageReader =
			MovingImageReaderType::New();

		//itk::DICOMImageIO2Factory::RegisterOneFactory();
		/*FILE *fp;
		fp=fopen("test.dat","w");
		fclose(fp); */

		std::string a= "Image1.hdr";
		std::string b="Image2.hdr";
		
		fixedImageReader->SetFileName(a.c_str());
		movingImageReader->SetFileName(b.c_str());

		FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
		MovingImageType::Pointer movingImage = movingImageReader->GetOutput();




		registration->SetFixedImage(  fixedImage   );
		registration->SetMovingImage(   movingImageReader->GetOutput()   );

		fixedImageReader->Update();

		FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();

		registration->SetFixedImageRegion( fixedRegion );

		movingImageReader->Update();
		std::cout << fixedImage->GetBufferedRegion()<<std::endl;
		std::cout << movingImage->GetBufferedRegion()<<std::endl;

		//  Software Guide : BeginLatex
		//
		//  Here we define the parameters of the BSplineDeformableTransform grid.  We
		//  arbitrarily decide to use a grid with $5 \times 5$ nodes within the image.
			//  The reader should note that the BSpline computation requires a
			//  finite support region ( 1 grid node at the lower borders and 2
			//  grid nodes at upper borders). Therefore in this example, we set
			//  the grid size to be $8 \times 8$ and place the grid origin such that
			//  grid node (1,1) coinicides with the first pixel in the fixed image.
			// 
			//  \index{BSplineDeformableTransform}
			//
			//  Software Guide : EndLatex 


			// Software Guide : BeginCodeSnippet
			typedef TransformType::RegionType RegionType;
		RegionType bsplineRegion;
		RegionType::SizeType   gridSizeOnImage;
		RegionType::SizeType   gridBorderSize;
		RegionType::SizeType   totalGridSize;

		gridSizeOnImage.Fill( 5 );
		gridBorderSize.Fill( 3 );    // Border for spline order = 3 ( 1
lower, 2 upper )
		totalGridSize = gridSizeOnImage + gridBorderSize;
		std::cout	<< "totalGridSize: " << totalGridSize << std::endl;

		bsplineRegion.SetSize( totalGridSize );

		typedef TransformType::SpacingType SpacingType;
		SpacingType spacing =
			fixedImage->GetSpacing();//spacing[0]=1.5625;spacing[1]=1.5625;
		std::cout << "spacing: "<< spacing<<std::endl;
		//spacing.Fill(1);
		//  fixedImage->SetSpacing(spacing);
		//  movingImage->SetSpacing(spacing);

		typedef TransformType::OriginType OriginType;
		OriginType origin =
			fixedImage->GetOrigin();//origin[0]=-9.0569601058959961;origin[1]=-191.00799560546875;
		std::cout << "origin: "<< origin <<std::endl;
		//origin.Fill(0);
		//	fixedImage->SetOrigin(origin);
		//	movingImage->SetOrigin(origin);

		FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();

		for(unsigned int r=0; r<ImageDimension; r++)
		{
			spacing[r] *= floor( static_cast<double>(fixedImageSize[r] - 1)  / 
				static_cast<double>(gridSizeOnImage[r] - 1) );
			origin[r]  -=  spacing[r]; 
		}

		transform->SetGridSpacing( spacing );
		transform->SetGridOrigin( origin );
		transform->SetGridRegion( bsplineRegion );


		typedef TransformType::ParametersType     ParametersType;

		const unsigned int numberOfParameters =
			transform->GetNumberOfParameters();

		std::cout << "numberOfParameters:" << numberOfParameters << std::endl;

		ParametersType parameters( numberOfParameters );

		parameters.Fill( 0.0 );

		transform->SetParameters( parameters );
		//  Software Guide : EndCodeSnippet



		//  Software Guide : BeginLatex
		//  
		//  We now pass the parameters of the current transform as the initial
		//  parameters to be used when the registration process starts. 
		//
		//  Software Guide : EndLatex 

		// Software Guide : BeginCodeSnippet
		registration->SetInitialTransformParameters( transform->GetParameters() );
		// Software Guide : EndCodeSnippet

		std::cout << "Intial Parameters = " << std::endl;
		std::cout << transform->GetParameters() << std::endl;

		//  Software Guide : BeginLatex
		//  
		//  Next we set the parameters of the LBFGS Optimizer. 
		//
		//  Software Guide : EndLatex 

		// Software Guide : BeginCodeSnippet
		optimizer->SetGradientConvergenceTolerance( 0.05 );
		optimizer->SetLineSearchAccuracy( 0.9 );
		optimizer->SetDefaultStepLength( 1.5 );
		optimizer->TraceOn();
		// optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
		//I think it is too big to waste time
		optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
		// Software Guide : EndCodeSnippet



		// Add a time probe
		itk::TimeProbesCollectorBase collector;

		std::cout << std::endl << "Starting Registration" << std::endl;

		try 
		{ 
			collector.Start( "Registration" );
			registration->StartRegistration(); 
			collector.Stop( "Registration" );
		} 
		catch( itk::ExceptionObject & err ) 
		{ 
			std::cerr << "ExceptionObject caught !" << std::endl; 
			std::cerr << err << std::endl; 
			return -1;
		} 

		OptimizerType::ParametersType finalParameters = 
			registration->GetLastTransformParameters();

		std::cout << "Last Transform Parameters" << std::endl;
		std::cout << finalParameters << std::endl;


		// Report the time taken by the registration
		collector.Report();

		//  Software Guide : BeginLatex
		//  
		//  Let's execute this example using the rat lung images from the
previous examples.
			//  
			// \begin{itemize}
			// \item \code{RatLungSlice1.mha}
			// \item \code{RatLungSlice2.mha}
			// \end{itemize}
			//
			//
			//  Software Guide : EndLatex 

			// Software Guide : BeginCodeSnippet
			transform->SetParameters( finalParameters );
		// Software Guide : EndCodeSnippet


		typedef itk::ResampleImageFilter< 
			MovingImageType, 
			FixedImageType >    ResampleFilterType;

		ResampleFilterType::Pointer resample = ResampleFilterType::New();

		resample->SetTransform( transform );
		resample->SetInput( movingImageReader->GetOutput() );

		resample->SetSize(    fixedImage->GetLargestPossibleRegion().GetSize() );


		// typedef TransformType::SpacingType SpacingType;
		SpacingType spacing1 = fixedImage->GetSpacing(); //not used
		// spacing1.Fill(1);

		// typedef TransformType::OriginType OriginType;
		OriginType origin1 = fixedImage->GetOrigin(); //not used
		// origin1.Fill(0);
		//   resample->SetOutputOrigin(  origin1 );
		//  resample->SetOutputSpacing( spacing1 );

		resample->SetOutputOrigin(  fixedImage->GetOrigin() );
		resample->SetOutputSpacing( fixedImage->GetSpacing() );
		resample->SetDefaultPixelValue( 100 );

		typedef  unsigned  short  OutputPixelType;

		typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;

		typedef itk::CastImageFilter< 
			FixedImageType,
			OutputImageType > CastFilterType;

		typedef itk::ImageFileWriter< OutputImageType >  WriterType;


		WriterType::Pointer      writer =  WriterType::New();
		CastFilterType::Pointer  caster =  CastFilterType::New();


		writer->SetFileName( "output.hdr" );


		caster->SetInput( resample->GetOutput() );
		writer->SetInput( caster->GetOutput()   );


		try
		{
			writer->Update();
		}
		catch( itk::ExceptionObject & err ) 
		{ 
			std::cerr << "ExceptionObject caught !" << std::endl; 
			std::cerr << err << std::endl; 
			return -1;
		} 



		typedef itk::SquaredDifferenceImageFilter< 
			FixedImageType, 
			FixedImageType, 
			OutputImageType > DifferenceFilterType;

		DifferenceFilterType::Pointer difference = DifferenceFilterType::New();

		WriterType::Pointer writer2 = WriterType::New();
		writer2->SetInput( difference->GetOutput() );  


		// Compute the difference image between the 
		// fixed and resampled moving image.
		if( 1 )
		{
			difference->SetInput1( fixedImageReader->GetOutput() );
			difference->SetInput2( resample->GetOutput() );
			writer2->SetFileName( "Diff_fixed_resampled.png" );
			try
			{
				writer2->Update();
			}
			catch( itk::ExceptionObject & err ) 
			{ 
				std::cerr << "ExceptionObject caught !" << std::endl; 
				std::cerr << err << std::endl; 
				return -1;
			} 
		}


		// Compute the difference image between the 
		// fixed and moving image before registration.
		if( 1)
		{
			writer2->SetFileName( "Diff_fixed_moving.png" );
			difference->SetInput1( fixedImageReader->GetOutput() );
			difference->SetInput2( movingImageReader->GetOutput() );
			try
			{
				writer2->Update();
			}
			catch( itk::ExceptionObject & err ) 
			{ 
				std::cerr << "ExceptionObject caught !" << std::endl; 
				std::cerr << err << std::endl; 
				return -1;
			} 
		}



		// Generate the explicit deformation field resulting from 
		// the registration.

		typedef itk::Vector< float, ImageDimension >  VectorType;
		typedef itk::Image< VectorType, ImageDimension >  DeformationFieldType;

		DeformationFieldType::Pointer field = DeformationFieldType::New();
		field->SetRegions( fixedRegion );
		field->SetOrigin( fixedImage->GetOrigin() );
		field->SetSpacing( fixedImage->GetSpacing() );
		std::cout << "fixedRegion " << fixedRegion << std::endl;	
		std::cout << "fixedImage->GetOrigin(): " << fixedImage->GetOrigin()<<
			std::endl;
		std::cout << "fixedImage->GetSpacing(): " << fixedImage->GetSpacing()
			<< std::endl;
		//  origin1.Fill(0);
		//  spacing1.Fill(1);
		//  field->SetOrigin(origin);
		//  field->SetSpacing(spacing);
		field->Allocate();

		typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;
		FieldIterator fi( field, fixedRegion );

		fi.GoToBegin();

		TransformType::InputPointType  fixedPoint;
		TransformType::OutputPointType movingPoint;
		DeformationFieldType::IndexType index;

		VectorType displacement;
		int n = 0;
		while( ! fi.IsAtEnd() )
		{
			index = fi.GetIndex();
			field->TransformIndexToPhysicalPoint( index, fixedPoint );
			movingPoint = transform->TransformPoint( fixedPoint );
			displacement[0] = movingPoint[0] - fixedPoint[0];
			displacement[1] = movingPoint[1] - fixedPoint[1];
			fi.Set( displacement );
			++fi;
			n++;
		}
		std::cout << "Iterator Number: " << n << std::endl;



		typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
		FieldWriterType::Pointer fieldWriter = FieldWriterType::New();

		fieldWriter->SetInput( field );

		if( 1 )
		{
			fieldWriter->SetFileName( "DeformationField.mhd" );
			try
			{
				fieldWriter->Update();
			}
			catch( itk::ExceptionObject & excp )
			{
				std::cerr << "Exception thrown " << std::endl;
				std::cerr << excp << std::endl;
				return EXIT_FAILURE;
			}
		}

		return 0;
	}



Is  there something I'm doign wrong here ??

regards
--Sachin
On Sat, 2 Oct 2004 18:33:42 -0400, Sachin Jambawalikar
<sachinjam at gmail.com> wrote:
> Hi Luis ,
>  Even when I try running from command line   the exe
> deformabelregistration4 and deformableregistration6 (after compiling
> itk with examples and no modifications )  I get the same errors for my
> datasets . But If  I use  Brain web simulated datasets they work .
> is it because of the size of my data  512 x 512 x 60 (0.68x0.68x3)
> for 3D  where as the brain data had size  181x 217x 180 ( 1x1x1)  or
> is it because of the contrast enhancement between the images.
> 
> You can download  two of my volumes  in meta  image format (I tried
> using dicom Images 2D with just 2 2D slices  but same result ) from:
> 
> http://www.ic.sunysb.edu/Stu/sjambawa/b_images/volumes.tar.gz
> 
> thanks ,
> 
> Best  Regards
> 
> --Sachin
> 
> 
> 
> 
> On Fri, 01 Oct 2004 15:00:52 -0400, Luis Ibanez <luis.ibanez at kitware.com> wrote:
> >
> > Hi Sachin,
> >
> > Did you modified the code in any way before
> > applying it to your data ?
> >
> > If so, please post to the list the *exact* .cxx
> > file that you used.
> >
> > It seems that the problem is related to the optimizer
> > setting, and not to the number of parameters.
> >
> > In particular we suspect of the lines: 265-269.
> >
> >    optimizer->SetGradientConvergenceTolerance( 0.05 );
> >    optimizer->SetLineSearchAccuracy( 0.9 );
> >    optimizer->SetDefaultStepLength( 1.5 );
> >
> > Since the error message is comming from the optimizer
> > and it is related to the convergence criteria.
> >
> >    Thanks
> >
> >      Luis
> >
> > --------------------------
> >
> >
> > Sachin Jambawalikar wrote:
> >
> > > Hi all,
> > >
> > > I am using the example code DeformableRegistration4 to register   2D
> > > slices of dynamci contrast breast mri data  .
> > >
> > > I get the following exception:
> > >
> > >
> > > ImageRegion (0033BF04)
> > >   Dimension: 2
> > >   Index: [0, 0]
> > >   Size: [512, 512]
> > >
> > >
> > > ImageRegion (0033DEA4)
> > >   Dimension: 2
> > >   Index: [0, 0]
> > >   Size: [512, 512]
> > >
> > >
> > > totalGridSize: [8, 8]
> > > spacing: [0.683594, 0.683594]
> > > origin: [-192.401, -149.008]
> > > numberOfParameters:128
> > > Intial Parameters =
> > > [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
> > > , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> > > 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> > >  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
> > > , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
> > >
> > > Starting Registration
> > >  F = 18707.5, GNORM = 198.347
> > > *************************************************
> > >    I   NFN    FUNC        GNORM       STEPLENGTH
> > > IFLAG= -1. LINE SEARCH FAILED. SEE DOCUMENTATION OF ROUTINE MCSRCH ERROR RETURN
> > >  OF LINE SEARCH: INFO=3 POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE  INCORRECT OR
> > > INCORRECT TOLERANCESvnl_lbfgs: ** EEEK **
> > > ExceptionObject caught !
> > >
> > > itk::ExceptionObject (00F2FB80)
> > > Location: "Unknown"
> > > File: g:\itk_vtk\itk\insighttoolkit-1.6.0\code\common\itkBSplineDeformableTransf
> > > orm.txx
> > > Line: 232
> > > Description: itk::ERROR: BSplineDeformableTransform(003337D0): Mismatched betwee
> > > n parameters size and region size
> > >
> > >
> > >
> > > Am i doing something wrong
> > >
> > > Regards
> > > --Sachin
> > > _______________________________________________
> > > Insight-users mailing list
> > > Insight-users at itk.org
> > > http://www.itk.org/mailman/listinfo/insight-users
> > >
> > >
> >
> >
>


More information about the Insight-users mailing list