[Insight-users] Image interpolation problem

[Luis Ibanez]

Hi Marcus,

By "non-uniformly" spaced 3D image do you mean an
"unstructured grid" ?

Depending on how your image data has been sample
you may find more efficient ways of doing this
interpolation with other ITK classes.

The Kernel splines, as you already noticed, are
not a good solution for large number of sampling
points. They will result in long computation times
and large memory consumption.

Please describe in detail the type of sampling
grid that you used for your "non-uniform" image.


   Thanks


      Luis


----------------------
Marcus Alley wrote:

> 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
> ****************************************************************
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>