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);
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)];
TMB = [];
F = 0;
B = tformarray(A,T,R,TDA,TDB,TSB,TMB,F);
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);
