Repository: DIAGNijmegen/rse-panimg (language: python)

Issue:
.e2e file support
Placeholder issue for July cycle pitch: https://github.com/DIAGNijmegen/rse-roadmap/issues/63


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index cf25a01..ef255c7 100644
--- a/README.md
+++ b/README.md
@@ -61,6 +61,7 @@ To access the help test you can use `panimg -h`.
 | `.dcm`                              | `.mha`  | `dicom`    | <sup>[1](#footnote1)</sup> |
 | `.nii`                              | `.mha`  | `nifti`    |                            |
 | `.nii.gz`                           | `.mha`  | `nifti`    |                            |
+| `.e2e`                              | `.mha`  | `oct`      | <sup>[2](#footnote2)</sup> |
 | `.fds`                              | `.mha`  | `oct`      | <sup>[2](#footnote2)</sup> |
 | `.fda`                              | `.mha`  | `oct`      | <sup>[2](#footnote2)</sup> |
 | `.png`                              | `.mha`  | `fallback` | <sup>[3](#footnote3)</sup> |
diff --git a/panimg/image_builders/oct.py b/panimg/image_builders/oct.py
index 58a3434..b3d85da 100644
--- a/panimg/image_builders/oct.py
+++ b/panimg/image_builders/oct.py
@@ -1,6 +1,14 @@
 from collections import defaultdict
 from pathlib import Path
-from typing import DefaultDict, Iterable, Iterator, List, Set, Tuple, Union
+from typing import (
+    DefaultDict,
+    Iterable,
+    Iterator,
+    List,
+    Set,
+    Tuple,
+    Union,
+)
 
 import SimpleITK
 import numpy as np
@@ -42,15 +50,20 @@ class OctDimensions(BaseModel):
 
 def _create_itk_images(
     *,
-    file: Path,
+    file: Tuple[Path, str],
     oct_volumes: Iterable[OCTVolumeWithMetaData],
     fundus_images: Iterable[FundusImageWithMetaData],
     oct_slice_size: OctDimensions,
 ) -> Iterator[SimpleITKImage]:
+
+    filepath = file[0]
+    filetype = file[1]
+
     for volume in oct_volumes:
         eye_choice = LATERALITY_TO_EYE_CHOICE[volume.laterality]
+
         yield _create_itk_oct_volume(
-            file=file,
+            file=filepath,
             volume=volume.volume,
             oct_slice_size=oct_slice_size,
             eye_choice=eye_choice,
@@ -58,7 +71,7 @@ def _create_itk_images(
     for image in fundus_images:
         eye_choice = LATERALITY_TO_EYE_CHOICE[image.laterality]
 
-        if file.suffix != ".e2e":
+        if filetype != "e2e":
             img_array = image.image.astype(np.uint8)
             img_array = img_array[:, :, ::-1]
             is_vector = True
@@ -67,7 +80,7 @@ def _create_itk_images(
             is_vector = False
 
         yield _create_itk_fundus_image(
-            file=file,
+            file=filepath,
             image=img_array,
             eye_choice=eye_choice,
             is_vector=is_vector,
@@ -101,9 +114,10 @@ def _create_itk_oct_volume(
 
 
 def _create_itk_fundus_image(
-    *, file: Path, image: np.ndarray, eye_choice: EyeChoice, is_vector: bool
+    *, file: Path, image: np.ndarray, eye_choice: EyeChoice, is_vector: bool,
 ) -> SimpleITKImage:
     img = SimpleITK.GetImageFromArray(image, isVector=is_vector)
+
     return SimpleITKImage(
         image=img,
         name=file.stem + "_fundus" + file.suffix,
@@ -141,43 +155,50 @@ def _get_image(
     Iterable[OCTVolumeWithMetaData],
     Iterable[FundusImageWithMetaData],
     OctDimensions,
+    Tuple[Path, str],
 ]:
-    if file.suffix == ".fds":
-        fds_img = FDS(file)
-        oct_slice_size = _extract_slice_size(img=fds_img)
-        return (
-            [fds_img.read_oct_volume()],
-            [fds_img.read_fundus_image()],
-            oct_slice_size,
-        )
-    elif file.suffix == ".fda":
-        fda_img = FDA(file)
-        oct_slice_size = _extract_slice_size(img=fda_img)
-        return (
-            [fda_img.read_oct_volume()],
-            [fda_img.read_fundus_image()],
-            oct_slice_size,
-        )
-    elif file.suffix == ".e2e":
-        # support for E2E files is disabled until we have a proper E2E file to
-        # test with
-        raise ValueError
-        e2e_img = E2E(file)
-        # We were unable to retrieve slice size information from the e2e files.
-        # The following default values are taken from:
-        # https://bitbucket.org/uocte/uocte/wiki/Heidelberg%20File%20Format
-        oct_slice_size = OctDimensions(
-            extent_x_mm=6, resolution_y_mm=0.0039, extent_z_mm=4.5
-        )
-        # Note that the return types from oct_converter are different
-        # for e2e files
-        return (
-            e2e_img.read_oct_volume(),
-            e2e_img.read_fundus_image(),
-            oct_slice_size,
-        )
-    else:
-        raise ValueError
+    with open(file, "rb") as f:
+        header = f.read(7)
+
+        if b"FDS" in header:
+            fds_img = FDS(file)
+            file_info = file, "fds"
+            oct_slice_size = _extract_slice_size(img=fds_img)
+            return (
+                [fds_img.read_oct_volume()],
+                [fds_img.read_fundus_image()],
+                oct_slice_size,
+                file_info,
+            )
+        elif b"FDA" in header:
+            fda_img = FDA(file)
+            file_info = file, "fda"
+            oct_slice_size = _extract_slice_size(img=fda_img)
+            return (
+                [fda_img.read_oct_volume()],
+                [fda_img.read_fundus_image()],
+                oct_slice_size,
+                file_info,
+            )
+        elif b"CMD" in header:
+            e2e_img = E2E(file)
+            file_info = file, "e2e"
+            # We were unable to retrieve slice size information from the e2e files.
+            # The following default values are taken from:
+            # https://bitbucket.org/uocte/uocte/wiki/Heidelberg%20File%20Format
+            oct_slice_size = OctDimensions(
+                extent_x_mm=6, resolution_y_mm=0.0039, extent_z_mm=4.5
+            )
+            # Note that the return types from oct_converter are different
+            # for e2e files
+            return (
+                e2e_img.read_oct_volume(),
+                e2e_img.read_fundus_image(),
+                oct_slice_size,
+                file_info,
+            )
+        else:
+            raise ValueError
 
 
 def image_builder_oct(*, files: Set[Path]) -> Iterator[SimpleITKImage]:
@@ -203,10 +224,12 @@ def image_builder_oct(*, files: Set[Path]) -> Iterator[SimpleITKImage]:
 
     for file in files:
         try:
-            oct_volumes, fundus_images, oct_slice_size = _get_image(file=file)
+            oct_volumes, fundus_images, oct_slice_size, file_out = _get_image(
+                file=file
+            )
 
             yield from _create_itk_images(
-                file=file,
+                file=file_out,
                 oct_volumes=oct_volumes,
                 # Skip fundus_images until we decide how to handle these
                 fundus_images=[],  # fundus_images,


Test patch:
diff --git a/panimg/contrib/oct_converter/readers/e2e.py b/panimg/contrib/oct_converter/readers/e2e.py
index d0f7273..f9619ac 100644
--- a/panimg/contrib/oct_converter/readers/e2e.py
+++ b/panimg/contrib/oct_converter/readers/e2e.py
@@ -1,3 +1,6 @@
+import re
+import struct
+
 import numpy as np
 from construct import (
     PaddedString,
@@ -6,7 +9,6 @@
     Int32sn,
     Int32un,
     Int8un,
-    Array,
 )
 from panimg.contrib.oct_converter.image_types import (
     OCTVolumeWithMetaData,
@@ -35,34 +37,7 @@ def __init__(self, filepath):
         self.filepath = Path(filepath)
         if not self.filepath.exists():
             raise FileNotFoundError(self.filepath)
-        self.header_structure = Struct(
-            "magic" / PaddedString(12, "ascii"),
-            "version" / Int32un,
-            "unknown" / Array(10, Int16un),
-        )
-        self.main_directory_structure = Struct(
-            "magic" / PaddedString(12, "ascii"),
-            "version" / Int32un,
-            "unknown" / Array(10, Int16un),
-            "num_entries" / Int32un,
-            "current" / Int32un,
-            "prev" / Int32un,
-            "unknown3" / Int32un,
-        )
-        self.sub_directory_structure = Struct(
-            "pos" / Int32un,
-            "start" / Int32un,
-            "size" / Int32un,
-            "unknown" / Int32un,
-            "patient_id" / Int32un,
-            "study_id" / Int32un,
-            "series_id" / Int32un,
-            "slice_id" / Int32sn,
-            "unknown2" / Int16un,
-            "unknown3" / Int16un,
-            "type" / Int32un,
-            "unknown4" / Int32un,
-        )
+
         self.chunk_structure = Struct(
             "magic" / PaddedString(12, "ascii"),
             "unknown" / Int32un,
@@ -92,82 +67,52 @@ def __init__(self, filepath):
             "unknown2" / Int8un,
         )
 
-    def read_oct_volume(self):
-        """ Reads OCT data.
+    def find_data_chunks(self, f):
+        data = f.read()
+        # find all 'MDbData' chunks
+        regex_pattern = re.compile(b"MDbData")
+        matches = regex_pattern.finditer(data)
+        positions = []
+        for match in matches:
+            positions.append(match.start())
+        return positions
+
+    def extract_laterality_data(self, f):
+        raw = f.read(20)
+        try:
+            laterality_data = self.lat_structure.parse(raw)
+            if laterality_data.laterality == 82:
+                laterality = "R"
+            elif laterality_data.laterality == 76:
+                laterality = "L"
+        except:
+            laterality = None
+        return laterality
 
+    def read_oct_volume(self):
+        """ Reads oct data.
             Returns:
                 obj:OCTVolumeWithMetaData
         """
         with open(self.filepath, "rb") as f:
-            raw = f.read(36)
-            header = self.header_structure.parse(raw)
 
-            raw = f.read(52)
-            main_directory = self.main_directory_structure.parse(raw)
+            chunk_positions = self.find_data_chunks(f)
 
-            # traverse list of main directories in first pass
-            directory_stack = []
-
-            current = main_directory.current
-            while current != 0:
-                directory_stack.append(current)
-                f.seek(current)
-                raw = f.read(52)
-                directory_chunk = self.main_directory_structure.parse(raw)
-                current = directory_chunk.prev
-
-            # traverse in second pass and  get all subdirectories
-            chunk_stack = []
             volume_dict = {}
-            for position in directory_stack:
-                f.seek(position)
-                raw = f.read(52)
-                directory_chunk = self.main_directory_structure.parse(raw)
-
-                for ii in range(directory_chunk.num_entries):
-                    raw = f.read(44)
-                    chunk = self.sub_directory_structure.parse(raw)
-                    volume_string = f"{chunk.patient_id}_{chunk.study_id}_{chunk.series_id}"
-                    if volume_string not in volume_dict.keys():
-                        volume_dict[volume_string] = chunk.slice_id / 2
-                    elif chunk.slice_id / 2 > volume_dict[volume_string]:
-                        volume_dict[volume_string] = chunk.slice_id / 2
-
-                    if chunk.start > chunk.pos:
-                        chunk_stack.append([chunk.start, chunk.size])
-
-            # initalise dict to hold all the image volumes
-            volume_array_dict = {}
-            volume_array_dict_additional = (
-                {}
-            )  # for storage of slices not caught by extraction
-            for volume, num_slices in volume_dict.items():
-                if num_slices > 0:
-                    # num_slices + 1 here due to evidence that a slice was being missed off the end in extraction
-                    volume_array_dict[volume] = [0] * int(num_slices + 1)
-
-            # traverse all chunks and extract slices
-            for start, pos in chunk_stack:
+            for start in chunk_positions:
                 f.seek(start)
                 raw = f.read(60)
                 chunk = self.chunk_structure.parse(raw)
 
                 if chunk.type == 11:  # laterality data
-                    raw = f.read(20)
-                    try:
-                        laterality_data = self.lat_structure.parse(raw)
-                        if laterality_data.laterality == 82:
-                            self.laterality = "R"
-                        elif laterality_data.laterality == 76:
-                            self.laterality = "L"
-                    except:
-                        self.laterality = None
+                    self.laterality = self.extract_laterality_data(f)
 
                 if chunk.type == 1073741824:  # image data
-                    raw = f.read(20)
-                    image_data = self.image_structure.parse(raw)
-
                     if chunk.ind == 1:  # oct data
+                        volume_string = f"{chunk.patient_id}_{chunk.study_id}_{chunk.series_id}"
+                        # read data
+                        raw = f.read(20)
+                        image_data = self.image_structure.parse(raw)
                         all_bits = [
                             f.read(2)
                             for i in range(
@@ -180,40 +125,25 @@ def read_oct_volume(self):
                         image = np.array(raw_volume).reshape(
                             image_data.width, image_data.height
                         )
-                        image = 256 * pow(image, 1.0 / 2.4)
-                        volume_string = f"{chunk.patient_id}_{chunk.study_id}_{chunk.series_id}"
-                        if volume_string in volume_array_dict.keys():
-                            volume_array_dict[volume_string][
-                                int(chunk.slice_id / 2) - 1
-                            ] = image
-                        else:
-                            # try to capture these additional images
-                            if (
-                                volume_string
-                                in volume_array_dict_additional.keys()
-                            ):
-                                volume_array_dict_additional[
-                                    volume_string
-                                ].append(image)
-                            else:
-                                volume_array_dict_additional[volume_string] = [
-                                    image
-                                ]
-                            # print('Failed to save image data for volume {}'.format(volume_string))
+                        normalized_float = pow(image, 1.0 / 2.4) / pow(
+                            2, 1.0 / 2.4
+                        )
+                        image = (normalized_float * (256 * 256 - 1)).astype(
+                            np.uint16
+                        )
+                        if volume_string not in volume_dict:
+                            volume_dict[volume_string] = {}
+                        volume_dict[volume_string][
+                            int(chunk.slice_id / 2)
+                        ] = image
 
             oct_volumes = []
-            for key, volume in volume_array_dict.items():
-                oct_volumes.append(
-                    OCTVolumeWithMetaData(
-                        volume=volume,
-                        patient_id=key,
-                        laterality=self.laterality,
-                    )
-                )
-            for key, volume in volume_array_dict_additional.items():
+            for key, volume in volume_dict.items():
+                slice_order = sorted([id for id in volume.keys()])
+                ordered_volume = [volume[id] for id in slice_order]
                 oct_volumes.append(
                     OCTVolumeWithMetaData(
-                        volume=volume,
+                        volume=ordered_volume,
                         patient_id=key,
                         laterality=self.laterality,
                     )
@@ -228,55 +158,17 @@ def read_fundus_image(self):
                 obj:FundusImageWithMetaData
         """
         with open(self.filepath, "rb") as f:
-            raw = f.read(36)
-            header = self.header_structure.parse(raw)
-
-            raw = f.read(52)
-            main_directory = self.main_directory_structure.parse(raw)
-
-            # traverse list of main directories in first pass
-            directory_stack = []
-
-            current = main_directory.current
-            while current != 0:
-                directory_stack.append(current)
-                f.seek(current)
-                raw = f.read(52)
-                directory_chunk = self.main_directory_structure.parse(raw)
-                current = directory_chunk.prev
-
-            # traverse in second pass and  get all subdirectories
-            chunk_stack = []
-            for position in directory_stack:
-                f.seek(position)
-                raw = f.read(52)
-                directory_chunk = self.main_directory_structure.parse(raw)
-
-                for ii in range(directory_chunk.num_entries):
-                    raw = f.read(44)
-                    chunk = self.sub_directory_structure.parse(raw)
-                    if chunk.start > chunk.pos:
-                        chunk_stack.append([chunk.start, chunk.size])
-
+            chunk_positions = self.find_data_chunks(f)
             # initalise dict to hold all the image volumes
             image_array_dict = {}
-
             # traverse all chunks and extract slices
-            for start, pos in chunk_stack:
+            for start in chunk_positions:
                 f.seek(start)
                 raw = f.read(60)
                 chunk = self.chunk_structure.parse(raw)
 
                 if chunk.type == 11:  # laterality data
-                    raw = f.read(20)
-                    try:
-                        laterality_data = self.lat_structure.parse(raw)
-                        if laterality_data.laterality == 82:
-                            self.laterality = "R"
-                        elif laterality_data.laterality == 76:
-                            self.laterality = "L"
-                    except:
-                        self.laterality = None
+                    self.laterality = self.extract_laterality_data(f)
 
                 if chunk.type == 1073741824:  # image data
                     raw = f.read(20)
diff --git a/tests/resources/oct/e2e_minimized.E2E b/tests/resources/oct/e2e_minimized.E2E
new file mode 100644
index 0000000..be34a90
Binary files /dev/null and b/tests/resources/oct/e2e_minimized.E2E differ
diff --git a/tests/test_oct.py b/tests/test_oct.py
index 342f361..4b1eaaa 100644
--- a/tests/test_oct.py
+++ b/tests/test_oct.py
@@ -70,6 +70,33 @@
                 "eye_choice": "U",
             },
         ),
+        (
+            # Minimized .e2e OCT file was created by taking example OCT file from RUMC
+            # and downsizing:
+            #   - OCT volume (chunk type = 1073741824, chunk ind = 1) to 512x496x1
+            #   - Fundus image (chunk type = 1073741824, chunk ind = 0) to 2x3x1
+            #   - deleting all MDbData chunks other than the OCT slice, the fundus
+            #     image and the laterality information (chunk type = 11)
+            RESOURCE_PATH / "oct/e2e_minimized.E2E",
+            {
+                "width": 2,
+                "height": 3,
+                "depth": None,
+                "voxel_width_mm": None,
+                "voxel_height_mm": None,
+                "voxel_depth_mm": None,
+                "eye_choice": "U",
+            },
+            {
+                "width": 512,
+                "height": 496,
+                "depth": 1,
+                "voxel_width_mm": 0.01171875,
+                "voxel_height_mm": 0.0039,
+                "voxel_depth_mm": 4.5,
+                "eye_choice": "U",
+            },
+        ),
     ),
 )
 def test_image_builder_oct(