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

1 view (last 30 days)
Ana on 23 Oct 2014
Commented: Ana on 24 Oct 2014
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...
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);

Alex Taylor on 23 Oct 2014
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?
##### 2 CommentsShowHide 1 older comment
Ana on 24 Oct 2014
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...

### Categories

Find more on DICOM Format in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!