reverted sitk back to pydicom. Fixed failing test.

pyalfe/data_structure.py

@@ -10,7 +10,7 @@
 from pathlib import Path
 
 from bids import BIDSLayout
-import SimpleITK as sitk
+import pydicom
 
 from strenum import StrEnum
 
@@ -46,7 +46,7 @@ class Tissue(IntEnum):
     CEREBELLUM = 6
 
 
-Series = collections.namedtuple('Series', ['series_uid', 'series_path', 'instances'])
+Series = collections.namedtuple('Series', ['series_path', 'instances'])
 
 
 class PipelineDataDir(ABC):
@@ -577,10 +577,41 @@ class PatientDicomDataDir:
     def __init__(self, dicom_dir) -> None:
         self.dicom_dir = dicom_dir
 
+    def is_dicom_file(self, file):
+        """Function to check if file is a dicom file.
+        The first 128 bytes are preamble the next 4 bytes should
+        contain DICM otherwise it is not a dicom file. Adapted from
+        https://github.com/icometrix/dicom2nifti/blob/main/dicom2nifti/common.py
+
+        Parameters
+        ----------
+        file: str
+            file to check for the DICM header block
+
+        Return
+        ------
+        bool
+            True if it is a dicom file otherwise False.
+        """
+        file_stream = open(file, 'rb')
+        file_stream.seek(128)
+        data = file_stream.read(4)
+        file_stream.close()
+        if data == b'DICM':
+            return True
+        try:
+            dicom_headers = pydicom.read_file(
+                file, defer_size="1 KB", stop_before_pixels=True, force=True
+            )
+            if dicom_headers is not None:
+                return True
+        except RuntimeError:
+            pass
+        return False
+
     def get_all_dicom_series_instances(self, accession):
-        """This method returns a dictionary where the keys are the series IDs
-        available for a patient-accession (study) and the values are the list
-        of all the instances in the series.
+        """This method returns a list of Series each containing a series path
+        and all of the instances in the series.
 
         Parameters
         ----------
@@ -589,21 +620,20 @@ def get_all_dicom_series_instances(self, accession):
 
         Returns
         -------
-        dict[str, list[str]]
-            A dictionary with series IDs as keys and
-            list of instances as values.
+        list[Series]
+            A list of Series (series_path, instances)
         """
-        ret = {}
+        ret = []
         accession_path = os.path.join(self.dicom_dir, accession)
+
         for path in glob.glob(os.path.join(accession_path, '*')):
-            series_uids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(path)
-            if len(series_uids) == 0:
-                continue
-            series_uid = series_uids[0]
-            instances = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(path, series_uid)
-            ret[series_uid] = Series(
-                series_uid=series_uid, series_path=path, instances=instances
-            )
+            instances = []
+            for file in glob.glob(os.path.join(path, '*')):
+                if os.path.isfile(file) and self.is_dicom_file(file):
+                    instances.append(file)
+            if len(instances) > 0:
+                ret.append(Series(series_path=path, instances=instances))
+
         return ret
 
     def get_series(self, accession, series_uid):

pyalfe/pipeline.py

@@ -78,12 +78,11 @@ def __init__(
             t1_registration,
             tissue_segmentation,
             resampling,
-            quantification
+            quantification,
         ]
         super().__init__(steps)
 
 
-
 class DicomProcessingPipelineRunner(PipelineRunner):
     """The Dicom processing pipeline runner
 
@@ -92,6 +91,7 @@ class DicomProcessingPipelineRunner(PipelineRunner):
     dicom_processing: DicomProcessing
         The dicom processing object.
     """
+
     def __init__(self, dicom_processing: DicomProcessing):
         steps = [dicom_processing]
-        super().__init__(steps)
+        super().__init__(steps)

pyalfe/tasks/dicom_processing.py

@@ -12,10 +12,21 @@
     Orientation,
 )
 from pyalfe.tasks import Task
-from pyalfe.utils.dicom import extract_dicom_meta, get_max_echo_series_crc, ImageMeta
+from pyalfe.utils.dicom import (
+    extract_image_meta,
+    get_max_echo_series_crc,
+    ImageMeta,
+)
 from pyalfe.utils.technique import detect_modality, detect_orientation
 
 
+def none_to_inf(val):
+    if val is not None:
+        return val
+    else:
+        return math.inf
+
+
 class DicomProcessing(Task):
     """This task classifies raw dicom images into their modality and converts
     them to NIfTI.
@@ -91,12 +102,6 @@ def get_best(image_list: list[ImageMeta]) -> ImageMeta:
             The image meta for the image with the smallest slice thickness.
         """
 
-        def none_to_inf(val):
-            if val is not None:
-                return val
-            else:
-                return math.inf
-
         return min(
             image_list, key=lambda image_meta: none_to_inf(image_meta.slice_thickness)
         )
@@ -105,48 +110,66 @@ def none_to_inf(val):
     def select_orientation(
         orientation_dict: dict[Orientation, ImageMeta]
     ) -> tuple[Orientation, ImageMeta]:
-        selected_orientation = None
-        selected_image = None
+        """This function selects the orientation from multiple orientations
+        choices. If Axial is available it chooses Axial. Otherwise it chooses
+        the orientation with the smallest slice thickness.
+
+        Parameters
+        ----------
+        orientation_dict : dict[Orientation, ImageMeta]
+            A dictionaly mapping orientaitons to image meta data.
+
+        Returns
+        -------
+        tuple[Orientation, ImageMeta]
+            The selected orientation and the coresponding image meta data.
+        """
 
         if Orientation.AXIAL in orientation_dict:
             selected_orientation = Orientation.AXIAL
-            selected_image = orientation_dict[Orientation.AXIAL]
         else:
-            min_thickness = float('inf')
-            for orientation in [Orientation.SAGITTAL, Orientation.CORONAL]:
-                if orientation in orientation_dict:
-                    image_meta = orientation_dict[orientation]
-                    if min_thickness > image_meta['slice_thickness']:
-                        min_thickness = image_meta['slice_thickness']
-                        selected_orientation = orientation
-                        selected_image = image_meta
-        return selected_orientation, selected_image
+            selected_orientation = min(
+                orientation_dict,
+                key=lambda orientation: none_to_inf(
+                    orientation_dict[orientation].slice_thickness
+                ),
+            )
+        return selected_orientation, orientation_dict[selected_orientation]
 
     def run(self, accession: str) -> None:
+        """Runs the dicom processing task for a given accession (study number).
+
+        Parameters
+        ----------
+        accession : str
+            The accession (study number).
+        """
         classified = defaultdict(defaultdict(list).copy)
         converted, conversion_failed = defaultdict(str), defaultdict(str)
         skipped = []
-        all_series = self.dicom_dir.get_all_dicom_series_instances(accession)
+        series_uid_map = {}
 
-        for series in all_series.values():
+        for series in self.dicom_dir.get_all_dicom_series_instances(accession):
+            dicom_meta = extract_image_meta(series.instances[0])
 
-            dicom_meta = extract_dicom_meta(
-                series.instances[0], series_uid=series.series_uid
-            )
             modality = detect_modality(
                 dicom_meta.series_desc,
                 dicom_meta.seq,
                 dicom_meta.contrast_agent,
                 dicom_meta.te,
             )
-            orientation = detect_orientation(dicom_meta.patient_orientation_vector)
 
+            orientation = detect_orientation(dicom_meta.patient_orientation_vector)
             if modality is None or orientation is None:
                 skipped.append(dicom_meta)
                 continue
 
             classified[modality][orientation].append(dicom_meta)
 
+            # we need this map so we can find the series path and instances
+            # for conversion NIfTI
+            series_uid_map[dicom_meta.series_uid] = series
+
         for modality, orientations in classified.items():
             orientation_best = defaultdict(Modality)
 
@@ -157,11 +180,13 @@ def run(self, accession: str) -> None:
             nifti_path = self.pipeline_dir.get_input_image(
                 accession=accession, modality=modality
             )
+
             selected_series_uid = selected_image_meta.series_uid
-            dcm_series_dir = all_series[selected_series_uid].series_path
+            dcm_series_dir = series_uid_map[selected_series_uid].series_path
             max_echo_series_crc = get_max_echo_series_crc(
-                all_series[selected_series_uid].instances
+                series_uid_map[selected_series_uid].instances
             )
+
             conversion_status = self.dicom2nifti(
                 dcm_series_dir=dcm_series_dir,
                 nifti_path=nifti_path,
@@ -172,4 +197,5 @@ def run(self, accession: str) -> None:
                 self.logger.warning(f'failed to convert {dcm_series_dir} to nifti.')
                 conversion_failed[modality] = selected_image_meta
                 continue
+
             converted[modality] = selected_image_meta