Merge pull request #35 from amithjkamath/main

Feature to read RT Dose files

pyradise/data/__init__.py

@@ -1,5 +1,5 @@
 from .annotator import Annotator
-from .image import Image, ImageProperties, IntensityImage, SegmentationImage
+from .image import Image, ImageProperties, IntensityImage, DoseImage, SegmentationImage
 from .modality import Modality
 from .organ import Organ, OrganAnnotatorCombination
 from .subject import Subject

pyradise/data/image.py

@@ -16,7 +16,7 @@
 
 TransformInfo = TypeVar("TransformInfo")
 
-__all__ = ["Image", "IntensityImage", "SegmentationImage", "ImageProperties"]
+__all__ = ["Image", "IntensityImage", "SegmentationImage", "DoseImage", "ImageProperties"]
 
 
 class ImageProperties:
@@ -744,6 +744,34 @@ def __str__(self) -> str:
         return f"SegmentationImage: {self.organ.get_name()} / {self.annotator.get_name()}"
 
 
+class DoseImage(IntensityImage):
+    """A dose image class to specialize for properties of RTDose volumes.
+
+    Args:
+        image (Union[sitk.Image, itk.Image]): The image data as :class:`itk.Image` or :class:`SimpleITK.Image`.
+        modality (Union[Modality, str]): The image :class:`~pyradise.data.modality.Modality` or the modality's name.
+    """
+
+    def __init__(self, image: Union[sitk.Image, itk.Image], modality: Union[Modality, str], scaling_value: float) -> None:
+
+        # Handle situation where RTDose intensity images are 4D - with a singleton dimensions.
+        if image.GetDimension() == 4:
+            if image.GetSize()[0] == 1:
+                image = image[0, :, :, :]
+            elif image.GetSize()[1] == 1:
+                image = image[:, 0, :, :]
+            elif image.GetSize()[2] == 1:
+                image = image[:, :, 0, :]
+            elif image.GetSize()[3] == 1:
+                image = image[:, :, :, 0]
+
+        # See here for why scaling is needed: https://dicom.innolitics.com/ciods/rt-dose/rt-dose/3004000e
+        image = sitk.Cast(image, sitk.sitkFloat32)
+        image = image * scaling_value
+
+        super().__init__(image, modality)
+
+
 # Preparation for next release
 # class DoseImage(Image):
 #     """A dose image class including a :class:`~pyradise.data.taping.TransformTape`.

pyradise/data/subject.py

@@ -3,7 +3,7 @@
 from warnings import warn
 
 from .annotator import Annotator
-from .image import Image, IntensityImage, SegmentationImage
+from .image import Image, IntensityImage, DoseImage, SegmentationImage
 from .modality import Modality
 from .organ import Organ
 
@@ -385,7 +385,7 @@ def get_images_by_type(self, image_type: type) -> List[Image]:
         Returns:
             List[Image]: A list of all images of the specified type.
         """
-        if image_type == IntensityImage:
+        if image_type == IntensityImage or image_type == DoseImage:
             return self.intensity_images
         elif image_type == SegmentationImage:
             return self.segmentation_images

pyradise/fileio/crawling.py

@@ -14,8 +14,8 @@
 from .modality_config import ModalityConfiguration
 from .series_info import (DicomSeriesImageInfo, DicomSeriesInfo,
                           DicomSeriesRegistrationInfo, DicomSeriesRTSSInfo,
-                          FileSeriesInfo, IntensityFileSeriesInfo,
-                          SegmentationFileSeriesInfo)
+                          DicomSeriesDoseInfo, FileSeriesInfo,
+                          IntensityFileSeriesInfo, SegmentationFileSeriesInfo)
 
 __all__ = ["Crawler", "SubjectFileCrawler", "DatasetFileCrawler", "SubjectDicomCrawler", "DatasetDicomCrawler"]
 
@@ -451,6 +451,28 @@ def _get_rtss_files(paths: Tuple[str, ...]) -> Tuple[str, ...]:
 
         return tuple(rtss_files)
 
+
+    @staticmethod
+    def _get_rtdose_files(paths: Tuple[str, ...]) -> Tuple[str, ...]:
+        """Get all DICOM RTDOSE files in the subject directory.
+
+        Args:
+            paths (Tuple[str, ...]): The DICOM file paths to check if they specify a DICOM RTDOSE file.
+
+        Returns:
+            Tuple[str, ...]: The DICOM RTDOSE file paths.
+        """
+        valid_sop_class_uid = "1.2.840.10008.5.1.4.1.1.481.2"  # RT Structure Set Storage
+
+        rtdose_files = []
+        for path in paths:
+            dataset = load_dataset_tag(path, (Tag(0x0008, 0x0016),))
+
+            if dataset.get("SOPClassUID", None) == valid_sop_class_uid:
+                rtdose_files.append(path)
+
+        return tuple(rtdose_files)
+
     @staticmethod
     def _generate_image_infos(image_paths: Tuple[Tuple[str, ...], ...]) -> Tuple[DicomSeriesImageInfo]:
         """Generate the :class:`~pyradise.fileio.series_info.DicomSeriesImageInfo` entries for the DICOM file paths
@@ -504,7 +526,7 @@ def _generate_rtss_info(rtss_paths: Tuple[str, ...]) -> Tuple[DicomSeriesRTSSInf
             rtss_paths (Tuple[str, ...]): The DICOM RTSS file paths.
 
         Returns:
-            Tuple[DicomSeriesRTStructureSetInfo, ...]: AThe retrieved
+            Tuple[DicomSeriesRTStructureSetInfo, ...]: The retrieved
              :class:`~pyradise.fileio.series_info.DicomSeriesRTStructureSetInfo` entries.
         """
         infos = []
@@ -515,6 +537,27 @@ def _generate_rtss_info(rtss_paths: Tuple[str, ...]) -> Tuple[DicomSeriesRTSSInf
 
         return tuple(infos)
 
+
+    @staticmethod
+    def _generate_rtdose_info(rtdose_paths: Tuple[str, ...]) -> Tuple[DicomSeriesImageInfo]:
+        """Generate the :class:`~pyradise.fileio.series_info.DicomSeriesImageInfo` entries for the DICOM file
+        paths specified.
+
+        Args:
+            rtdose_paths (Tuple[str, ...]): The DICOM RTDOSE file paths.
+
+        Returns:
+            Tuple[DicomSeriesImageInfo, ...]: The retrieved
+             :class:`~pyradise.fileio.series_info.DicomSeriesImageInfo` entries.
+        """
+        infos = []
+
+        for path in rtdose_paths:
+            rtdose_info = DicomSeriesDoseInfo(path)
+            infos.append(rtdose_info)
+
+        return tuple(infos)
+
     def _export_modality_config(self, infos: Tuple[DicomSeriesInfo, ...]) -> None:
         """Export the retrieved :class:`~pyradise.fileio.modality_config.ModalityConfiguration` to a file.
 
@@ -647,16 +690,20 @@ def execute(self) -> Tuple[DicomSeriesInfo, ...]:
         remaining_paths = tuple(set(remaining_paths) - set(registration_paths))
 
         rtss_paths = self._get_rtss_files(remaining_paths)
+        remaining_paths = tuple(set(remaining_paths) - set(rtss_paths))
+
+        rtdose_paths = self._get_rtdose_files(remaining_paths)
 
         # generate the series infos
         image_infos = self._generate_image_infos(image_paths)
         registration_infos = self._generate_registration_infos(registration_paths, image_infos)
         rtss_infos = self._generate_rtss_info(rtss_paths)
+        rtdose_infos = self._generate_rtdose_info(rtdose_paths)
 
         # apply the modality config and write it to disk if requested
-        self._apply_modality_config(image_infos)
+        self._apply_modality_config(image_infos + rtdose_infos)
 
-        return image_infos + registration_infos + rtss_infos
+        return image_infos + registration_infos + rtss_infos + rtdose_infos
 
 
 class DatasetDicomCrawler(Crawler):

pyradise/fileio/dicom_conversion.py

@@ -24,14 +24,14 @@
 from pydicom.uid import (PYDICOM_IMPLEMENTATION_UID, ImplicitVRLittleEndian,
                          generate_uid)
 
-from pyradise.data import (IntensityImage, Modality, Organ, SegmentationImage,
+from pyradise.data import (IntensityImage, Modality, Organ, SegmentationImage, DoseImage,
                            Subject, str_to_modality)
 from pyradise.utils import (chunkify, convert_to_itk_image,
                             get_slice_direction, get_slice_position,
                             get_spacing_between_slices, load_dataset,
                             load_dataset_tag, load_datasets)
 
-from .series_info import (DicomSeriesImageInfo, DicomSeriesRegistrationInfo,
+from .series_info import (DicomSeriesImageInfo, DicomSeriesDoseInfo, DicomSeriesRegistrationInfo,
                           DicomSeriesRTSSInfo, RegistrationInfo, SeriesInfo)
 
 __all__ = [
@@ -2594,7 +2594,10 @@ def convert(self) -> Tuple[IntensityImage, ...]:
 
             # if no registration info is available, the image is added as is
             if reg_info is None:
-                image_ = IntensityImage(image, info.modality)
+                if isinstance(info, DicomSeriesDoseInfo):
+                    image_ = DoseImage(image, info.modality, info.scaling_value)
+                else:
+                    image_ = IntensityImage(image, info.modality)
                 image_.add_data({"SeriesInstanceUID": info.series_instance_uid})
                 images.append(image_)
 
@@ -2610,7 +2613,10 @@ def convert(self) -> Tuple[IntensityImage, ...]:
                     )
 
                 image = self._transform_image(image, reg_info.transform, is_intensity=True)
-                image_ = IntensityImage(image, info.modality)
+                if isinstance(info, DicomSeriesDoseInfo):
+                    image_ = DoseImage(image, info.modality, info.scaling_value)
+                else:
+                    image_ = IntensityImage(image, info.modality)
                 image_.add_data({"SeriesInstanceUID": info.series_instance_uid})
                 images.append(image_)
 

pyradise/fileio/series_info.py

@@ -21,6 +21,7 @@
     "SegmentationFileSeriesInfo",
     "DicomSeriesInfo",
     "DicomSeriesImageInfo",
+    "DicomSeriesDoseInfo",
     "DicomSeriesRegistrationInfo",
     "DicomSeriesRTSSInfo",
     "ReferenceInfo",
@@ -417,6 +418,23 @@ def update(self) -> None:
         self._is_updated = True
 
 
+class DicomSeriesDoseInfo(DicomSeriesImageInfo):
+    """A :class:`DicomSeriesDoseInfo` class for DICOM Dose images. In addition to the information provided by the
+    :class:`DicomSeriesImageInfo` class, this class contains a flag to indicate the image is a Dose volume.
+
+    Args:
+        paths (Tuple[str, ...]): The paths to the DICOM image files to load.
+    """
+
+    def __init__(self, paths: Tuple[str, ...]) -> None:
+        super().__init__(paths)
+        scaling_tag = [Tag(0x3004,0x000E)]
+        dataset = load_dataset_tag(self.path[0], scaling_tag)
+
+        self.scaling_value = str(dataset.get("DoseGridScaling", 1.0))
+        self.is_dose_image = True
+
+
 # noinspection PyUnresolvedReferences
 @dataclass
 class ReferenceInfo: