[Insight-users] Transform point

[Cristina Ciavarro]

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