[Insight-users] Transform point
Cristina Ciavarro
cristinaciavarro at hotmail.com
Tue May 6 05:44:59 EDT 2008
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
More information about the Insight-users
mailing list