<reponame>vespa-mrs/vespa
# Python modules
import os

# 3rd party modules
import pydicom

# Our modules
import vespa.common.util.dir_list as dir_list

 

# MIN & MAX PIXEL are limited by what can be stuffed into a byte, because
# we turn the DICOM data into RGB data with R == G == B for each pixel.
_MAX_PIXEL = 255
_MIN_PIXEL = 0


class DicomFileInfo(object):
    """
    A convenience container for a small subset of attributes from a DICOM file.

    The advantages of this over a pydicom Dataset object are --
    1) This uses Pythonic names
    2) All attributes are guaranteed to exist (although they may be empty)
       so no need to test "if 'foo' in dataset".
    3) This offers flexible sorting. Lists of DicomFileInfo objects will sort on
       the attributes listed in the class-level attribute 'sort_attributes'.
       To sort on e.g. study_id and series_number:

           files = get_all_files(path)
           DicomFileInfo.sort_attributes = ('study_id', 'series_number')
           files.sort()

       One cannot sort on arbitrary dataset objects, only on instance attributes.

    """
    sort_attributes = ('patient_name','study_id','series_number','instance_number')

    def __init__(self, dataset, filename=""):
        """
        The dataset param must be a pydicom dataset. The filename is
        stored in this object for convenience but is otherwise unused.

        'slice_location' is calculated from ImagePositionPatient (0020 0037)
        and ImageOrientationPatient (0020 0037) if present according to the
        algorithm outlined at great length by <NAME> here:
          http://www.itk.org/pipermail/insight-users/2003-September/004762.html
          http://www.cmake.org/pipermail/insight-users/2003-September/004762.html
        If those links die you can beseech Google:
          http://www.google.com/search?q=%22slightly-rotated+coronal+acquisition%22

        """
        self.dataset            = dataset
        self.filename           = filename
        self.patient_name       = dataset.PatientName       # 0010 0010
        self.patient_id         = dataset.PatientID         # 0010 0020
        self.study_instance_uid = dataset.StudyInstanceUID  # 0020 000D
        self.study_id           = ''                        # 0020 0010
        self.study_description  = ''                        # 0008 1030
        self.image_id           = 0                         # 0054 0400
        self.instance_number    = 0                         # 0020 0013
        self.series_number      = 0                         # 0020 0011
        self.series_description = ''                        # 0008 103E
        self.series_instance_uid = ''                       # 0020 000E
        self.echo_time          = 0.001                     # 0018 0081
        self.acquisition_time   = 0                         # 0008 0032

        # Here I populate the not-always-present values
        if 'ImageType' in dataset:         self.image_type          = dataset.ImageType         # 0008 0008
        if 'StudyID' in dataset:           self.study_id            = dataset.StudyID
        if 'ImageID' in dataset:           self.image_id            = dataset.ImageID
        if 'SeriesNumber' in dataset:      self.series_number       = dataset.SeriesNumber
        if 'StudyDescription' in dataset:  self.study_description   = dataset.StudyDescription
        if 'SeriesInstanceUID' in dataset: self.series_instance_uid = dataset.SeriesInstanceUID
        if 'SeriesDescription' in dataset: self.series_description  = dataset.SeriesDescription
        if 'InstanceNumber' in dataset:    self.instance_number     = dataset.InstanceNumber
        if 'EchoTime' in dataset:          self.echo_time           = dataset.EchoTime
        if 'AcquisitionTime' in dataset:   self.acquisition_time    = dataset.AcquisitionTime

        self.slice_location = 0
        if ('ImagePositionPatient'     in dataset)  and \
            ('ImageOrientationPatient' in dataset)  and \
            dataset.ImagePositionPatient            and \
            (len(dataset.ImageOrientationPatient) >= 6):

            o = dataset.ImageOrientationPatient
            slice_normals = [ (o[1] * o[5]) - (o[2] * o[4]),
                              (o[2] * o[3]) - (o[0] * o[5]),
                              (o[0] * o[4]) - (o[1] * o[3]), ]

            self.slice_location = sum([a * b for a, b in zip(slice_normals, dataset.ImagePositionPatient)])

    # bjs June 2020, Py3 no longer supports __cmp__() method, have to do all comparisons below

    def __lt__(self, other):
        """ Compares based on the class attribute sort_attributes"""
        item1, item2 = self.get_cmp_lists(other)
        return item1 < item2

    def __gt__(self, other):
        """ Compares based on the class attribute sort_attributes"""
        item1, item2 = self.get_cmp_lists(other)
        return item1 > item2

    def __eq__(self, other):
        """ Compares based on the class attribute sort_attributes"""
        item1, item2 = self.get_cmp_lists(other)
        return item1 == item2

    def __ne__(self, other):
        """ Compares based on the class attribute sort_attributes"""
        item1, item2 = self.get_cmp_lists(other)
        return item1 != item2

    def get_cmp_lists(self, other):
        item1 = [getattr(self, item) for item in self.sort_attributes]
        item2 = [getattr(other, item) for item in self.sort_attributes]
        return item1, item2


    def get_image(self):
        """Returns a wx.Image containing the image data in the dataset.
        Only a small subset of the image types described in the DICOM
        specification are supported. Specifically --
        - the dataset's PixelData attribute must be populated
        - the image must be grayscale
        - the data must be 8, 16, or 32-bit
        - the transfer syntax UID must indicate raw mode (big and little
          endian are both OK)

        If this method encounters something it doesn't understand, it raises
        a NotImplementedError.
        """
        import wx

        # I create a shortcut reference for the dataset
        ds = self.dataset

        if 'PixelData' not in ds:
            return None
        else:
            # The Transfer Syntax UID is prefixed with a noise string
            # (1.2.840.10008.1.) that's followed by the parts that matter. In
            # this case, '2' and '2.1' ==> raw little endian and '2.2' ==>
            # raw big endian.
            # ref: DICOM spec, Part 6: Data Dictionary, Annex A, Table A-1
            RAW_UIDS = ('1.2.840.10008.1.2', '1.2.840.10008.1.2.1',
                        '1.2.840.10008.1.2.2')

            # Turn the pydicom string into a list of integers
            pixels = [ord(c) for c in ds.PixelData]

            if ds.TransferSyntaxUID not in RAW_UIDS:
                raise NotImplementedError('Unsupported: TransferSyntaxUID (0002,0010) == %d' % ds.TransferSyntaxUID)
            else:
                # Data is in raw mode
                little_endian = (ds.TransferSyntaxUID != '1.2.840.10008.1.2.2')

                if ds.BitsAllocated not in (8, 16, 32):
                    raise NotImplementedError('Unsupported: BitsAllocated (0028,0100) == %d' % ds.BitsAllocated)

                # dataset.BitsAllocated is 8, 16, or 32
                bytes_per_pixel = ds.BitsAllocated / 8

                if bytes_per_pixel == 1:
                    # No need to massage the data
                    pass
                elif bytes_per_pixel == 2:
                    # Combine the byte sets into ints
                    for i in range(0, len(pixels) // 2):
                        if little_endian:
                            pixels[i] = (pixels[(i * 2) + 1] << 8) |   \
                                         pixels[ i * 2     ]
                        else:
                            # big endian
                            pixels[i] = (pixels[(i * 2)    ] << 8) |   \
                                         pixels[(i * 2) + 1]

                    # The second half of the data is now trash
                    pixels = pixels[:len(pixels) // 2]
                elif bytes_per_pixel == 4:
                    for i in range(0, len(pixels) // 4):
                        if little_endian:
                            pixels[i] = (pixels[(i * 4) + 3] << 24) |  \
                                        (pixels[(i * 4) + 2] << 16) |  \
                                        (pixels[(i * 4) + 1] <<  8) |  \
                                         pixels[(i * 4)    ]
                        else:
                            # big endian
                            pixels[i] = (pixels[(i * 4)    ] << 24) |  \
                                        (pixels[(i * 4) + 1] << 16) |  \
                                        (pixels[(i * 4) + 2] <<  8) |  \
                                         pixels[(i * 4) + 3]

                    # All but the first quarter of the data is now trash
                    pixels = pixels[:len(pixels) // 4]

            # Now pixels is a list of ints and we can forget about bits per
            # pixel and endianness.

            if ('PhotometricInterpretation' in ds) and \
               (not ds.PhotometricInterpretation.startswith('MONOCHROME')):
                # We only handle monochrome images
                raise NotImplementedError('Unsupported: PhotometricInterpretation (0028,0004) == %d' % ds.PhotometricInterpretation)
            else:
                # If Photometric Interpretation (0028,0004) isn't present,
                # I *assume* MONOCRHOME1 or MONOCRHOME2.
                if ('WindowCenter' in ds) and ('WindowWidth' in ds):
                    # Normalize pixels based on the algorithm presented in the
                    # DICOM spec (see the function for details)
                    normalize_pixels(pixels, ds.WindowCenter, ds.WindowWidth)

                if ('PhotometricInterpretation' in ds) and \
                   (ds.PhotometricInterpretation == 'MONOCHROME1'):
                    # MONOCHROME1 ==> low values are bright and high values
                    # are dark. In other words, it is the opposite of what
                    # most images libraries expect.
                    pixels = [_MAX_PIXEL - n for n in pixels]

                # Turn the pixel values into a string in preparation for
                # passing them to the image library.
                # wx.Image wants RGB data which is why I create three copies
                # of each character.
                pixels = "".join([chr(int(round(n))) * 3 for n in pixels])

                image = wx.EmptyImage(ds.Rows, ds.Columns)
                image.SetData(pixels)

            return image



def normalize_pixels(pixels, window_center, window_width):
    """
    Applies a linear conversion to the pixels based on the Window
    Center (0028,1050) and Window Width (0028,1051) values of the dataset.
    Pixels must be an iterable (e.g. list); it is altered in-place.
    The pixels returned are floating point values between _MIN_PIXEL and
    _MAX_PIXEL.

    This implements the algorithm presented in the DICOM specification
    Part 3: Information Object Definitions, Annex C.11.2.1.2.
    ftp://medical.nema.org/medical/dicom/2008/08_03pu.pdf

    Here's the pseudocode from the DICOM spec:
        if      (x <= c - 0.5 - (w-1)/2), then y = ymin
        else if (x >  c - 0.5 + (w-1)/2), then y = ymax
        else y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax - ymin) + ymin

    The Python code below is the same as the pseudocode with some loop
    invariants calculated here to improve performance.

    """
    min_threshold = (window_center - 0.5 - ((window_width - 1) / 2))
    max_threshold = (window_center - 0.5 + ((window_width - 1) / 2))
    output_range = _MAX_PIXEL - _MIN_PIXEL

    window_width  -= 1.0
    window_center -= 0.5

    for i, value in enumerate(pixels):
        if value <= min_threshold:
            pixels[i] = _MIN_PIXEL
        elif value > max_threshold:
            pixels[i] = _MAX_PIXEL
        else:
            pixels[i] = ((((value - window_center) / window_width) + 0.5) * output_range) + _MIN_PIXEL


def is_dicom(filename):
    """
    Returns True if the file in question is a DICOM file, else False. 
    
    - a DICOM file starts with 128 reserved bytes followed by 'DICM'.
    - ref: DICOM spec, Part 10: Media Storage and File Format for Media 
           Interchange, 7.1 DICOM FILE META INFORMATION 
    """
    if os.path.isfile(filename):
        s = open(filename, 'rb').read(132)
        if isinstance(s,(bytes, bytearray)):
            try:
                s = s.decode('utf-8')
            except:
                # bugfix - was trying to read PNG and found a char utf8 did not like
                try:
                    s = s.decode('utf-16')
                except:
                    return False
        return s.endswith('DICM')
    else:
        return False


def get_all_files(path):
    """
    Gets all DICOM files in the directory indicated by the path and
    returns them as an unsorted list of DicomFileInfo instances.
    To sort the files, call the sort() method of the returned list.
    """
    return [f for f in get_files(path)]


def get_files(path, manufacturer=['SIEMENS',]):
    """
    Return DicomFileInfo instances (or a subclass thereof) for the DICOM files
    in the directory indicated by the path. This is a Python generator so they
    are returned one by one.
    
    Example, use it in a loop like so:

        for dicom_file in get_files(the_path):
            do_something(dicom_file)

    See also get_all_files().

    """
    filenames = dir_list.DirList(path, dir_list.FilterFile).fq_names

    for filename in filenames:
        if is_dicom(filename):
            dataset = pydicom.dcmread(filename)
            dataset.decode()                        # change strings to unicode
            df = DicomFileInfo(dataset, filename)
            yield df
        #else:
            # Not a DICOM file; ignore it.


# Python modules
import os

# 3rd party modules
import pydicom

# Our modules
import vespa.common.util.dir_list as dir_list

 

# MIN & MAX PIXEL are limited by what can be stuffed into a byte, because
# we turn the DICOM data into RGB data with R == G == B for each pixel.
_MAX_PIXEL = 255
_MIN_PIXEL = 0


class DicomFileInfo(object):
    """
    A convenience container for a small subset of attributes from a DICOM file.

    The advantages of this over a pydicom Dataset object are --
    1) This uses Pythonic names
    2) All attributes are guaranteed to exist (although they may be empty)
       so no need to test "if 'foo' in dataset".
    3) This offers flexible sorting. Lists of DicomFileInfo objects will sort on
       the attributes listed in the class-level attribute 'sort_attributes'.
       To sort on e.g. study_id and series_number:

           files = get_all_files(path)
           DicomFileInfo.sort_attributes = ('study_id', 'series_number')
           files.sort()

       One cannot sort on arbitrary dataset objects, only on instance attributes.

    """
    sort_attributes = ('patient_name','study_id','series_number','instance_number')

    def __init__(self, dataset, filename=""):
        """
        The dataset param must be a pydicom dataset. The filename is
        stored in this object for convenience but is otherwise unused.

        'slice_location' is calculated from ImagePositionPatient (0020 0037)
        and ImageOrientationPatient (0020 0037) if present according to the
        algorithm outlined at great length by <NAME> here:
          http://www.itk.org/pipermail/insight-users/2003-September/004762.html
          http://www.cmake.org/pipermail/insight-users/2003-September/004762.html
        If those links die you can beseech Google:
          http://www.google.com/search?q=%22slightly-rotated+coronal+acquisition%22

        """
        self.dataset            = dataset
        self.filename           = filename
        self.patient_name       = dataset.PatientName       # 0010 0010
        self.patient_id         = dataset.PatientID         # 0010 0020
        self.study_instance_uid = dataset.StudyInstanceUID  # 0020 000D
        self.study_id           = ''                        # 0020 0010
        self.study_description  = ''                        # 0008 1030
        self.image_id           = 0                         # 0054 0400
        self.instance_number    = 0                         # 0020 0013
        self.series_number      = 0                         # 0020 0011
        self.series_description = ''                        # 0008 103E
        self.series_instance_uid = ''                       # 0020 000E
        self.echo_time          = 0.001                     # 0018 0081
        self.acquisition_time   = 0                         # 0008 0032

        # Here I populate the not-always-present values
        if 'ImageType' in dataset:         self.image_type          = dataset.ImageType         # 0008 0008
        if 'StudyID' in dataset:           self.study_id            = dataset.StudyID
        if 'ImageID' in dataset:           self.image_id            = dataset.ImageID
        if 'SeriesNumber' in dataset:      self.series_number       = dataset.SeriesNumber
        if 'StudyDescription' in dataset:  self.study_description   = dataset.StudyDescription
        if 'SeriesInstanceUID' in dataset: self.series_instance_uid = dataset.SeriesInstanceUID
        if 'SeriesDescription' in dataset: self.series_description  = dataset.SeriesDescription
        if 'InstanceNumber' in dataset:    self.instance_number     = dataset.InstanceNumber
        if 'EchoTime' in dataset:          self.echo_time           = dataset.EchoTime
        if 'AcquisitionTime' in dataset:   self.acquisition_time    = dataset.AcquisitionTime

        self.slice_location = 0
        if ('ImagePositionPatient'     in dataset)  and \
            ('ImageOrientationPatient' in dataset)  and \
            dataset.ImagePositionPatient            and \
            (len(dataset.ImageOrientationPatient) >= 6):

            o = dataset.ImageOrientationPatient
            slice_normals = [ (o[1] * o[5]) - (o[2] * o[4]),
                              (o[2] * o[3]) - (o[0] * o[5]),
                              (o[0] * o[4]) - (o[1] * o[3]), ]

            self.slice_location = sum([a * b for a, b in zip(slice_normals, dataset.ImagePositionPatient)])

    # bjs June 2020, Py3 no longer supports __cmp__() method, have to do all comparisons below

    def __lt__(self, other):
        """ Compares based on the class attribute sort_attributes"""
        item1, item2 = self.get_cmp_lists(other)
        return item1 < item2

    def __gt__(self, other):
        """ Compares based on the class attribute sort_attributes"""
        item1, item2 = self.get_cmp_lists(other)
        return item1 > item2

    def __eq__(self, other):
        """ Compares based on the class attribute sort_attributes"""
        item1, item2 = self.get_cmp_lists(other)
        return item1 == item2

    def __ne__(self, other):
        """ Compares based on the class attribute sort_attributes"""
        item1, item2 = self.get_cmp_lists(other)
        return item1 != item2

    def get_cmp_lists(self, other):
        item1 = [getattr(self, item) for item in self.sort_attributes]
        item2 = [getattr(other, item) for item in self.sort_attributes]
        return item1, item2


    def get_image(self):
        """Returns a wx.Image containing the image data in the dataset.
        Only a small subset of the image types described in the DICOM
        specification are supported. Specifically --
        - the dataset's PixelData attribute must be populated
        - the image must be grayscale
        - the data must be 8, 16, or 32-bit
        - the transfer syntax UID must indicate raw mode (big and little
          endian are both OK)

        If this method encounters something it doesn't understand, it raises
        a NotImplementedError.
        """
        import wx

        # I create a shortcut reference for the dataset
        ds = self.dataset

        if 'PixelData' not in ds:
            return None
        else:
            # The Transfer Syntax UID is prefixed with a noise string
            # (1.2.840.10008.1.) that's followed by the parts that matter. In
            # this case, '2' and '2.1' ==> raw little endian and '2.2' ==>
            # raw big endian.
            # ref: DICOM spec, Part 6: Data Dictionary, Annex A, Table A-1
            RAW_UIDS = ('1.2.840.10008.1.2', '1.2.840.10008.1.2.1',
                        '1.2.840.10008.1.2.2')

            # Turn the pydicom string into a list of integers
            pixels = [ord(c) for c in ds.PixelData]

            if ds.TransferSyntaxUID not in RAW_UIDS:
                raise NotImplementedError('Unsupported: TransferSyntaxUID (0002,0010) == %d' % ds.TransferSyntaxUID)
            else:
                # Data is in raw mode
                little_endian = (ds.TransferSyntaxUID != '1.2.840.10008.1.2.2')

                if ds.BitsAllocated not in (8, 16, 32):
                    raise NotImplementedError('Unsupported: BitsAllocated (0028,0100) == %d' % ds.BitsAllocated)

                # dataset.BitsAllocated is 8, 16, or 32
                bytes_per_pixel = ds.BitsAllocated / 8

                if bytes_per_pixel == 1:
                    # No need to massage the data
                    pass
                elif bytes_per_pixel == 2:
                    # Combine the byte sets into ints
                    for i in range(0, len(pixels) // 2):
                        if little_endian:
                            pixels[i] = (pixels[(i * 2) + 1] << 8) |   \
                                         pixels[ i * 2     ]
                        else:
                            # big endian
                            pixels[i] = (pixels[(i * 2)    ] << 8) |   \
                                         pixels[(i * 2) + 1]

                    # The second half of the data is now trash
                    pixels = pixels[:len(pixels) // 2]
                elif bytes_per_pixel == 4:
                    for i in range(0, len(pixels) // 4):
                        if little_endian:
                            pixels[i] = (pixels[(i * 4) + 3] << 24) |  \
                                        (pixels[(i * 4) + 2] << 16) |  \
                                        (pixels[(i * 4) + 1] <<  8) |  \
                                         pixels[(i * 4)    ]
                        else:
                            # big endian
                            pixels[i] = (pixels[(i * 4)    ] << 24) |  \
                                        (pixels[(i * 4) + 1] << 16) |  \
                                        (pixels[(i * 4) + 2] <<  8) |  \
                                         pixels[(i * 4) + 3]

                    # All but the first quarter of the data is now trash
                    pixels = pixels[:len(pixels) // 4]

            # Now pixels is a list of ints and we can forget about bits per
            # pixel and endianness.

            if ('PhotometricInterpretation' in ds) and \
               (not ds.PhotometricInterpretation.startswith('MONOCHROME')):
                # We only handle monochrome images
                raise NotImplementedError('Unsupported: PhotometricInterpretation (0028,0004) == %d' % ds.PhotometricInterpretation)
            else:
                # If Photometric Interpretation (0028,0004) isn't present,
                # I *assume* MONOCRHOME1 or MONOCRHOME2.
                if ('WindowCenter' in ds) and ('WindowWidth' in ds):
                    # Normalize pixels based on the algorithm presented in the
                    # DICOM spec (see the function for details)
                    normalize_pixels(pixels, ds.WindowCenter, ds.WindowWidth)

                if ('PhotometricInterpretation' in ds) and \
                   (ds.PhotometricInterpretation == 'MONOCHROME1'):
                    # MONOCHROME1 ==> low values are bright and high values
                    # are dark. In other words, it is the opposite of what
                    # most images libraries expect.
                    pixels = [_MAX_PIXEL - n for n in pixels]

                # Turn the pixel values into a string in preparation for
                # passing them to the image library.
                # wx.Image wants RGB data which is why I create three copies
                # of each character.
                pixels = "".join([chr(int(round(n))) * 3 for n in pixels])

                image = wx.EmptyImage(ds.Rows, ds.Columns)
                image.SetData(pixels)

            return image



def normalize_pixels(pixels, window_center, window_width):
    """
    Applies a linear conversion to the pixels based on the Window
    Center (0028,1050) and Window Width (0028,1051) values of the dataset.
    Pixels must be an iterable (e.g. list); it is altered in-place.
    The pixels returned are floating point values between _MIN_PIXEL and
    _MAX_PIXEL.

    This implements the algorithm presented in the DICOM specification
    Part 3: Information Object Definitions, Annex C.11.2.1.2.
    ftp://medical.nema.org/medical/dicom/2008/08_03pu.pdf

    Here's the pseudocode from the DICOM spec:
        if      (x <= c - 0.5 - (w-1)/2), then y = ymin
        else if (x >  c - 0.5 + (w-1)/2), then y = ymax
        else y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax - ymin) + ymin

    The Python code below is the same as the pseudocode with some loop
    invariants calculated here to improve performance.

    """
    min_threshold = (window_center - 0.5 - ((window_width - 1) / 2))
    max_threshold = (window_center - 0.5 + ((window_width - 1) / 2))
    output_range = _MAX_PIXEL - _MIN_PIXEL

    window_width  -= 1.0
    window_center -= 0.5

    for i, value in enumerate(pixels):
        if value <= min_threshold:
            pixels[i] = _MIN_PIXEL
        elif value > max_threshold:
            pixels[i] = _MAX_PIXEL
        else:
            pixels[i] = ((((value - window_center) / window_width) + 0.5) * output_range) + _MIN_PIXEL


def is_dicom(filename):
    """
    Returns True if the file in question is a DICOM file, else False. 
    
    - a DICOM file starts with 128 reserved bytes followed by 'DICM'.
    - ref: DICOM spec, Part 10: Media Storage and File Format for Media 
           Interchange, 7.1 DICOM FILE META INFORMATION 
    """
    if os.path.isfile(filename):
        s = open(filename, 'rb').read(132)
        if isinstance(s,(bytes, bytearray)):
            try:
                s = s.decode('utf-8')
            except:
                # bugfix - was trying to read PNG and found a char utf8 did not like
                try:
                    s = s.decode('utf-16')
                except:
                    return False
        return s.endswith('DICM')
    else:
        return False


def get_all_files(path):
    """
    Gets all DICOM files in the directory indicated by the path and
    returns them as an unsorted list of DicomFileInfo instances.
    To sort the files, call the sort() method of the returned list.
    """
    return [f for f in get_files(path)]


def get_files(path, manufacturer=['SIEMENS',]):
    """
    Return DicomFileInfo instances (or a subclass thereof) for the DICOM files
    in the directory indicated by the path. This is a Python generator so they
    are returned one by one.
    
    Example, use it in a loop like so:

        for dicom_file in get_files(the_path):
            do_something(dicom_file)

    See also get_all_files().

    """
    filenames = dir_list.DirList(path, dir_list.FilterFile).fq_names

    for filename in filenames:
        if is_dicom(filename):
            dataset = pydicom.dcmread(filename)
            dataset.decode()                        # change strings to unicode
            df = DicomFileInfo(dataset, filename)
            yield df
        #else:
            # Not a DICOM file; ignore it.
# Python modules # 3rd party modules # Our modules # MIN & MAX PIXEL are limited by what can be stuffed into a byte, because # we turn the DICOM data into RGB data with R == G == B for each pixel. A convenience container for a small subset of attributes from a DICOM file.      The advantages of this over a pydicom Dataset object are --     1) This uses Pythonic names     2) All attributes are guaranteed to exist (although they may be empty)        so no need to test "if 'foo' in dataset".     3) This offers flexible sorting. Lists of DicomFileInfo objects will sort on        the attributes listed in the class-level attribute 'sort_attributes'.        To sort on e.g. study_id and series_number:             files = get_all_files(path)            DicomFileInfo.sort_attributes = ('study_id', 'series_number')            files.sort()         One cannot sort on arbitrary dataset objects, only on instance attributes. The dataset param must be a pydicom dataset. The filename is         stored in this object for convenience but is otherwise unused.          'slice_location' is calculated from ImagePositionPatient (0020 0037)         and ImageOrientationPatient (0020 0037) if present according to the         algorithm outlined at great length by <NAME> here:           http://www.itk.org/pipermail/insight-users/2003-September/004762.html           http://www.cmake.org/pipermail/insight-users/2003-September/004762.html         If those links die you can beseech Google:           http://www.google.com/search?q=%22slightly-rotated+coronal+acquisition%22 # 0010 0010 # 0010 0020 # 0020 000D # 0020 0010 # 0008 1030 # 0054 0400 # 0020 0013 # 0020 0011 # 0008 103E # 0020 000E # 0018 0081 # 0008 0032 # Here I populate the not-always-present values # 0008 0008 # bjs June 2020, Py3 no longer supports __cmp__() method, have to do all comparisons below Compares based on the class attribute sort_attributes Compares based on the class attribute sort_attributes Compares based on the class attribute sort_attributes Compares based on the class attribute sort_attributes Returns a wx.Image containing the image data in the dataset.         Only a small subset of the image types described in the DICOM         specification are supported. Specifically --         - the dataset's PixelData attribute must be populated         - the image must be grayscale         - the data must be 8, 16, or 32-bit         - the transfer syntax UID must indicate raw mode (big and little           endian are both OK)          If this method encounters something it doesn't understand, it raises         a NotImplementedError. # I create a shortcut reference for the dataset # The Transfer Syntax UID is prefixed with a noise string # (1.2.840.10008.1.) that's followed by the parts that matter. In # this case, '2' and '2.1' ==> raw little endian and '2.2' ==> # raw big endian. # ref: DICOM spec, Part 6: Data Dictionary, Annex A, Table A-1 # Turn the pydicom string into a list of integers # Data is in raw mode # dataset.BitsAllocated is 8, 16, or 32 # No need to massage the data # Combine the byte sets into ints # big endian # The second half of the data is now trash # big endian # All but the first quarter of the data is now trash # Now pixels is a list of ints and we can forget about bits per # pixel and endianness. # We only handle monochrome images # If Photometric Interpretation (0028,0004) isn't present, # I *assume* MONOCRHOME1 or MONOCRHOME2. # Normalize pixels based on the algorithm presented in the # DICOM spec (see the function for details) # MONOCHROME1 ==> low values are bright and high values # are dark. In other words, it is the opposite of what # most images libraries expect. # Turn the pixel values into a string in preparation for # passing them to the image library. # wx.Image wants RGB data which is why I create three copies # of each character. Applies a linear conversion to the pixels based on the Window     Center (0028,1050) and Window Width (0028,1051) values of the dataset.     Pixels must be an iterable (e.g. list); it is altered in-place.     The pixels returned are floating point values between _MIN_PIXEL and     _MAX_PIXEL.      This implements the algorithm presented in the DICOM specification     Part 3: Information Object Definitions, Annex C.11.2.1.2.     ftp://medical.nema.org/medical/dicom/2008/08_03pu.pdf      Here's the pseudocode from the DICOM spec:         if      (x <= c - 0.5 - (w-1)/2), then y = ymin         else if (x >  c - 0.5 + (w-1)/2), then y = ymax         else y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax - ymin) + ymin      The Python code below is the same as the pseudocode with some loop     invariants calculated here to improve performance. Returns True if the file in question is a DICOM file, else False.           - a DICOM file starts with 128 reserved bytes followed by 'DICM'.     - ref: DICOM spec, Part 10: Media Storage and File Format for Media             Interchange, 7.1 DICOM FILE META INFORMATION # bugfix - was trying to read PNG and found a char utf8 did not like Gets all DICOM files in the directory indicated by the path and     returns them as an unsorted list of DicomFileInfo instances.     To sort the files, call the sort() method of the returned list. Return DicomFileInfo instances (or a subclass thereof) for the DICOM files     in the directory indicated by the path. This is a Python generator so they     are returned one by one.          Example, use it in a loop like so:          for dicom_file in get_files(the_path):             do_something(dicom_file)      See also get_all_files(). # change strings to unicode #else: # Not a DICOM file; ignore it.