Shifting problem in output results when using affine tformarray interpolation with DICOM images

I have two 3D DICOM images: PT (2x2x2mm voxel size) and CT (1x1x2mm voxel size), and I want to resample the first one (lower resolution) to the second one (higher resoltion) by using the function tformarray. I first calculate the matrices which characterize both images: T_CT and T_PT, which are 4x4 matrices where the first three rows correspond to the voxel size on each axis (img_XX.dspa(i) for i=1 --> x, i=2 --> y, and i=3 --> z) and the following three rows correspond to the image corner in DICOM coordinates (CENTER of the corner voxel, img_XX.dpcj(i) for each axis j=x,y,z). Then I used them to calculate the matrix of the transformation T from T_PT/T_CT in maketform.

T_CT = eye(4);
T_CT(1,1) =  img_CT.dspa(1);
T_CT(2,2) =  img_CT.dspa(2);
T_CT(3,3) = -img_CT.dspa(3);
T_CT(4,1) =  img_CT.dpcx(1);
T_CT(4,2) =  img_CT.dpcy(1);
T_CT(4,3) =  img_CT.dpcz(1);

T_PT = eye(4);
T_PT(1,1) =  img_PT.dspa(1);
T_PT(2,2) =  img_PT.dspa(2);
T_PT(3,3) = -img_PT.dspa(3);
T_PT(4,1) =  img_PT.dpcx(1);
T_PT(4,2) =  img_PT.dpcy(1);
T_PT(4,3) =  img_PT.dpcz(1);

A = PT;
T = maketform('affine',T_PT/T_CT); %transformation matrix
R = makeresampler('linear','fill'); 
TDA = [1,2,3];
TDB = [1,2,3];
TSB = [img_CT.ddim(1),img_CT.ddim(2),img_CT.ddim(3)]; % dimensions of CT image
TMB = [];
F = 0; % outside of target
B = tformarray(A,T,R,TDA,TDB,TSB,TMB,F); % resampled PT
clear PT; clear A;
PTonCT = B; clear B;

By doing this I get a perfect interpolation of the PT image to the higher resolution. However, it is slightly shifted in the three directions. At a first glance, the shift seems to be between 1 and 1/2 voxel. I tried by shifting the DICOM corner by the voxel lenght in each direction (SEE the two new matrices T_CT and T_PT below) and I got much better results, which makes me think about a problem of indexation of something like that. I mean, matlab considers starting at voxel (1,1,1) instead of (0,0,0) for example. However, I think the shift is still not totally compensated. Any idea of where this problem is coming from? Am I missunderstanding the input data and I must introduce the corner of the corner voxel instead of the center like in DICOM convention or either it is really a problem of first index? Where can I have more information about what this function does and the format of the input data? I already looked at mathworks but it is not very clear...

Thanks in advance

T_CT = eye(4);
T_CT(1,1) =  img_CT.dspa(1);
T_CT(2,2) =  img_CT.dspa(2);
T_CT(3,3) = -img_CT.dspa(3);
T_CT(4,1) = -img_CT.dspa(1) + img_CT.dpcx(1);
T_CT(4,2) = -img_CT.dspa(2) + img_CT.dpcy(1);
T_CT(4,3) = -img_CT.dspa(3) + img_CT.dpcz(1);

T_PT = eye(4);
T_PT(1,1) =  img_PT.dspa(1);
T_PT(2,2) =  img_PT.dspa(2);
T_PT(3,3) = -img_PT.dspa(3);
T_PT(4,1) = -img_PT.dspa(1) + img_PT.dpcx(1);
T_PT(4,2) = -img_PT.dspa(2) + img_PT.dpcy(1);
T_PT(4,3) = -img_PT.dspa(3) + img_PT.dpcz(1);

Answer 1

[Alex Taylor]

In the default coordinate system used for images in MATLAB, the center of the first pixel is (1,1,1), which implies that the UL corner is at (0.5,0.5,0.5). This is a common cause of similar bugs. I haven't worked through the logic of your code, but I would start there.

Any reason you aren't using imwarp for this?

[Ana]

Alex,thanks for your comment. This is exactly what I am doing, using the fist pixel (1,1,1) as the center of the corner pixel. So I dont understand then why I get good results only when I shift a distance equal to the pixel size in the 3 directions. About imwarp, I've never use it. Always used tformarry, but I could try.Do you think I would get rid of this problem with imwarp?

[Ana]

Hi again! Apparently imwarp is only for 2D and I using 3D images. Just for info, I tried to use also interp3 and I got again the results with a shift. However, this time the shift was even bigger...