[Insight-users] Image interpolation problem

[Marcus Alley]

Hi,

I'm new to itk, and as such I'm having trouble figuring out how to
solve my particular problem.  I have non-uniformly spaced 3D image
data from which I would like to extract interpolated data at various
points.

I've attempted to solve this using the itk kernel based transforms,
and I've included a simple example program at the end of this note.
The program works quite well for small image dimensions, but the
problem occurs if I set any one of the dimensions (nx, ny, nz) to
something more appropriate (for my problem) like 256.  In this case
the ComputeWMatrix() either crashes with a segmentation fault or else
takes far too long to be really useful for my application.

I would appreciate any advice on better ways to approach the problem.
I'm happy to clarify anything that I might have left out.

Thanks very much,
Marc

vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
Example:

#include "itkThinPlateSplineKernelTransform.h"

using namespace std;

int main()
{
  typedef itk::ThinPlateSplineKernelTransform<float, 4> TransformType;
  const int nx = 4;
  const int ny = 4;
  const int nz = 4;

  // Target landmark 
  TransformType::PointSetType::Pointer tLandmarks = 
    TransformType::PointSetType::New();
  tLandmarks->GetPoints()->Reserve(nx*ny*nz);

  // Source landmark set
  TransformType::PointSetType::Pointer sLandmarks = 
    TransformType::PointSetType::New();
  sLandmarks->GetPoints()->Reserve(nx*ny*nz);
  
  // Iterators
  TransformType::PointsIterator tIt = tLandmarks->GetPoints()->Begin();
  TransformType::PointsIterator sIt = sLandmarks->GetPoints()->Begin();

  for (int x = 0; x < nx; x++) {
    for (int y = 0; y < ny; y++) {
      for (int z = 0; z < nz; z++) {
	// Source points
	TransformType::InputPointType ptS;
	ptS[0] = x;
	ptS[1] = y; 
	ptS[2] = z; 
	ptS[3] = 0;
	sIt.Value() = ptS;
	sIt++;
	// Target points
	TransformType::InputPointType ptT;
	ptT[0] = 2*x;
	ptT[1] = 2*y; 
	ptT[2] = 2*z; 
	ptT[3] = 2*ptT[0] + 3*ptT[1];
	tIt.Value() = ptT;
	tIt++;
      }
    }
  }

  // Create the transform
  TransformType::Pointer mytrans = TransformType::New();
  mytrans->SetSourceLandmarks(sLandmarks);
  mytrans->SetTargetLandmarks(tLandmarks);
  mytrans->ComputeWMatrix();
  
  sIt = sLandmarks->GetPoints()->Begin();
  TransformType::PointsIterator endIt = sLandmarks->GetPoints()->End();
  while (sIt != endIt) {
    TransformType::InputPointType pt0 = sIt.Value();
    TransformType::InputPointType pt1 = mytrans->TransformPoint(pt0);
    cout << pt0 << " warps to: " << pt1 << endl;
    sIt++;
  }

  cout << endl;
  
  return 0;
}
  
****************************************************************
Marcus Alley, PhD		        | mtalley at stanford.edu
Research Associate			| office: 650.498.5581
Department of Radiology                 | fax: 650.723.5795
Radiological Sciences Laboratory        
Stanford University School of Medicine  
Lucas MRS Center, Room P-272
Stanford, California, 94305-5488
****************************************************************