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

1 次查看(过去 30 天)
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);

回答(1 个)

Alex Taylor
Alex Taylor 2014-10-23
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 个评论
Ana
Ana 2014-10-23
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
Ana 2014-10-24
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...

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 DICOM Format 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by