How imwarp transfer points in Matlab?

Question

I am using Matlab to transform an image to target image. I have geometric transformation(tform).

for example this is my 'tform':

    1.0235    0.0022   -0.0607         0
   -0.0276    1.0002    0.0089         0
   -0.0170   -0.0141    1.1685         0
   12.8777    5.0311  -70.0325    1.0000

In Matlab2013, it is possible to do it easily using imwarp:

%nii = 3D MR Image
I = nii.img; 
dii=nii.hdr.dime.pixdim(2:4);
Rfixed=imref3d(size(I),dii(2),dii(1),dii(3));    
new_img= imwarp(old_img, Rfixed, tform, 'OutputView', Rfixed);

the result is perfect using imwarp (red lung in the image)

I need to know how imwarp is working, Then I wrote my own function

function [new_img] = aff3d(old_img, tform, range_x, range_y, range_z)

   [U, V, W] = ndgrid(range_x, range_y, range_z);
   xyz = [reshape(U,[],1)';reshape(V,[],1)';reshape(W,[],1)'];
   xyz = [xyz; ones(1,size(xyz,2))];


   uvw = tform.T * xyz;
   % avoid homogeneous coordinate  
   uvw = uvw(1:3,:)';


   xi = reshape(uvw(:,1), length(range_x),length(range_y),length(range_z));
   yi = reshape(uvw(:,2), length(range_x),length(range_y),length(range_z));
   zi = reshape(uvw(:,3), length(range_x),length(range_y),length(range_z));

   old_img = single(old_img);
   new_img = interp3(old_img,yi,xi,zi,'linear');

   ii = find(isnan(new_img));
   if(~isempty(ii))
      new_img(ii) = 0;
   end
end

the result of my function (more info) is not match with imwarp output ( red lung is not locating in a correct place), anybody can help me?

Answer 1

As was suggested by Ander, try multiplying by the inverse transformation:

Tinv = tform.invert();
TinvMatrix = Tinv.T;

So your code would become:

function [new_img] = aff3d(old_img, tform, range_x, range_y, range_z)
   [U, V, W] = ndgrid(range_x, range_y, range_z);
   xyz = [reshape(U,[],1)';reshape(V,[],1)';reshape(W,[],1)'];
   xyz = [xyz; ones(1,size(xyz,2))];

   tformInv = invert(tform);
   uvw = tformInv.T * xyz;
   % avoid homogeneous coordinate  
   uvw = uvw(1:3,:)';
   xi = reshape(uvw(:,1), length(range_x),length(range_y),length(range_z));
   yi = reshape(uvw(:,2), length(range_x),length(range_y),length(range_z));
   zi = reshape(uvw(:,3), length(range_x),length(range_y),length(range_z));
   old_img = single(old_img);
   new_img = interp3(old_img,yi,xi,zi,'linear');
   ii = find(isnan(new_img));
   if(~isempty(ii))
      new_img(ii) = 0;
   end
end

In your code, you are interpolating within the old_img to try to find the new_img which has been warped. This implies that what you want to do is use the inverse mapping that maps from the output image space to the input image space. You appear to be interpolating your old image using the forward mapping of points, which is incorrect.

http://blogs.mathworks.com/steve/2006/04/28/spatial-transforms-forward-mapping/ http://blogs.mathworks.com/steve/2006/05/05/spatial-transformations-inverse-mapping/

I would use the above links to review forward vs. inverse mapping. IMWARP uses inverse mapping.

Part of the reason for confusion is that when people think of geometric transformations, they generally think in terms of the forward mapping of how points map from the old_image to the new_image. For this reason, the "T" property of the affine3d transformation is phrased in terms of the forward mapping.

When geometric transformations need to be implemented in software, it's much easier/better to implement things in terms of the inverse mapping. That is what imwarp does, and that's why you need to invert the transformation while you are trying to reproduce the imwarp behavior. If you read the blog post on inverse mapping, this algorithm is exactly what IMWARP is doing.

The only wrinkle you will need to work through is what IMWARP does in non-default coordinate systems (using non-default spatial referencing objects) in the case where the WorldLimits are not evenly divisible by the discrete grid of pixels. This behavior is arbitrary, there is no "right" behavior. The IMWARP behavior is to to honor the requested resolution (PixelExtentInWorld) and to slightly adjust the world limits in this case.