[Insight-users] 3D bspline deformable registration works for
example data but gives exception for breast mri data
Luis Ibanez
luis.ibanez at kitware.com
Tue Oct 12 09:11:48 EDT 2004
Hi Sachin,
We were able to reproduce the error that you reported.
Could you please log a bug entry in the bug tracker:
http://www.itk.org/Bug/
You can create your own account by just providing your
email address as user Id.
Please let us know if you encounter any problem while
logging this bug entry.
Thanks
Luis
---------------------------
Sachin Jambawalikar wrote:
> 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
>>>>
>>>>
>>>
>>>
> _______________________________________________
> 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