fix(vol_reader.py): get number of read layers from file context inste…

…ad of assuming 17 layers.

src/eyepy/io/__init__.py

@@ -253,6 +253,87 @@ def import_duke_mat(path: Union[str, Path]) -> EyeVolume:
     return volume
 
 
+def import_dukechiu2_mat(path: Union[str, Path]) -> EyeVolume:
+    """Import an OCT volume from the Duke dataset (Chiu_BOE_2014).
+
+    The dataset is available at https://people.duke.edu/~sf59/Chiu_BOE_2014_dataset.htm
+    OCT volumes are stored as .mat files which are parsed by this function and returned as
+    EyeVolume object.
+
+    Args:
+        path: Path to the .mat file
+
+    Returns:
+        Parsed data as EyeVolume object
+    """
+    import scipy.io as sio
+
+    loaded = sio.loadmat(path)
+    volume = np.moveaxis(loaded['images'], -1, 0)
+
+    layer_versions = [
+        'manualLayers1', 'manualLayers2', 'automaticLayersDME',
+        'automaticLayersNormal'
+    ]
+    fluid_versions = ['manualFluid1', 'manualFluid2', 'automaticFluidDME']
+    all_layer_maps = {
+        x: np.moveaxis(loaded[x], -1, 1).astype(np.float32)
+        for x in layer_versions
+    }
+    # Manual annotations are instances of fluid regions, we convert to binary here
+    all_pixel_maps = {
+        x: (np.moveaxis(loaded[x], -1, 0) > 0).astype(bool)
+        for x in fluid_versions
+    }
+
+    bscan_meta = [
+        EyeBscanMeta(
+            start_pos=(0, 0.123 * i),
+            end_pos=(0.01133 * (volume.shape[2] - 1), 0.123 * i),
+            pos_unit='mm',
+        ) for i in range(volume.shape[0] - 1, -1, -1)
+    ]
+    meta = EyeVolumeMeta(
+        scale_x=
+        0.01133,  # This value is the average of min and max values given in the paper.
+        scale_y=0.00387,
+        scale_z=
+        0.123,  # This value is the average of min and max values given in the paper.
+        scale_unit='mm',
+        bscan_meta=bscan_meta,
+    )
+
+    volume = EyeVolume(data=volume, meta=meta)
+    names = {
+        7: 'BM',
+        6: 'OS/RPE',
+        5: 'ISM/ISE',
+        4: 'OPL/ONL',
+        3: 'INL/OPL',
+        2: 'IPL/INL',
+        1: 'NFL/GCL',
+        0: 'ILM'
+    }
+
+    for layer_version in layer_versions:
+        layer_maps = all_layer_maps[layer_version]
+        for i, height_map in enumerate(layer_maps):
+            volume.add_layer_annotation(
+                height_map,
+                name=f'{layer_version}_{names[i]}',
+            )
+
+    for fluid_version in fluid_versions:
+        pixel_map = all_pixel_maps[fluid_version]
+        for i, height_map in enumerate(pixel_map):
+            volume.add_pixel_annotation(
+                pixel_map,
+                name=fluid_version,
+            )
+
+    return volume
+
+
 def import_retouch(path: Union[str, Path]) -> EyeVolume:
     """Import an OCT volume from the Retouch dataset.
 

src/eyepy/io/he/vol_reader.py

@@ -57,11 +57,11 @@
     '__empty' / cs.Padding(168) * 'Spare bytes for future use.',
     # cs.Probe(cs.this.bscan_hdr.sizeof(cs.this)),
     'layer_segmentations' / Segmentations *
-    'Layer segmentations for 17 Layers. If layers are not segmented the data is empty',
-    '__empty' /
-    cs.Padding(cs.this.bscan_hdr_size - 256 - 17 * 4 * cs.this._.size_x) *
+    'Layer segmentations for num_seg Layers. If layers are not segmented the data is empty',
+    '__empty' / cs.Padding(cs.this.bscan_hdr_size - 256 -
+                           cs.this.num_seg * 4 * cs.this._.size_x) *
     "B-scan data starts after 'bscan_hdr_size' bytes from the start of the"
-    'B-scan. This header contains 256 bytes for general information, then 17 '
+    'B-scan. This header contains 256 bytes for general information, then "num_seg" '
     'layer segmentations each composed size_x * 4 bytes. Then we have padding'
     'to fill the bscan_hdr_size.',
     'data' / Bscan * 'B-scan data',
@@ -200,14 +200,27 @@ def volume(self) -> EyeVolume:
 
         layer_height_maps = self.layers
         for name, i in SEG_MAPPING.items():
+            if i >= layer_height_maps.shape[0]:
+                logger.warning(
+                    'The volume contains less layers than expected. The naming might not be correct.'
+                )
+                break
             volume.add_layer_annotation(layer_height_maps[i], name=name)
 
         return volume
 
     @property
     def layers(self):
-        layers = np.stack(
-            [b.layer_segmentations for b in self.parsed_file.bscans], axis=0)
+        la = [b.layer_segmentations for b in self.parsed_file.bscans]
+        n_layers = np.unique([len(l) for l in la])
+        if len(n_layers) > 1:
+            max_layers = np.max(n_layers)
+            la = [
+                np.pad(l, ((0, max_layers - len(l)), (0, 0)),
+                       'constant',
+                       constant_values=np.nan) for l in la
+            ]
+        layers = np.stack(la, axis=0)
         layers[layers >= 3.0e+38] = np.nan
         # Currently the shape is (n_bscans, n_layers, width). Swap the first two axes
         # to get (n_layers, n_bscans, width)
@@ -302,7 +315,8 @@ def bytes(self):
                  grid_type1=0,
                  grid_offset1=0,
                  prog_id='unknown',
-                 localizer=skimage.util.img_as_ubyte(self.volume.localizer.data),
+                 localizer=skimage.util.img_as_ubyte(
+                     self.volume.localizer.data),
                  bscans=self._bscan_dicts))
 
     @property

src/eyepy/io/utils.py

@@ -20,6 +20,7 @@
 
 
 def _get_meta_attr(meta_attr):
+
     def prop_func(self):
         return getattr(self.meta, meta_attr)
 
@@ -134,6 +135,7 @@ def _compute_localizer_oct_transform(
 
 
 def get_date_adapter(construct, epoch, second_frac):
+
     class DateAdapter(cs.Adapter):
 
         def _decode(self, obj, context, path):
@@ -173,7 +175,7 @@ class SegmentationsAdapter(cs.Adapter):
     def _decode(self, obj, context, path):
         return np.ndarray(buffer=obj,
                           dtype='float32',
-                          shape=(17, context._.size_x))
+                          shape=(context.num_seg, context._.size_x))
 
     def _encode(self, obj, context, path):
         return obj.tobytes()
@@ -182,7 +184,8 @@ def _encode(self, obj, context, path):
 IntDate = get_date_adapter(cs.Int64ul, datetime(1601, 1, 1), 1e-7)
 FloatDate = get_date_adapter(cs.Float64l, datetime(1899, 12, 30), 60 * 60 * 24)
 Localizer = LocalizerAdapter(cs.Bytes(cs.this.size_x_slo * cs.this.size_y_slo))
-Segmentations = SegmentationsAdapter(cs.Bytes(17 * cs.this._.size_x * 4))
+Segmentations = SegmentationsAdapter(
+    cs.Bytes(cs.this.num_seg * cs.this._.size_x * 4))
 Bscan = BscanAdapter(cs.Bytes(cs.this._.size_y * cs.this._.size_x * 4))