fix(eyepy.core): set axis in all layer height maps to (n_bscans, width)

Layer in the sample were not correctly oriented due to a different axis order in the HeXMLReader

pyproject.toml

@@ -3,7 +3,7 @@ name = "eyepie"
 packages = [
     {include = "eyepy", from = "src"},
 ]
-version = "0.6.8"
+version = "0.7.0"
 description = "The Python package for working with ophthalmological data."
 authors = ["Olivier Morelle <oli4morelle@gmail.com>"]
 license = "MIT"

src/eyepy/__init__.py

@@ -5,6 +5,7 @@
 __email__ = "oli4morelle@gmail.com"
 __version__ = "0.7.0"
 
+from eyepy.core import drusen
 from eyepy.core import EyeBscan
 from eyepy.core import EyeBscanLayerAnnotation
 from eyepy.core import EyeBscanMeta

src/eyepy/core/__init__.py

@@ -8,3 +8,4 @@
 from .eyemeta import EyeEnfaceMeta
 from .eyemeta import EyeVolumeMeta
 from .eyevolume import EyeVolume
+from .utils import drusen

src/eyepy/core/annotations.py

@@ -104,16 +104,16 @@ class EyeVolumeVoxelAnnotation:
     """ """
 
     def __init__(
-        self,
-        volume: EyeVolume,
-        # Type hint for an optional boolean numpy array
-        data: Optional[npt.NDArray[np.bool_]] = None,
-        meta: Optional[dict] = None,
-        radii=(1.5, 2.5),
-        n_sectors=(1, 4),
-        offsets=(0, 45),
-        center=None,
-        **kwargs,
+            self,
+            volume: EyeVolume,
+            # Type hint for an optional boolean numpy array
+            data: Optional[npt.NDArray[np.bool_]] = None,
+            meta: Optional[dict] = None,
+            radii=(1.5, 2.5),
+            n_sectors=(1, 4),
+            offsets=(0, 45),
+            center=None,
+            **kwargs,
     ):
         """
 
@@ -130,7 +130,9 @@ def __init__(
         self.volume = volume
 
         if data is None:
-            self.data = np.full(self.volume.shape, fill_value=False, dtype=bool)
+            self.data = np.full(self.volume.shape,
+                                fill_value=False,
+                                dtype=bool)
         else:
             self.data = data
 
@@ -149,8 +151,7 @@ def __init__(
                 "n_sectors": n_sectors,
                 "offsets": offsets,
                 "center": center,
-            }
-        )
+            })
 
         if "name" not in self.meta:
             self.meta["name"] = "Voxel Annotation"
@@ -296,7 +297,8 @@ def plot(
             divider = make_axes_locatable(ax)
             cax = divider.append_axes("right", size="5%", pad=0.05)
             plt.colorbar(
-                cm.ScalarMappable(colors.Normalize(vmin=vmin, vmax=vmax), cmap=cmap),
+                cm.ScalarMappable(colors.Normalize(vmin=vmin, vmax=vmax),
+                                  cmap=cmap),
                 cax=cax,
             )
 
@@ -343,32 +345,21 @@ def quantification(self):
         return self._quantification
 
     def _quantify(self):
-        enface_voxel_size_ym3 = (
-            self.volume.localizer.scale_x
-            * 1e3
-            * self.volume.localizer.scale_y
-            * 1e3
-            * self.volume.scale_y
-            * 1e3
-        )
-        oct_voxel_size_ym3 = (
-            self.volume.scale_x
-            * 1e3
-            * self.volume.scale_z
-            * 1e3
-            * self.volume.scale_y
-            * 1e3
-        )
+        enface_voxel_size_ym3 = (self.volume.localizer.scale_x * 1e3 *
+                                 self.volume.localizer.scale_y * 1e3 *
+                                 self.volume.scale_y * 1e3)
+        oct_voxel_size_ym3 = (self.volume.scale_x * 1e3 * self.volume.scale_z *
+                              1e3 * self.volume.scale_y * 1e3)
 
         enface_projection = self.enface
 
         results = {}
         for name, mask in self.masks.items():
-            results[f"{name} [mm³]"] = (
-                (enface_projection * mask).sum() * enface_voxel_size_ym3 / 1e9
-            )
+            results[f"{name} [mm³]"] = ((enface_projection * mask).sum() *
+                                        enface_voxel_size_ym3 / 1e9)
 
-        results["Total [mm³]"] = enface_projection.sum() * enface_voxel_size_ym3 / 1e9
+        results["Total [mm³]"] = enface_projection.sum(
+        ) * enface_voxel_size_ym3 / 1e9
         results["Total [OCT voxels]"] = self.projection.sum()
         results["OCT Voxel Size [µm³]"] = oct_voxel_size_ym3
         results["Laterality"] = self.volume.laterality
@@ -451,10 +442,8 @@ def plot_quantification(
         mask_img = np.zeros(self.volume.localizer.shape, dtype=float)[region]
         visible = np.zeros_like(mask_img)
         for mask_name in self.masks.keys():
-            mask_img += (
-                self.masks[mask_name][region]
-                * self.quantification[mask_name + " [mm³]"]
-            )
+            mask_img += (self.masks[mask_name][region] *
+                         self.quantification[mask_name + " [mm³]"])
             visible += self.masks[mask_name][region]
 
         if vmin is None:
@@ -466,7 +455,8 @@ def plot_quantification(
             divider = make_axes_locatable(ax)
             cax = divider.append_axes("right", size="5%", pad=0.05)
             plt.colorbar(
-                cm.ScalarMappable(colors.Normalize(vmin=vmin, vmax=vmax), cmap=cmap),
+                cm.ScalarMappable(colors.Normalize(vmin=vmin, vmax=vmax),
+                                  cmap=cmap),
                 cax=cax,
             )
 
@@ -513,11 +503,11 @@ def data(self):
         Returns:
 
         """
-        return self.eyevolumelayerannotation.data[-(self.index + 1), :]
+        return self.eyevolumelayerannotation.data[self.index, :]
 
     @data.setter
     def data(self, value):
-        self.eyevolumelayerannotation.data[-(self.index + 1), :] = value
+        self.eyevolumelayerannotation.data[self.index, :] = value
 
     @property
     def knots(self):
@@ -554,7 +544,9 @@ def __init__(
         self.enface = enface
 
         if data is None:
-            self.data = np.full(self.enface.shape, fill_value=False, dtype=bool)
+            self.data = np.full(self.enface.shape,
+                                fill_value=False,
+                                dtype=bool)
         else:
             self.data = data
 

src/eyepy/core/eyevolume.py

@@ -569,7 +569,7 @@ def add_layer_annotation(self, height_map=None, meta=None, **kwargs):
         """
 
         Args:
-            height_map: Height in shape (n_Bscans, Bscan_width)
+            height_map: Height in shape (n_Bscans, Bscan_width) The first index refers to the bottom most B-scan
             meta: name, current_color, and knots
             **kwargs:
 

src/eyepy/core/utils.py

@@ -135,8 +135,8 @@ def drusen(rpe_height, bm_height, volume_shape, minimum_height=2):
     # Create drusen map
     drusen_map = np.zeros(volume_shape, dtype=bool)
     # Exclude normal RPE and RPE from the drusen area.
-    rpe = np.flip(np.rint(rpe_height + 1).astype(int), axis=0)
-    irpe = np.flip(np.rint(idealrpe).astype(int), axis=0)
+    rpe = np.rint(rpe_height + 1).astype(int)
+    irpe = np.rint(idealrpe).astype(int)
     for sli in range(drusen_map.shape[0]):
         for col in range(drusen_map.shape[2]):
             if not rpe[sli, col] == -9223372036854775808:

src/eyepy/io/he/e2e_reader.py

@@ -226,7 +226,7 @@ def layers(self):
         for i in idset:
             layer = np.full((n_bscans, size_x), np.nan)
             for sl in range(n_bscans):
-                layer[i, :] = layer_segmentations[i, sl]
+                layer[sl, :] = layer_segmentations[i, sl]
 
             layers[i] = layer
 

src/eyepy/io/he/vol_reader.py

@@ -204,7 +204,7 @@ def volume(self) -> EyeVolume:
     @property
     def layers(self):
         return np.stack(
-            [b.layer_segmentations for b in self.parsed_file.bscans], axis=1)
+            [b.layer_segmentations for b in self.parsed_file.bscans], axis=0)
 
     @property
     def localizer_meta(self) -> EyeEnfaceMeta:

src/eyepy/io/he/xml_reader.py

@@ -317,15 +317,15 @@ def volume(self) -> EyeVolume:
         data = np.stack(bscans, axis=0)
 
         layer_height_maps = {
-            name: np.full((data.shape[2], data.shape[0]),
+            name: np.full((data.shape[0], data.shape[2]),
                           np.nan,
                           dtype=np.float32)
             for name in layer_heights
         }
 
         for name, heights in layer_heights.items():
             for index, layer_height in heights:
-                layer_height_maps[name][:, index] = layer_height
+                layer_height_maps[name][index, :] = layer_height
 
         localizer = self.localizer
         volume_meta = self.meta