DICOM affine matrix transformation from image space to patient space in Matlab

Question

From the nifti header its easy to get the affine matrix. However in the DICOM header there are lots of entries, but its unclear to me which entries describe the transformation of which parameter to which new space.

I have found a tutorial which is quite detailed, but I cant find the entries they refer to. Also, that tutorial is written for Python, not Matlab. It lists these header entries:

Entries needed:
Image Position (0020,0032)
Image Orientation (0020,0037)
Pixel Spacing (0028,0030)

I cant find these if I load the header with dicominfo() . Maybe they are vendor specific or maybe they are nested away somewhere in the struct. Also the Pixel Spacing they refer to consists of two values, so I think their tutorial will only work for single slice transformations. More header entries about slice thickness and slicegap would be needed. Its also not easy to calculate the correct transformation for the z coordinates.

Does anybody know how to find these entries or how to transform image coordinates to patient coordinates with other information from a DICOM header? I use Matlab.

Answer 1

Ok so they were nested away in what might be a vendor specific entry of the struct. When loaded in Matlab, the name of the nest was inf.PerFrameFunctionalGroupsSequence.Item_X., then the framenumber, and then some more nesting which was more straightforward/self explanatory so I wont need to add it here. But search for the entries you need there. The slice spacing is called SpacingBetweenSlices (or slice thickness in the single slice case), the pixel spacing is called PixelSpacing and then there are ImagePositionPatient for the translation and ImageOrientationPatient for the rotation. Below is the code I wrote when following the steps from the nipy link below.

What happens is you load the direction cosines in a rotation matrix to align the basis vectors and you load the the pixel spacing and slice spacing in a matrix to scale the basis vectors and you load the image position to translate the new coordinate system. Finding the directoin cosines for the z direction takes some calculations because dicom apparently was designed for 2d images. In the single slice case the z direction cosines is the unit vector orthogonal to the x and y direction cosines (the cross product between the two) and in the multi slice case you can calculate it from all the differences in the translations between the slcies. After this you still want to apply the transformation which is also not immediately straightforward.

%load the header
inf = dicominfo(filename, 'dictionary', yourvendorspecificdictfilehere);

nSl = double(inf.MRSeriesNrOfSlices);
nY = double(inf.Height);
nX = double(inf.Width);
T1 = double(inf.PerFrameFunctionalGroupsSequence.Item_1.PlanePositionSequence.Item_1.ImagePositionPatient);

%load pixel spacing / scaling / resolution
RowColSpacing = double(inf.PerFrameFunctionalGroupsSequence.Item_1.PixelMeasuresSequence.Item_1.PixelSpacing);
%of inf.PerFrameFunctionalGroupsSequence.Item_1.PrivatePerFrameSq.Item_1.Pixel_Spacing;
dx = double(RowColSpacing(1));
dX = [1; 1; 1].*dx;%cols
dy = double(RowColSpacing(2));
dY = [1; 1; 1].*dy;%rows
dz = double(inf.SpacingBetweenSlices);%inf.PerFrameFunctionalGroupsSequence.Item_1.PrivatePerFrameSq.Item_1.SliceThickness; %thickness of spacing?
dZ = [1; 1; 1].*dz;

%directional cosines per basis vector
dircosXY = double(inf.PerFrameFunctionalGroupsSequence.Item_1.PlaneOrientationSequence.Item_1.ImageOrientationPatient);
dircosX = dircosXY(1:3);
dircosY = dircosXY(4:6);
if nSl == 1;
    dircosZ = cross(dircosX,dircosY);%orthogonal to other two direction cosines!
else
    N = nSl;%double(inf.NumberOfFrames);
    TN = double(-eval(['inf.PerFrameFunctionalGroupsSequence.Item_',sprintf('%d', N),'.PlanePositionSequence.Item_1.ImagePositionPatient']));
    dircosZ = ((T1-TN)./nSl)./dZ;
end

%all dircos together
dimensionmixing = [dircosX dircosY dircosZ];

%all spacing together
dimensionscaling = [dX dY dZ];

%mixing and spacing of dimensions together
R = dimensionmixing.*dimensionscaling;%maps from image basis to patientbasis

%offset and R together
A = [[R T1];[0 0 0 1]];

%you probably want to switch X and Y
%(depending on how you load your dicom into a matlab array)
Aold = A;
A(:,1) = Aold(:,2);
A(:,2) = Aold(:,1);

This results in this affine formula:

So basically I followed this tutorial. What was the biggest struggle was getting the Z direction and the translation correct. Also finding identifying and converting the correct entries was not straightforward for me. I do think my answer adds something to that tutorial though, because it was pretty hard to find the entries they refer to and now I wrote some Matlab code getting the affine matrix from a DICOM header. Before using the found affine matrix you also might need to find the Z coordinates for all of your frames, which might not be trivial if your dataset has more than four dimensions (dicomread puts all higher dimensions in one big fourth dimension)

-Edit- Corrected Z direction and translation of the transformation