affinetransform - Interpolating or resampling algorithm in matlab for transformed 3D image, preferably sinc interpolation

Question

I have one 3D dataset and one 2D dataset which is a slice through the first volume. They are at different scales, resolutions and in a different coordinate system, but of both I know the affine transformation to world coordinates. I then think I know how to apply those, but how do I get an image back from those transformed coordinates again using sinc interpolation? I would like to learn how to perform this/how this works. The first comments below already pointed me at existing functions inside matlab that perform linear interpolation, but I would also like to know how to do this myself so I can use sinc interpolation (and others).

I can round the coordinates and get values for those, which would be nearest neighbor interpolation. I would like to lose as little information as possible irrespective of computation time, I think I should use sinc interpolation then. When I have the transformed coordinates, how do I make (for example sinc) an interpolation algorithm?

for example:

    %%get data
A = rand(200,250,250); % I want to get the slice from this that corresponds to B
B = rand(1200,1200); % I want to get the data from A that corresponds to the same pixels

    %%get coordinates for A
siza = size(A); clear xx;clear yy;clear zz;
[xx,yy,zz] = meshgrid(1:siza(1), 1:siza(2), 1:siza(3));
coora = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];
    %%get coordinates for B
sizb = size(B); clear xx;clear yy;clear zz;
[xx,yy] = meshgrid(1:sizb(1),1:sizb(2)); zz = zeros(size(xx));
coorb = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];

    %%define affine transformation matrices
T3d = [-0.02  0.02  1    -88 ;
       -0.98  0    -0.02  130;
        0     0.98 -0.02 -110;
        0     0     0     1  ];
T2d = [-0.2   0     0     126;
        0     0.2  -0.2  -131;
        0     0     2     43 ;
        0     0     0     1  ];

    %%transform A coordinates to world coordinates and world coordinates to B coordinates
cooraInBref = T3d*inv(T2d)*coora;
aslice = zeros(size(B));

    %% then nearest neighbor would go something like this (dont exactly know how to do this either):
cooraInBround = round(cooraInBref);
for idx = 1:length(coorb);
    if cooraInBround(3,idx) == 0
        aslice(cooraInBround(1,idx),cooraInBround(2,idx)) = ...;% dont know how to do this
    end
end
    %% how would I implement sinc interpolation instead of rounding the transformed coordinates

Related questions that dont quite help me further:

Matlab 3D data interpolation

Nearest-neighbor interpolation algorithm in MATLAB

"Nearest neighbor"-like interpolation in MATLAB

How to apply an affine transformation (4x4 matrix) to ndgrid/meshgrid results?

Interpolating 2D Matrix Data

scattered data interpolation

Python/PIL affine transformation

How do I rotate a 3D matrix by 90 degrees counterclockwise?

Resizing a 3D image (and resampling)

image transformations

Ready to use functions from matlab as pointed out in comments by chappjc anonsubmitter85 and cape code:

http://mathworks.com/help/matlab/math/interpolating-scattered-data.html

http://mathworks.com/help/matlab/ref/griddata.html?refresh=true

My other SE question I use to get the affine matrix

Using scatteredInterpolant as suggested in the comments and the linked question took almost ten minutes and resulted in a slice with all NaN`s, am now trying other methods.

Answer 1

I now start the other way around: get the 2d points to sample and then transform those to the 3d volume and then use interp3 to calculate the values. Its pretty quick and works well. This code only works for getting a single slice, but i think you can easily adapt it to get a whole transformed volume. I still dont know how to do sinc interpolation though.

%% transformation from one image to the other
Affine = inv(T2d)*T3d

%% get coordinates for B
sizb = size(B); clear xx;clear yy;clear zz;
[xx,yy] = meshgrid(1:sizb(1),1:sizb(2)); 
zz = ones(size(xx));
coorb = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];

%% transformed coordinates
coorb_t = Affine*coorb;
idxX = reshape(coorb_t(:,1), sizb(1), sizb(2), 1);
idxY = reshape(coorb_t(:,2), sizb(1), sizb(2), 1);
idxZ = reshape(coorb_t(:,3), sizb(1), sizb(2), 1);

%% interpolate
Asliced = interp3(A, idxX, idxY, idxZ, 'cubic');

Still not sure if I should have used zeros or ones for Z.