多模态DICOM变换矩阵

Question

我在MATLAB中有这个函数，它获取与同一患者的两种不同磁共振成像模式相对应的两个体积的DICOM信息结构，并输出M，即从info1描述的体积的图像坐标转换为info2描述的体积的图像坐标的仿射变换矩阵：

function [M,Rot] = GetTransformMatrix(info1, info2)
% This function calculates the 4x4 transform and 3x3 rotation matrix 
% between two image coordinate system. 
% M=Tipp*R*S*T0;
% Tipp:translation
% R:rotation
% S:pixel spacing
% T0:translate to center(0,0,0) if necessary
% info1: dicominfo of 1st coordinate system
% info2: dicominfo of 2nd coordinate system
% Rot: rotation matrix between coordinate system

[Mdti,Rdti] = TransMatrix(info1);
[Mtf,Rtf] = TransMatrix(info2);
% First we transform into patient coordinates by multiplying by Mdti, and
% then we convert again into image coordinates of the second volume by
% multiplying by inv(Mtf)
M =  inv(Mtf) * Mdti;
Rot = inv(Rtf) * Rdti;
end

function [M,R] = TransMatrix(info)
%This function calculates the 4x4 transform matrix from the image
%coordinates to patient coordinates. 
ipp=info.ImagePositionPatient;
iop=info.ImageOrientationPatient;
ps=info.PixelSpacing;
Tipp=[1 0 0 ipp(1); ...
      0 1 0 ipp(2); ...
      0 0 1 ipp(3); ...
      0 0 0 1];
r=iop(1:3);  c=iop(4:6); s=cross(r',c');
R = [r(1) c(1) s(1) 0; ...
     r(2) c(2) s(2) 0; ...
     r(3) c(3) s(3) 0; ...
     0 0 0 1];
if info.MRAcquisitionType=='3D' % 3D turboflash
    S = [ps(2) 0 0 0; 0 ps(1) 0 0; 0 0 info.SliceThickness 0 ; 0 0 0 1];
else % 2D epi dti
    S = [ps(2) 0 0 0; ...
         0 ps(1) 0 0; ...
         0 0 info.SpacingBetweenSlices 0; ...
         0 0 0 1];
end
T0 = [ 1 0 0 0; ...
       0 1 0 0; ...
       0 0 1 0; ...
       0 0 0 1];
M = Tipp * R * S * T0;
end

为了将ADC图像坐标转换为T2图像坐标，我将以以下方式获得转换矩阵：

[M,Rot] = GetTransformMatrix(ADCHeader,T2Header);
M = M';

然后我获得了MATLAB仿射变换对象：

tform = affine3d(M);

要随后向T2卷注册ADCVolume，请执行以下操作：

[ADCVolume_reg,~] = imwarp(ADCVolume,tform,'Interp','cubic','FillValues',0,'OutputView',imref3d(size(T2Volume)));

当ADC和T2之间的DICOM参数info.SpacingBetweenSlices相同时，这会产生正确的结果，尽管当它不同时，我得到了错误的注册。

这可能是GetTransformMatrix函数的错误，或者是我传递给imwarp的参数错误？

Answer 1

毕竟，代码是正确的。唯一的问题是，我加载的头文件并不对应于所有卷的第一个片(因为我使用的加载器加载了带有第一个InstanceNumber (0020,0013)的文件头，这并不总是与第一个片相对应)，这使得我得到了错误的ImagePositionPatient值(0020,0032)。

使用此加载程序后，所有问题都得到了解决：

http://www.mathworks.com/matlabcentral/fileexchange/23455-dicom-directory--of-slices--to-3d-volume-image

其按照SliceLocation (0020,1041)对报头进行排序，该参数始终与片段编号相关。