[Insight-users] Extraction position of points after TPS registration

[Agata Krasoń]

Hi All ;)

I  have a question.
I am trying to extract position of landmarks after TPS transformation.
I can't receive position of points after transformation.
I don't know what's wrong with this code ?
Could anyone look at this code ?
I would appreciate for any help please.
I need to receive points after TPS registration.

Here is my code :



#include "itkImageRegistrationMethod.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkOrientedImage.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkImage.h"
#include "itkVector.h"
#include "itkResampleImageFilter.h"
#include "itkLandmarkBasedTransformInitializer.h"
#include "itkRigid2DTransform.h"
#include "itkBSplineDeformableTransform.h"
#include "itkCenteredTransformInitializer.h"
#include "itkVersorRigid3DTransform.h"
#include "itkAffineTransform.h"
#include "itkBSplineDeformableTransform.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkBSplineResampleImageFunction.h"
#include "itkBSplineDecompositionImageFilter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
#include "itkTransformFileReader.h"
#include "itkCommand.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkThinPlateSplineKernelTransform.h"
#include "itkVolumeSplineKernelTransform.h"
#include "itkElasticBodySplineKernelTransform.h"
#include "itkPoint.h"
#include "itkPointSet.h"
#include "itkTransformFileWriter.h"



class CommandIterationUpdate : public itk::Command
{
public:
  typedef  CommandIterationUpdate   Self;
  typedef  itk::Command             Superclass;
  typedef itk::SmartPointer<Self>   Pointer;
  itkNewMacro( Self );
protected:
  CommandIterationUpdate() {};
public:
  typedef itk::RegularStepGradientDescentOptimizer  OptimizerType;
  typedef   const OptimizerType *                   OptimizerPointer;

  void Execute(itk::Object *caller, const itk::EventObject & event)
    {
    Execute( (const itk::Object *)caller, event);
    }

  void Execute(const itk::Object * object, const itk::EventObject & event)
    {
    OptimizerPointer optimizer =
      dynamic_cast< OptimizerPointer >( object );
    if( !(itk::IterationEvent().CheckEvent( &event )) )
      {
      return;
      }
    std::cout << optimizer->GetCurrentIteration() << "   ";
    std::cout << optimizer->GetValue() << "   ";
    std::cout << std::endl;
    }
};


using namespace std;


int main(int argc, char* argv[])
{

double tab[30] =
{
 // Tracker
-81.29,-31.07, -770.58,
    -83.11, -21.26, -822.64,
    -93.45, -32.44, -858.72,
    -68.08, -126.89,-813.07,
    -61.04, 75.74, -808.36,

 // Image
    140.6, 230.7, -30.5,
    140.2, 231.7, -71.1,
    144.8, 235.9, -116.1,
    45.8, 220.2, -66.7,
    231.6, 211.3, -66.1
 };

  const     unsigned int   Dimension = 3;
  typedef   unsigned char  PixelType;
  typedef   double CoordinateRepType;
  typedef   itk::ThinPlateSplineKernelTransform<
CoordinateRepType,Dimension>     TPSTransformType;
  typedef   itk::Point< CoordinateRepType,Dimension >  PointType;
  typedef   std::vector< PointType >                   PointArrayType;
  typedef   TPSTransformType::PointSetType      PointSetType;
  typedef   PointSetType::Pointer            PointSetPointer;
  typedef   PointSetType::PointIdentifier  PointIdType;


  // Landmarks correspondances may be associated with the
SplineKernelTransforms
  // via Point Set containers. Let us define containers for the landmarks.
  PointSetType::Pointer sourceLandMarks = PointSetType::New();
  PointSetType::Pointer targetLandMarks = PointSetType::New();
  PointType trackerPoint;     PointType imagePoint;
  PointSetType::PointsContainer::Pointer sourceLandMarkContainer =
                                   sourceLandMarks->GetPoints();
  PointSetType::PointsContainer::Pointer targetLandMarkContainer =
                                   targetLandMarks->GetPoints();




  // 1 Landmark
  trackerPoint[0] = tab[0];
  trackerPoint[1] = tab[1];
  trackerPoint[2] = tab[2];
  imagePoint[0] = tab[15];
  imagePoint[1] = tab[16];
  imagePoint[2] = tab[17];
  sourceLandMarkContainer->InsertElement( 0,trackerPoint);
  targetLandMarkContainer->InsertElement(0,imagePoint);



  // 2 Landmark
  trackerPoint[0] = tab[3];
  trackerPoint[1] = tab[4];
  trackerPoint[2] = tab[5];
  imagePoint[0] = tab[18];
  imagePoint[1] = tab[19];
  imagePoint[2] = tab[20];
  sourceLandMarkContainer->InsertElement(1,trackerPoint);
  targetLandMarkContainer->InsertElement(1,imagePoint);

  // 3 Landmark
  trackerPoint[0] = tab[6];
  trackerPoint[1] = tab[7];
  trackerPoint[2] = tab[8];
  imagePoint[0] = tab[21];
  imagePoint[1] = tab[22];
  imagePoint[2] = tab[23];
  sourceLandMarkContainer->InsertElement( 2,trackerPoint);
  targetLandMarkContainer->InsertElement(2,imagePoint);

  // 4 Landmark
  trackerPoint[0] = tab[9];
  trackerPoint[1] = tab[10];
  trackerPoint[2] = tab[11];
  imagePoint[0] = tab[24];
  imagePoint[1] = tab[25];
  imagePoint[2] = tab[26];
  sourceLandMarkContainer->InsertElement( 3,trackerPoint);
  targetLandMarkContainer->InsertElement(3,imagePoint);

  // 5 Landmark
  trackerPoint[0] = tab[12];
  trackerPoint[1] = tab[13];
  trackerPoint[2] = tab[14];
  imagePoint[0] = tab[27];
  imagePoint[1] = tab[28];
  imagePoint[2] = tab[29];

  sourceLandMarkContainer->InsertElement(4,trackerPoint);
  targetLandMarkContainer->InsertElement(4,imagePoint);

  TPSTransformType::Pointer tps = TPSTransformType::New();
  tps->SetSourceLandmarks(sourceLandMarks);
  tps->SetTargetLandmarks(targetLandMarks);
  tps->ComputeWMatrix();
  tps->UpdateParameters();




    PointSetType::Pointer transformed_points = PointSetType::New();
  PointSetType::PointsContainer::Pointer transformed_points_container =
transformed_points->GetPoints();
  TPSTransformType::OutputPointType trackerPointNewPosition;

  for(unsigned int i = 0; i < sourceLandMarks->GetNumberOfPoints(); i++){
*    PointType src_pnt = sourceLandMarks->GetPoint(i);  // Here It doesn't
compile  ?*
    trackerPointNewPosition = tps->TransformPoint(trackerPoint);
    transformed_points_container->InsertElement(i,trackerPointNewPosition);
  }



  return EXIT_SUCCESS;


}

agatte
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20120815/9da7eb34/attachment.htm>