Work in progress on the DrusenFinder

eyepy/core/drusen.py

@@ -1,5 +1,54 @@
+from abc import ABC, abstractmethod
+
 import numpy as np
 from scipy.interpolate import interp1d
+from scipy.ndimage import label
+
+
+class DrusenFinder(ABC):
+    @abstractmethod
+    def find(self):
+        """ A function which returns a boolean map indicating drusen"""
+        raise NotImplementedError()
+
+    @abstractmethod
+    def filter(self, drusen_map):
+        """ A function which filters drusen based on the drusen properties """
+        raise NotImplementedError()
+
+
+class DefaultDrusenFinder(DrusenFinder):
+    def __init__(self, degree=3, iterations=5, outlier_threshold=5,
+                 poly_fit_type="regularized", minimum_height=2,
+                 minimum_area=0, minimum_volume=0, voxel_size=(1, 1, 1)):
+        """
+
+        Parameters
+        ----------
+        degree: Maximum degree of the Polynomial used for finding the normal RPE
+        iterations: Number of iterations in the outlier removal
+        outlier_threshold: Threshold for RPE Pixel to be considered in the next iteration
+        poly_fit_type: If "regularized", regularize the normal RPE fitting.
+        minimum_height: Minimum drusen  height.
+        minimum_area: Minimum enface area. Currently not used
+        minimum_volume: Minimum drusen volume. Currently not used
+        voxel_size: (height, width, depth) length of the voxel in µm.
+        """
+        self.degree = degree
+        self.iterations = iterations
+        self.outlier_threshold = outlier_threshold
+        self.poly_fit_type = poly_fit_type
+        self.minimum_height = minimum_height
+        self.minimum_volume = minimum_volume
+        self.minimum_area = minimum_area
+        self.voxel_size = voxel_size
+
+    def filter(self, drusen_map):
+        return filter_by_height(drusen_map, self.minimum_height, self.voxel_size)
+
+    def find(self, rpe_height, bm_height, scan_shape):
+        return drusen(rpe_height, bm_height, scan_shape, degree=self.degree, iterations=self.iterations,
+                      outlier_threshold=self.outlier_threshold, poly_fit_type=self.poly_fit_type)
 
 
 def drusen(rpe_height, bm_height, scan_shape, degree=3, iterations=5,
@@ -22,11 +71,8 @@ def drusen(rpe_height, bm_height, scan_shape, degree=3, iterations=5,
     # Exclude normal RPE and RPE from the drusen area.
     for col in range(drusen_map.shape[1]):
         drusen_map[(rpe_height + 1).astype(int)[col]: normal_rpe_height.astype(int)[col], col] = 1
-    # drusen_ranges = [np.s_[start, stop] for start, stop in zip(normal_rpe_height.astype(int) + 1, rpe_height.astype(int))]
-    # drusen_map[drusen_ranges, :] = 1
     return drusen_map
 
-
 def interpolate_layer(layer_height, kind="nearest"):
     nans = np.isnan(layer_height)
     x = np.arange(layer_height.shape[0])[~nans]
@@ -40,7 +86,6 @@ def interpolate_layer(layer_height, kind="nearest"):
 def normal_rpe(rpe_height, bm_height, scan_shape, degree=3, iterations=5,
                outlier_threshold=5, poly_fit_type="regularized"):
     """ Estimate the normal RPE
-
     First the shift to make the BM a horizontal line is found. Then this shift
     is applied to the RPE. A third degree polynomial is fitted to the shifted
     RPE and the resulting normal RPE curve is shifted back into the original
@@ -101,3 +146,28 @@ def compute_regularized_fit(x, y, deg):
     weighted_average = np.average(result_matrix, axis=0,
                                   weights=[1., 1., 0.1 * 2, 0.1 ** 4])
     return weighted_average
+
+
+def filter_by_height(drusen_map, minimum_height=2):
+    if minimum_height == 0:
+        return drusen_map
+    connected_component_array, num_drusen = label(drusen_map)
+    height_map = np.sum(drusen_map, axis=0)
+    component_height_array = component_max_height(connected_component_array, height_map)
+
+    filtered_drusen = np.copy(drusen_map)
+    filtered_drusen[component_height_array <= minimum_height] = 0
+    return filtered_drusen
+
+
+def component_max_height(connected_component_array, height_map):
+    labels = np.unique(connected_component_array)
+    max_heights = np.zeros_like(connected_component_array)
+    for component_label in labels:
+        region = connected_component_array == component_label
+        max_heights[region] = np.max(height_map[np.sum(region, axis=0) > 0])
+    return max_heights
+
+
+def filter_by_area():
+    pass

eyepy/core/octbase.py

@@ -5,14 +5,16 @@
 from matplotlib import pyplot as plt
 
 from eyepy.core import config
-from eyepy.core.drusen import drusen
+from eyepy.core.drusen import DefaultDrusenFinder
 
 
 class Oct(ABC):
 
     @abstractmethod
-    def __init__(self):
+    def __init__(self, drusenfinder=DefaultDrusenFinder()):
+        self.drusenfinder = drusenfinder
         self._drusen = None
+        self._drusen_raw = None
 
     @abstractmethod
     def __getitem__(self, key):
@@ -58,6 +60,12 @@ def drusen(self):
             self._drusen = np.stack([x.drusen for x in self._bscans], axis=-1)
         return self._drusen
 
+    @property
+    def drusen_raw(self):
+        if self._drusen_raw is None:
+            self._drusen_raw = np.stack([x.drusen_raw for x in self._bscans], axis=-1)
+        return self._drusen_raw
+
     def plot(self, ax=None, bscan_positions=None, line_kwargs={"linewidth": 0.5, "color": "green"}):
         """ Plot Slo with localization of corresponding B-Scans"""
 
@@ -108,6 +116,11 @@ def plot_bscans(self, bs_range=range(0, 8), cols=4, layers=None,
 
 class Bscan:
 
+    @abstractmethod
+    def __init__(self):
+        self._drusen = None
+        self._drusen_raw = None
+
     @property
     @abstractmethod
     def scan(self):
@@ -128,12 +141,22 @@ def segmentation(self):
     def segmentation_raw(self):
         raise NotImplementedError()
 
+    @property
+    @abstractmethod
+    def drusen_raw(self):
+        """ Return drusen computed from the RPE and BM layer segmentation, not filtered """
+        if self._drusen_raw is None:
+            self._drusen_raw = self.drusenfinder.find(self.segmentation["RPE"], self.segmentation["BM"],
+                                                      self.scan.shape)
+        return self._drusen_raw
+
     @property
     @abstractmethod
     def drusen(self):
-        """Return drusen computed from the RPE and BM layer segmentation"""
-        return drusen(self.segmentation["RPE"], self.segmentation["BM"],
-                      self.scan.shape)
+        """ Return filtered drusen """
+        if self._drusen is None:
+            self._drusen = self.drusenfinder.filter(self.drusen_raw)
+        return self._drusen
 
     def plot(self, ax=None, layers=None, drusen=False, layers_kwargs=None, layers_color=None,
              annotation_only=False):