[Insight-users] Transform point
Luis Ibanez
luis.ibanez at kitware.com
Tue May 6 09:24:18 EDT 2008
Hi Cristina,
Thanks for the detailed and clear description of the problem.
The main suspect that comes to mind is that the collection of
input points may not be falling inside the BSpline grid.
When that happens, a null-deformation is applied to the points,
and only the Bulk transform takes effect in changing their
coordinates.
Potential reasons for this problem may be:
1) Lack of setting the Direction of the BSpline Grid
based on the direction of your fixed image
See for example, line 290 of the CVS version of
Insight/Examples/Registration/DeformableRegistration12.cxx
290: transform->SetGridDirection( gridDirection );
2) Accidental errors on the settings of origin and
spacing of the BSpline Grid.
3) It may be that the collection of input points are
in a coordinate system different from the Fixed
image coordinate system.
How did you collected these points ?
Item (1) is the most likely suspect.
What you may want to do, is to make sure that you
use ITK 3.6 (or the CVS version) and you may want
to take a look at the BSplineDeformableTransformInitializer
paper in the Insight Journal:
http://insight-journal.org/midas/handle.php?handle=1926/1338
If you use this class for initializing your BSplineDeformable
transform, then you will rule out items (1) and (2) from the
list above.
Please give it a try, and let us know what you find,
Thanks
Luis
--------------------------
Cristina Ciavarro wrote:
> Dear itk users,
> I have developped a code in which I give as INPUT InputPoints, a file that contains deformable coefficients estimated by deformable registration, a file that contains VersorX, versorY, versorZ, translationX, translationY, translationZ computed by rigid registration and as OUTPUT OutputPoints, thas is the new coordinates of InputPoints.
>
> The steps are:
> 1.Registration of an image: I write in two different text files rototranslation values and deformable coefficients.
> 2.Application of the transformation calculated from precedent registration to an image. The code is correct because the output image is identical to image of point (1).
> 3.Application of the transformation calculated from registration to a PointSet.
>
> I have a problem in (3): outputPoints are only rototranslated. There isn't the applcation of deformable coefficients. Infact, if I comment the row of code transformSpline->SetBulkTransform( bulkTransform ); outputPoints are identical to inputPoints.
>
> What is the problem?
> Thank you!
> Best regards,
> Cristina
>
> #if defined(_MSC_VER)
> #pragma warning ( disable : 4786 )
> #endif
> #include
> #include
> #include
> #include "itkPointSet.h"
> #include "itkTransform.h"
> #include "itkDefaultStaticMeshTraits.h"
> #include "itkMultiResolutionImageRegistrationMethod.h"
> #include "itkMattesMutualInformationImageToImageMetric.h"
> #include "itkBSplineInterpolateImageFunction.h"
> #include "itkImage.h"
> #include "itkBSplineDeformableTransform.h"
> #include "itkVersorRigid3DTransform.h"
> #include "itkCenteredTransformInitializer.h"
>
> int main( int argc, char *argv[] )
> {
> if( argc < 6 )
> {
> std::cerr << "Missing Parameters " << std::endl;
> std::cerr << "Usage: " << argv[0];
> std::cerr << " fixedPointsFile outputPointsfile ";
> std::cerr << " DeformCoeff RigidValue InputImage TestImage";
>
> return EXIT_FAILURE;
> }
> const unsigned int Dimension=3;
>
> typedef itk::PointSet< double,3,itk::DefaultStaticMeshTraits>PointSetType;
>
> PointSetType::Pointer fixedPointSet = PointSetType::New();
> PointSetType::Pointer outputPointSet = PointSetType::New();
>
> typedef PointSetType::PointType PointType;
> typedef PointSetType::PointsContainer PointsContainer;
>
> PointsContainer::Pointer fixedPointContainer = PointsContainer::New();
> PointsContainer::Pointer outputPointContainer = PointsContainer::New();
>
> PointType fixedPoint;
> PointType outputPoint;
>
> // Read the file containing coordinates of fixed points.
> std::ifstream fixedFile;
> fixedFile.open( argv[1] );
> if( fixedFile.fail() )
> {
> std::cerr << "Error opening points file with name : " << std::endl;
> std::cerr << argv[1] << std::endl;
> return 2;
> }
> unsigned int pointId = 0;
> fixedFile>> fixedPoint;
> while( !fixedFile.eof() )
> {
> fixedPointContainer->InsertElement( pointId, fixedPoint );
> fixedFile>> fixedPoint;
> pointId++;
> }
> fixedPointSet->SetPoints( fixedPointContainer );
> std::cout <<
> "Number of fixed Points = " << fixedPointSet->GetNumberOfPoints()
> << std::endl;
>
> const unsigned int ImageDimension = 3;
> typedef signed short PixelType;
>
> typedef itk::Image< PixelType, ImageDimension> FixedImageType;
> typedef itk::Image< PixelType, ImageDimension> MovingImageMarkerType;
>
> const unsigned int SpaceDimension = ImageDimension;
> const unsigned int SplineOrder = 3;
> typedef double CoordinateRepType;
>
> typedef itk::BSplineDeformableTransform<
> CoordinateRepType,
> SpaceDimension,
> SplineOrder> TransformSplineType;
>
>
> typedef itk::VersorRigid3DTransform BulkTransformType;
>
> typedef itk::CenteredTransformInitializer< BulkTransformType,
> FixedImageType, MovingImageMarkerType> TransformInitializerType;
>
>
> TransformSplineType::Pointer transformSpline = TransformSplineType::New();
> BulkTransformType::Pointer bulkTransform = BulkTransformType::New();
> TransformInitializerType::Pointer initializer = TransformInitializerType::New();
>
> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
> typedef itk::ImageFileReader< MovingImageMarkerType> MovingImageReaderMarkerType;
>
> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
> MovingImageReaderMarkerType::Pointer movingImageMarkerReader = MovingImageReaderMarkerType::New();
>
> fixedImageReader->SetFileName( argv[5] );
> movingImageMarkerReader->SetFileName( argv[6] );
> fixedImageReader->Update();
> movingImageMarkerReader->Update();
> //fixedImageReader->GetOutput()->Print(std::cout);
> //movingImageMarkerReader->GetOutput()->Print(std::cout);
>
>
> initializer->SetTransform( bulkTransform );
> initializer->SetFixedImage( fixedImageReader->GetOutput() );
> initializer->SetMovingImage( movingImageMarkerReader->GetOutput() );
> initializer->GeometryOn();
> initializer->InitializeTransform();
>
> //Read input transform file
> std::cout << std::endl;
> std::cout << "Reading input transformfile: " << argv[4] << std::endl;
> std::cout << std::endl;
> std::ifstream finTrans( argv[4] );
>
> //Set bulk transform parameters
> typedef BulkTransformType::ParametersType BulkParametersType;
> BulkParametersType bulkParameters( bulkTransform->GetNumberOfParameters() );
> std::cout << "GetNumberOfParameters bulk "<<
> bulkTransform->GetNumberOfParameters()<< std::endl;
>
> finTrans>> bulkParameters[0]>> bulkParameters[1]>> bulkParameters[2];
> finTrans>> bulkParameters[3]>> bulkParameters[4]>> bulkParameters[5];
> std::cout << bulkParameters<GetOutput();
> FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
>
>
> typedef TransformSplineType::RegionType RegionType;
> RegionType bsplineRegion;
> RegionType::SizeType gridSizeOnImage;
> RegionType::SizeType gridBorderSize;
> RegionType::SizeType totalGridSize;
>
> gridSizeOnImage.Fill( 17 );
> gridSizeOnImage[2]=6;
> gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2 upper )
> totalGridSize = gridSizeOnImage + gridBorderSize;
>
> bsplineRegion.SetSize( totalGridSize );
>
> typedef TransformSplineType::SpacingType SpacingType;
> SpacingType spacing = fixedImage->GetSpacing();
>
> typedef TransformSplineType::OriginType OriginType;
> OriginType origin = fixedImage->GetOrigin();;
>
> FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
>
> for(unsigned int r=0; r(fixedImageSize[r] - 1) /
> static_cast(gridSizeOnImage[r] - 1) );
> origin[r] -= spacing[r];
> }
>
> transformSpline->SetGridSpacing( spacing );
> transformSpline->SetGridOrigin( origin );
> transformSpline->SetGridRegion( bsplineRegion );
>
> std::cout << "spacing[0]= "<< spacing[0] << std::endl;
> std::cout << "spacing[1]= "<< spacing[1] << std::endl;
> std::cout << "spacing[2]= "<< spacing[2] << std::endl;
>
> std::cout << "origin[0]= "<< origin[0] << std::endl;
> std::cout << "origin[1]= "<< origin[1] << std::endl;
> std::cout << "origin[2]= "<< origin[2] << std::endl;
>
>
> //Read input transform file
> std::cout << std::endl;
> std::cout << "Reading input transformfile: " << argv[3] << std::endl;
> std::cout << std::endl;
> std::ifstream trasformazioneElastica( argv[3] );
>
>
> //Set deformable transform parameters
> typedef TransformSplineType::ParametersType ParametersType;
> ParametersType finalParameters( transformSpline->GetNumberOfParameters() );
> std::cout << "GetNumberOfParameters "<< transformSpline->GetNumberOfParameters()<< std::endl;
>
> for (int i=0; i< 10800; i++)
> {
>
> trasformazioneElastica>> finalParameters[i];
> }
>
> typedef BulkTransformType::CenterType CenterType;
> CenterType center;
> center[0]=0;
> center[1]=0;
> center[2]=-235;
> bulkTransform->SetCenter( center );
> bulkTransform->SetParameters( bulkParameters );
>
> transformSpline->SetBulkTransform( bulkTransform );
>
> transformSpline->SetParameters( finalParameters );
>
> for(int i = 0; iGetNumberOfPoints(); i++)
> {
> fixedPointSet->GetPoint( i, &fixedPoint );
> outputPoint = transformSpline->TransformPoint( fixedPoint );
> outputPointContainer->InsertElement( i, outputPoint );
> }
> outputPointSet->SetPoints(outputPointContainer);
>
> //Write the new pointset into file
> std::ofstream outputFile;
> outputFile.open( argv[2]);
> if( outputFile.fail() )
> {
> return -1;
> }
> for(int i = 0; iGetNumberOfPoints(); i++)
> {
> outputPointSet->GetPoint( i, & outputPoint );
> for(int j = 0; j
> _________________________________________________________________
> Una cena tra amici? Cerca il ristorante con Live Search Maps
> http://maps.live.it
> _______________________________________________
> 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