Calculating voxel resolution from DICOM headers

Hi all,

What is the most accurate method to take common DICOM tags and calculate voxel resolution? For instance, how do FOV and percent sampling impact the calculation of pixel spacing?

Many thanks!
Maggie

You can determine in plane resolution from the DICOM tag 0028,0030 Pixel Spacing Attribute, which is specified by a numeric pair - adjacent row spacing (delimiter) adjacent column spacing in mm (note the order, in NIfTI the order of the three spatial dimensions are column [i], row [j] and slice [k]).

In theory, one can determine distance between voxel centers in the slice direction using DICOM tag 0018,0088 Spacing Between Slices - which should be the sum of the SliceThickness (0018,0050) and the slice gap (slice gaps are common in MR to avoid interference between slices). Unfortunately, not all vendors used 0018,0088 consistently and therefore the robust method to determine distance between voxel centers in the slice direction is to calculate the distance between the 0020,0032 ImagePositionPatient of two neighboring slices. While this value is virtually always the same for all slices in a MRI sequence, in CT it is common to have variable slice thickness (e.g. for the head, more slices in the critical-for-life brain stem than the cortex).

Here is an example for the critical tags from a typical MRI scan:

(0028,0030) DS [1.75\1.75] # 10, 2 PixelSpacing
(0018,0050) DS [2.5] # 4, 1 SliceThickness
(0018,0088) DS [2.5] # 4, 1 SpacingBetweenSlices
(0020,0032) DS [-114.58726620674-104.09609556198-74.025070190429] # 50, 3 ImagePositionPatient

Thank you for the information. Very helpful!

Follow up, how does percent sampling impact voxel size? Or are the values from the pixel spacing dicom tag equivalent to voxel size (non-slice dimensions)?

I have a bit of code based on an equation for calculating voxel resolution that I can no longer find and I’m looking to determine the accuracy of the approach.

        fov_frequency = round(dicom_header.PixelSpacing[0] * dicom_header.Columns)

        # FOV_phase = FOV * percent phase FOV
        # Percent phase FOV (also known as RFOV) is the ratio of FOV dimension in phase
        # direction to the FOV dimension in the frequency direction.
        rfov = dicom_header.PercentPhaseFieldOfView / 100
        fov_phase = round(dicom_header.PixelSpacing[1] * dicom_header.Rows * rfov)

        # Percent sampling is the fraction of acquisition matrix lines acquired.
        percent_sampling = dicom_header.PercentSampling / 100

        if (
            dicom_header.InPlanePhaseEncodingDirection == "ROW"
            or dicom_header.InPlanePhaseEncodingDirection == "OTHER"
        ):

            acquisition_matrix_frequency = dicom_header.AcquisitionMatrix[1]
            acquisition_matrix_phase = round(
                dicom_header.AcquisitionMatrix[2] * rfov * percent_sampling
            )

        elif dicom_header.InPlanePhaseEncodingDirection == "COL":

            acquisition_matrix_frequency = dicom_header.AcquisitionMatrix[0]
            acquisition_matrix_phase = round(
                dicom_header.AcquisitionMatrix[3] * rfov * percent_sampling
            )

        pixel_size_frequency = fov_frequency / acquisition_matrix_frequency
        pixel_size_phase = fov_phase / acquisition_matrix_phase

        voxel_x = float("{:03.3f}".format(pixel_size_frequency))
        voxel_y = float("{:03.3f}".format(pixel_size_phase))
        voxel_z = float("{:03.3f}".format(dicom_header.SliceThickness))

Do you have any suggestions? Or stick to the pixel spacing for reporting voxel resolution?

Many thanks for the advice!
Maggie

@Maggie

I would suggest you test your algorithm on some DICOM images. Here are the results for DICOMs where interpolation, partial FOV and mosaics are acquired.

I would suggest the method to use would be:

pixdim_1 = dicom_header.PixelSpacing[1]
pixdim_2 = dicom_header.PixelSpacing[0]
pixdim_3 = dicom_header.SliceThickness

as long as you can assume that the vendor is accurately reporting the distance between voxel centers in SliceThickness (compare distance between the ImagePositionPatient of two adjacent slices for a more robust solution).

Your example does not appear robust in my sample data (and a sample dataset with rectangular voxels will reveal that you need to swap the order of pixel spacing):

>./pixdim.py PF_PAT_POS_BW_INTERP_TEST.MR.HEAD_HARMS.0029.0001.2017.09.20.16.16.44.39766.120172126.IMA
@margaretmahan pixdim 19.2 19.2 2.4
@Chris_Rorden pixdim 1.2000000476837 1.2000000476837 2.4000000953674
>./pixdim.py PF_PAT_POS_BW_INTERP_TEST.MR.HEAD_HARMS.0033.0001.2017.09.20.16.16.44.39766.120190570.IMA
@margaretmahan pixdim 19.2 19.2 2.4
@Chris_Rorden pixdim 2.4000000953674 2.4000000953674 2.4000000953674
>./pixdim.py PF_PAT_POS_BW_INTERP_TEST.MR.HEAD_HARMS.0035.0001.2017.09.20.16.16.44.39766.120198766.IMA
@margaretmahan pixdim 19.2 19.2 2.4
@Chris_Rorden pixdim 2.4000000953674 2.4000000953674 2.4000000953674

Here is the code to replicate this:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# python3 d_speed_size.py        : test speed/compression for folder 'corpus'
# python3 d_speed_size.py indir  : test speed/compression for folder 'indir'
import os
import sys
import pydicom as dcm


if __name__ == '__main__':
    """Compare speed and size for different compression tools

    Parameters
    ----------
    fnm : str
        name of DICOM file to inspect

    """
    if len(sys.argv) > 1:
        fnm = sys.argv[1]
    else:
        print('No filename provided')
        sys.exit()
    if not os.path.isfile(fnm):
        print('Unable to find "' + fnm +'"')
        sys.exit()
            
    dicom_header = dcm.read_file(fnm)
    # FOV_phase = FOV * percent phase FOV
    # Percent phase FOV (also known as RFOV) is the ratio of FOV dimension in phase
    # direction to the FOV dimension in the frequency direction.
    rfov = dicom_header.PercentPhaseFieldOfView / 100
    fov_phase = round(dicom_header.PixelSpacing[1] * dicom_header.Rows * rfov)
    fov_frequency =  round(dicom_header.PixelSpacing[0] * dicom_header.Columns)

    # Percent sampling is the fraction of acquisition matrix lines acquired.
    percent_sampling = dicom_header.PercentSampling / 100

    if (
        dicom_header.InPlanePhaseEncodingDirection == "ROW"
            or dicom_header.InPlanePhaseEncodingDirection == "OTHER"
    ):
        acquisition_matrix_frequency = dicom_header.AcquisitionMatrix[1]
        acquisition_matrix_phase = round(dicom_header.AcquisitionMatrix[2] * rfov * percent_sampling)
    elif dicom_header.InPlanePhaseEncodingDirection == "COL":
        acquisition_matrix_frequency = dicom_header.AcquisitionMatrix[0]
        acquisition_matrix_phase = round(
            dicom_header.AcquisitionMatrix[3] * rfov * percent_sampling
    )

    pixel_size_frequency = fov_frequency / acquisition_matrix_frequency
    pixel_size_phase = fov_phase / acquisition_matrix_phase

    voxel_x = float("{:03.3f}".format(pixel_size_frequency))
    voxel_y = float("{:03.3f}".format(pixel_size_phase))
    voxel_z = float("{:03.3f}".format(dicom_header.SliceThickness))
    print('@margaretmahan pixdim {0} {1} {2}'.format(voxel_x, voxel_y, voxel_z))
    pixdim_1 = dicom_header.PixelSpacing[1]
    pixdim_2 = dicom_header.PixelSpacing[0]
    pixdim_3 = dicom_header.SliceThickness
    print('@Chris_Rorden pixdim {0} {1} {2}'.format(pixdim_1, pixdim_2, pixdim_3))