[Insight-users] RE: using ThinPlate Splines transform.

[Luis Ibanez]

Hi Edvard,

The ThinPlateSplines should work on any dimension.
They are currently tested for 2D, 3D and 4D.

If your points are being mapped to a single slice
this may mean that your original target landmarks
are collapsed in 2D. You may want to verify the
third coordinate of all the landmarks that are
being feed into the ThinPlateSplineKernelTransform.


You may find useful to look at the code of the
application:

   InsightApplications/ThinPlateSplines
   http://www . itk . org/HTML/ThinPlateSplines . htm

Which maps points in 3D. (and visualize them
using VTK + FLTK).

Another place to look for code is the test in

Insight/Testing/Code/Common/itkSplineKernelTransformTest.cxx
http://www . itk . org/cgi-bin/cvsweb . cgi/Insight/Testing/Code/Common/itkSplineKernelTransformTest . cxx?cvsroot=Insight


Please let us know if you have further questions.


   Thanks


      Luis


-----------------------------
Edvard Antonin Dvorak wrote:
> Hi, I'm new to itk and have trawled through the postings for help with 
> the thin plate spline.
> I have written a 3d spline, having based it heavily on a downloaded 
> implementation of a 2d spline which works.
> Problem is that it does not seem to work. It seems to ignore the 3rd 
> dimension and concatenate all of the points to a single slice.
> 
> below is the code i am working with.
> I would appreciate any help
> 
> Yours sincerely
> Eddie
> 
> 
>   // typedefs for point
>   typedef itk::Point< double,3 >  PointType;
> 
>   // typedefs for vector
>   typedef std::vector< PointType >                   PointArrayType;
>      
>   // typedefs for spline transform
>   typedef itk::ThinPlateSplineKernelTransform< double,3 >  
> ThinPlateSplineTransformType;
> 
>   typedef ThinPlateSplineTransformType::Pointer  
> ThinPlateSplineTransformPointer;
> 
>   // typedefs for point sets acted on by spline transform
>   typedef ThinPlateSplineTransformType::PointSetType      PointSetType;
>   typedef PointSetType::Pointer                           PointSetPointer;
> 
>   ThinPlateSplineTransformPointer m_ThinPlateSplineTransformITK;
>   PointSetPointer   m_SourceLandMarks;  
>   PointSetPointer   m_TargetLandMarks;  
> 
> 
> void writeOffsetFile(ofstream& output,
>                      const int numPoints,
>                      const double* const xOffsets,
>                      const double* const yOffsets,
>                          const double* const zOffsets) {
>   for (int i = 0; i < numPoints; i++) {
>     output << xOffsets[i] << " " << yOffsets[i]<< " " << zOffsets[i]  << 
> endl;
>   }
>   if (output.fail()) {
>     cerr << "Error writing file." << endl;
>     exit(1);
>   }
> }
> 
> int main(int argc, char** argv)
> {
>   // parse command line arguments
>   cerr << "reading arguments..." << endl;
> 
>   int numFieldColumns   = atoi (argv[1]);
>   int numFieldRows              = atoi (argv[2]);
>   int numFieldSlices    = atoi (argv[3]);
>   int numQuieryPosX             = atoi (argv[4]);
>   int numQuieryIncX             = atoi (argv[5]);
>   int numQuieryPosY             = atoi (argv[6]);
>   int numQuieryIncY             = atoi (argv[7]);
>   int numQuieryPosZ             = atoi (argv[8]);
>   int numQuieryIncZ             = atoi (argv[9]);
>   char* landmarkFilename= argv[10];
>   char* displacementFieldFilename = argv[11];
> 
>   // create field generator and add query points
>   //
>   ThinPlateSplineDisplacementGenerator displacementGenerator; 
>  
>   // read landmark positions from file
>   cerr << "reading landmark positions..." << endl;
>   double* xLandmarkPositionsSource;
>   double* yLandmarkPositionsSource;
>   double* zLandmarkPositionsSource;
>   double* xLandmarkPositionsTarget;
>   double* yLandmarkPositionsTarget;
>   double* zLandmarkPositionsTarget;
>  
> 
>   ofstream input;
>   input.open(landmarkFilename);
> 
>   int numLandmarkPairs;
>   input >> numLandmarkPairs;
>   if (numLandmarkPairs < 0)
>   {
>     cerr << "Error reading file." << endl;
>     exit(1);
>   }
> 
>   xPositionsSource = new double[numLandmarkPairs];
>   yPositionsSource = new double[numLandmarkPairs];
>   zPositionsSource = new double[numLandmarkPairs];
>   xPositionsTarget = new double[numLandmarkPairs];
>   yPositionsTarget = new double[numLandmarkPairs];
>   zPositionsTarget = new double[numLandmarkPairs];
>  
>   if (xPositionsSource == 0
>       || yPositionsSource == 0
>           || zPositionsSource == 0
>       || xPositionsTarget == 0
>       || yPositionsTarget == 0
>           || zPositionsTarget == 0) {
>     cerr << "Error allocating memory." << endl;
>     exit(2);
>   }
> 
>   int currentLandmark = 0;
>   while (input >> xPositionsSource[currentLandmark]) {
>     input >> yPositionsSource[currentLandmark];
>         input >> zPositionsSource[currentLandmark];
>     input >> xPositionsTarget[currentLandmark];
>     input >> yPositionsTarget[currentLandmark];
>         input >> zPositionsTarget[currentLandmark];
>     currentLandmark++;
>   }
> 
>   if (currentLandmark != numLandmarkPairs) {
>     cerr << "Invalid number of landmark pairs." << endl;
>     exit(2);
>   }
> 
> 
>   // set landmarks
>   cerr << "setting landmarks..." << endl;
>   displacementGenerator.AddLandmarks(numLandmarks,
>                                      xLandmarkPositionsSource,
>                                      yLandmarkPositionsSource,
>                                          zLandmarkPositionsSource,
>                                      xLandmarkPositionsTarget,
>                                      yLandmarkPositionsTarget,
>                                          zLandmarkPositionsTarget);
>   delete [] xLandmarkPositionsSource;
>   delete [] yLandmarkPositionsSource;
>   delete [] zLandmarkPositionsSource;
>   delete [] xLandmarkPositionsTarget;
>   delete [] yLandmarkPositionsTarget;
>   delete [] zLandmarkPositionsTarget;
>  
>   //
>   // create array of query points (points to find displacement for)
>   //
>   cerr << "creating query points..." << endl;
>   int numQueries = (numFieldRows/numQuieryIncY) * 
> (numFieldColumns/numQuieryIncX) * (numFieldSlices/numQuieryIncZ);
>   double *xQueryPositions=new double[numQueries];
>   double *yQueryPositions=new double[numQueries];
>   double *zQueryPositions=new double[numQueries];
>   double *xDisplacements=new double[numQueries];
>   double *yDisplacements=new double[numQueries];
>   double *zDisplacements=new double[numQueries];
> 
>   if (xQueryPositions == 0
>       || yQueryPositions == 0
>           || zQueryPositions == 0
>       || xDisplacements == 0
>           || yDisplacements == 0
>       || zDisplacements == 0) {
>     cerr << "Can't allocate memory." << endl;   
>     return 2;
>   }
> 
>   // query points form regularly spaced grid
>   int arrayIndex = 0;
>   for (int z = numQuieryPosZ; z < numFieldSlices ;z+=numQuieryIncZ){
>           for (int x = numQuieryPosX; x < numFieldColumns; 
> x+=numQuieryIncX) {
>                   for (int y = numQuieryPosY; y < numFieldRows; 
> y+=numQuieryIncY) {
>                           xQueryPositions[arrayIndex] = x;
>                           yQueryPositions[arrayIndex] = y;
>                           zQueryPositions[arrayIndex] = z;
>                           ++arrayIndex;
>                   }
>           }
>   }
> 
>   //
>   // get displacement vectors
>   //
>   cerr << "computing displacement vectors..." << endl;
>   displacementGenerator.GetDisplacements(numQueries,
>                                          xQueryPositions,
>                                          yQueryPositions,
>                                          zQueryPositions,
>                                          xDisplacements,
>                                          yDisplacements,
>                                          zDisplacements);
>   cerr << "finished computing ..."<<endl;
> 
>   //
>   // write offsets to file
>   //
>   cerr << "writing displacement vectors..." << endl;
> 
>   ofstream output;
>   output.open(displacementFieldFilename);
>    for (int i = 0; i < numPoints; i++) {
>     output << xOffsets[i] << " " << yOffsets[i]<< " " << zOffsets[i]  << 
> endl;
>   }
>   if (output.fail()) {
>     cerr << "Error writing file." << endl;
>     exit(1);
>   }
> 
> 
>   input.close();
>   output.close();
>   return 1;
> }
> 
> 
> ----------------------------
> Edvard Antonin Dvorak
> edvorzak at medphys . ucl . ac . uk
> 02076796424
>