These are some notes on the algorithms that SPM uses to convert from DICOM to nifti. There are other notes in Siemens mosaic format.
The relevant SPM files are spm_dicom_headers.m, spm_dicom_dict.mat and spm_dicom_convert.m. These notes refer the version in SPM8, as of around January 2010.
This is obviously a Matlab .mat file. It contains variables group and element, and values, where values is a struct array, one element per (group, element) pair, with fields name and vr (the last a cell array).
Reads the given DICOM files into a struct. It looks like this was written by John Ahsburner (JA). Relevant fixes are:
When opening the DICOM file, SPM (subfunction readdicomfile)
The read_dicom subfunction reads a tag, then has a loop during which the tag is processed (by setting values into the return structure). At the end of the loop, it reads the next tag. The loop breaks when the current tag is empty, or is the item delimitation tag (group=FFFE, element=E00D).
After it has broken out of the loop, if the last tag was (FFFE, E00D) (item delimitation tag), and the tag length was not 0, then SPM sets the file pointer back by 4 bytes from the current position. JA comments that he didn’t find that in the standard, but that it seemed to be needed for the Philips Integra.
Tag lengths as read in read_tag subfunction. If current format is explicit (as in ‘explicit little endian’):
There’s a check for not-even tag length. If not even:
tag length of 13 set to tag length 10.
Written by John Ashburner and Jesper Andersson.
SPM makes a special case of Siemens ‘spectroscopy images’. These are images that have ‘SOPClassUID’ == ‘1.3.12.2.1107.5.9.1’ and the private tag of (29, 1210); for these it pulls out the affine, and writes a volume of ones corresponding to the acquisition planes.
For images that are not spectroscopy:
Next SPM distinguishes between Siemens mosaic format and standard DICOM.
Mosaic images are those with the Siemens private tag:
(0029, 1009) [CSA Image Header Version] OB: '20100114'
and a readable CSA header (see Siemens mosaic format), and with non-empty fields from that header of ‘AcquisitionMatrixText’, ‘NumberOfImagesInMosaic’, and with non-zero ‘NumberOfImagesInMosaic’. The rest are standard DICOM.
For converting mosaic format, see Siemens mosaic format. The rest of this page refers to standard (slice by slice) DICOMs.
Take first header, put as start of first volume. For each subsequent header:
Then, for each currently identified volume:
We now have a list of volumes, where each volume is a list of headers that may match.
For each volume:
For each volume, recalculate z coordinate as above. Calculate the z gaps. Subtract the mean of the z gaps from all z gaps. If the average of the (gap-mean(gap)) is greater than 1e-4, then print a warning that there are missing DICOM files.
This step happens if there were volumes with slices having the same z coordinate in the Second pass step above. The resort is on the set of DICOM headers that were in the volume, for which there were slices with identical z coordinates. We’ll call the list of headers that the routine is still working on - work_list.
This means - writing DICOM volumes from standard (slice by slice) DICOM datasets rather than Siemens mosaic format.
We need the (4,4) affine \(A\) going from voxel (array) coordinates in the DICOM pixel data, to mm coordinates in the DICOM patient coordinate system.
This section tries to explain how SPM achieves this, but I don’t completely understand their method. See Getting a 3D affine from a DICOM slice or list of slices for what I believe to be a simpler explanation.
First define the constants, matrices and vectors as in DICOM affine Definitions.
\(N\) is the number of slices in the volume.
Then define the following matrices:
For a volume with more than one slice (header), then \(a=1; b=1, c=N, d=1\). \(e, f, g\) are the values from \(T^N\), and \(h == 1\).
For a volume with only one slice (header) \(a=0, b=0, c=1, d=0\) and \(e, f, g, h\) are \(n_1 \Delta{s}, n_2 \Delta{s}, n_3 \Delta{s}, 0\).
The full transform appears to be \(A_{spm} = R L^{-1}\).
Now, SPM, don’t forget, is working in terms of Matlab array indexing, which starts at (1,1,1) for a three dimensional array, whereas DICOM expects a (0,0,0) start (see DICOM affine formula). In this particular part of the SPM DICOM code, somewhat confusingly, the (0,0,0) to (1,1,1) indexing is dealt with in the \(A\) transform, rather than the analyze_to_dicom transformation used by SPM in other places. So, the transform \(A_{spm}\) goes from (1,1,1) based voxel indices to mm. To get the (0, 0, 0)-based transform we want, we need to pre-apply the transform to take 0-based voxel indices to 1-based voxel indices:
This formula with the definitions above result in the single and multi slice formulae in 3D affine formulae.
See derivations/spm_dicom_orient.py for the derivations and some explanations.
Just apply scaling and offset from ‘RescaleSlope’ and ‘RescaleIntercept’ for each slice and write volume.