dicom.py

from datetime import datetime
from enum import Enum
import random
from typing import Tuple

import numpy as np
import pydicom


"""
Alias for the main dicom datatype.
"""
DICOM = pydicom.dataset.FileDataset

"""
Since DICOM images only allow images stored in short integer format,
the Siemens scanner software multiplies values by a factor before storing
so that no precision is lost.
The scale can be found in this private DICOM tag.
"""
DCM_TAG_PIXEL_SCALE_FACTOR = 0x00331038
"""
Maximum value that can be stored in an unsigned integer with 16 bist.
"""
UINT16_MAX = 2**16 - 1
"""
DICOM images always contain UIDs to indicate their uniqueness.
Thus, when a DICOM image is updated, UIDs have to be changed for
which the following prefix is used.
"""
UID_PREFIX = "1.2.826.0.1.3680043.2.521."
"""
Default extension of dicom files.
"""
EXTENSION_DICOM = ".dcm"


class DICOMTime(Enum):
    """
    Class for parsing dates and times defined in a DICOM header
    into a python datetime type. It maps the four datetime fields
    [Study, Series, Acquisition, Content] of the DICOM header into
    an enumeration type.

    Usage: DICOMTime.Series.to_datetime(dcm)
    """

    Study = 1
    Series = 2
    Acquisition = 3
    Content = 4

    def date_field(self) -> str:
        """
        Get the name of the date field according to this DICOMTime type.
        """
        return f"{self.name}Date"

    def time_field(self) -> str:
        """
        Get the name of the time field according to this DICOMTime type.
        """
        return f"{self.name}Time"

    def to_datetime(self, dicom: DICOM) -> datetime:
        """
        Get the datetime according to this DICOMTime type.
        """
        _date = dicom[self.date_field()].value
        _time = dicom[self.time_field()].value
        return datetime(
            year=int(_date[0:4]),
            month=int(_date[4:6]),
            day=int(_date[6:8]),
            hour=int(_time[0:2]),
            minute=int(_time[2:4]),
            second=int(_time[4:6]),
        )


def parse_age(patient_age: str) -> int:
    """
    Parse an age string as defined in the DICOM standard into an integer representing the age in years.

    :param patient_age: age string as defined in the DICOM standard
    :return: the age in years as a number
    """
    assert (
        type(patient_age) == str
    ), f"patient age needs to be a string and not {type(patient_age)}"
    assert (
        len(patient_age) == 4
    ), f"patient age [{patient_age}] has to be four characters long"
    _num, _format = patient_age[:3], patient_age[3]
    assert (
        _format == "Y"
    ), f"currently, only patient ages in years [Y] is supported, not [{_format}]"
    return int(_num)


def load_dcm(filename: str, direction: float = 0.0) -> Tuple[DICOM, np.ndarray]:
    """
    Load a DICOM image, the data as a numpy array and apply normalization of the Siemens SPECT/CT
    Scanner.

    :param filename: filename of the DICOM image
    :param direction: other than 0 changes the stacking order of slices to the desired direction
    :return: the dicom header and the scaled image as a numpy array
    """
    dcm = pydicom.dcmread(filename)
    image = dcm.pixel_array

    if DCM_TAG_PIXEL_SCALE_FACTOR in dcm:
        image = image / dcm[DCM_TAG_PIXEL_SCALE_FACTOR].value

    dcm, image = to_direction((dcm, image), direction=direction)
    return dcm, image


def load_dcm_img(
    filename: str,
    direction: float = 0.0,
) -> np.ndarray:
    """
    Load a DICOM image as a numpy array and apply normalization of the Siemens SPECT/CT
    Scanner.

    :param filename: filename of the DICOM image
    :param direction: other than 0 changes the stacking order of slices to the desired direction
    :return: the image scaled and loaded into a numpy array
    """
    _, image = load_dcm(filename, direction=direction)
    return image


def to_direction(
    dcm: Tuple[DICOM, np.ndarray], direction: float = 1
) -> Tuple[DICOM, np.ndarray]:
    """
    Change the stacking direction of slices in a DICOM image.
    The direction is specified by the sign of the [Spacing Between Slices](https://dicom.innolitics.com/ciods/nm-image/nm-reconstruction/00180088) parameter.

    :param dcm: a tuple of a DICOM header and an image
    :param direction: the desired stacking direction as a float (>0 positive, <0 negative, =0 stay as it is)
    :return: dicom header (with values updated) and image according to the desired stacking order
    """
    dcm, image = dcm
    spacing_between_slices = float(dcm.SpacingBetweenSlices)
    if spacing_between_slices * direction < 0:
        dcm.DetectorInformationSequence[0].ImagePositionPatient[2] *= (
            image.shape[0] * spacing_between_slices
        )
        dcm.SpacingBetweenSlices = -spacing_between_slices

        image = image[::-1]
    return dcm, image


def scale_image(image: np.ndarray, initial_scale=10000000) -> Tuple[np.ndarray, float]:
    """
    For memory efficiency, the Siemens SPECT/CT does not store images as floating point
    numbers, but as unsigned integers with 16 bits. In order to somewhat keep precision,
    the floating points are scaled with a factor of 10^x where x is chosen so in a way
    that keeps the numbers in range of uint16. This function replicated this process.

    :param image: an image in floating points format
    :param initial_scale: the initial scale which is reduced until the maximum number
    is smaller than the maximum uint16 number
    :return: the image scaled and converted to uint16 as well as the used scaling factor
    """
    if image.min() < 0.0:
        raise ValueError("Cannot scale images with negative values!")

    scale = initial_scale
    while (scale * image.max()) > UINT16_MAX:
        scale = scale / 10
    image = (image * scale).astype(np.uint16)
    return image, scale


def update_dcm(dcm: DICOM, image: np.ndarray) -> DICOM:
    """
    Update the image data in a DICOM file. This function scales the image, converts
    it to unsigned integers with 16 bits and updates the pixel data in the DICOM file.
    Additionally, other related tags in the DICOM header, such as image dimensions and
    maximum pixel values, are updated accordingly.
    Note that this function modifies the given DICOM file. If you want to keep the old
    one, you should copy it first.

    :param dcm: the DICOM file to be udpated
    :param image: the image put into the DICOM file
    :return: the updated DICOM file

    (600, 25), (1005, 4)
    """
    image, scale = scale_image(image)

    dcm.NumberOfFrames = image.shape[0]
    dcm.NumberOfSlices = image.shape[0]
    dcm.SliceVector = list(range(1, image.shape[0] + 1))
    dcm.Columns = image.shape[1]
    dcm.Rows = image.shape[2]
    dcm.PixelData = image.tobytes()
    dcm.WindowWidth = image.max()
    dcm.WindowCenter = image.max() / 2
    dcm.LargestImagePixelValue = image.max()
    dcm[DCM_TAG_PIXEL_SCALE_FACTOR].value = scale
    return dcm


def change_uid(dcm: DICOM) -> DICOM:
    """
    Change the UIDs (SeriesInstance and SOPInstance) in a DICOM header so that
    it becomes its own unique file. Note that this method does not guarantee
    that the UIDs are fully unique. Since the creation of UIDs is time dependent,
    this function should not be used to rapidly change many UIDs.

    :param dcm: the DICOM file to be udpated
    :return: the DICOM file with updated UIDs
    """
    dcm.SeriesInstanceUID = UID_PREFIX + str(
        random.randint(10000000000000, 99999999999999)
    )
    dcm.SOPInstanceUID = UID_PREFIX + str(
        random.randint(10000000000000, 99999999999999)
    )
    return dcm


def is_dicom(filename: str) -> bool:
    """
    Check if a file is a DICOM file or not.

    :param filename: the file to be checked
    :return: if it is a DICOM file or not
    """
    try:
        pydicom.dcmread(filename, stop_before_pixels=True)
        return True
    except:
        return False


if __name__ == "__main__":
    import argparse

    parser = argparse.ArgumentParser(
        description="Dump the header of a DICOM image",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "file", type=str, help="the DICOM file of which the header is dumped"
    )
    args = parser.parse_args()

    dcm = pydicom.dcmread(args.file)
    print(dcm)