2015-10-01 1 views
0

J'ai cette fonction dans MATLAB qui prend les structures d'information DICOM de deux volumes correspondant à deux modalités IRM différentes du même patient, et sort M, la matrice de transformation affine qui se transforme des coordonnées image du volume décrit par info1 dans l'image coordonnées du volume décrit par info2:Matrice de transformation DICOM multimodale

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 

Pour la conversion de l'image ADC coordonnées en coordonnées d'image T2, j'obtenir la matrice de transformation de la manière suivante:

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

Et puis-je obtenir l'objet de transformation affines Matlab:

tform = affine3d(M); 

Pour vous inscrire ensuite le ADCVolume au volume T2:

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

Ceci est ce qui donne le résultat correct lorsque les info.SpacingBetweenSlices des paramètres DICOM est la même chose entre ADC et T2, bien que je me trompe d'enregistrement quand il est différent.

Cela peut-il être une erreur avec la fonction GetTransformMatrix, ou est-ce que c'est sur les arguments que j'ai passés à imwarp?

Répondre

0

Après tout, le code est correct. Le seul problème était que les en-têtes que je chargeais ne correspondaient pas à la première tranche pour tous les volumes (parce que le chargeur que j'utilisais chargeait l'entête du fichier avec le premier InstanceNumber (0020, 0013), ce qui ne correspond à la première tranche), ce qui m'a fait obtenir la mauvaise valeur de ImagePositionPatient (0020,0032).

Tout se est résolu après avoir utilisé ce chargeur:

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

qui commande les en-têtes par SliceLocation (0020, 1041), un paramètre qui se rapporte toujours à trancher le numéro.