added drusen depth filter

eyepy/core/drusen.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 from scipy.interpolate import interp1d
-from scipy.ndimage import label
+import scipy.ndimage as ndimage
 
 
 class DrusenFinder(ABC):
@@ -158,7 +158,7 @@ def compute_regularized_fit(x, y, deg):
 def filter_by_height(drusen_map, minimum_height=2, voxel_size=(1, 1, 1)):
     if minimum_height == 0:
         return drusen_map
-    connected_component_array, num_drusen = label(drusen_map)
+    connected_component_array, num_drusen = ndimage.label(drusen_map)
     height_map = np.sum(drusen_map, axis=0)
     component_height_array = component_max_height(connected_component_array,
                                                   height_map)
@@ -168,6 +168,26 @@ def filter_by_height(drusen_map, minimum_height=2, voxel_size=(1, 1, 1)):
     return filtered_drusen.astype(bool)
 
 
+def filter_by_depth(drusen_map, minimum_depth=2):
+    filtered_drusen = np.copy(drusen_map)
+    if minimum_depth == 0:
+        return drusen_map
+    # get array where connected components get same label
+    connected_component_array, num_drusen = ndimage.label(drusen_map)
+    drusen_positions = ndimage.find_objects(connected_component_array)
+    # Go through each component, sum it along 2 axis and check max depth against threshold
+    for i, label in enumerate(range(1, num_drusen+1)):
+        druse = connected_component_array[drusen_positions[i]]
+        druse[druse==label] = 1
+        druse[druse!=label] = 0
+        drusen_depth = np.sum(druse, axis=2)
+        if np.max(drusen_depth) <= minimum_depth:
+            # Remove drusen for this label
+            filtered_drusen[connected_component_array==label] = False
+    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)

eyepy/core/octbase.py

@@ -260,7 +260,8 @@ def drusenfinder(self, drusenfinder):
         self._drusenfinder = drusenfinder
 
     def plot(self, ax=None, slo=True, drusen=False, bscan_region=False,
-             bscan_positions=None, masks=False, region=np.s_[...], alpha=1):
+             bscan_positions=None, masks=False, region=np.s_[...],
+             drusen_kwargs={}):
         """
 
         Parameters
@@ -284,7 +285,7 @@ def plot(self, ax=None, slo=True, drusen=False, bscan_region=False,
         if slo:
             self.plot_enface(ax=ax, region=region)
         if drusen:
-            self.plot_drusen(ax=ax, region=region, alpha=alpha)
+            self.plot_drusen(ax=ax, region=region, **drusen_kwargs)
         if bscan_positions is not None:
             self.plot_bscan_positions(ax=ax, bscan_positions=bscan_positions,
                                       region=region,
@@ -381,8 +382,7 @@ def plot_drusen(self, ax=None, region=np.s_[...], cmap="Reds",
             vmin = 1
 
         visible = np.zeros(drusen[region].shape)
-        visible[np.logical_and(vmin < drusen[region],
-                               drusen[region] < vmax)] = 1
+        visible[vmin < drusen[region]] = 1
 
         if cbar:
             divider = make_axes_locatable(ax)
@@ -391,7 +391,7 @@ def plot_drusen(self, ax=None, region=np.s_[...], cmap="Reds",
                 cm.ScalarMappable(colors.Normalize(vmin=vmin, vmax=vmax),
                                   cmap=cmap), cax=cax)
 
-        ax.imshow(drusen[region], alpha=visible * alpha, cmap=cmap, vmin=vmin,
+        ax.imshow(drusen[region], alpha=visible[region] * alpha, cmap=cmap, vmin=vmin,
                   vmax=vmax)
 
     def plot_enface_bscan(self, ax=None, n_bscan=0):
@@ -520,7 +520,7 @@ def plot(self, ax=None, layers=None, drusen=False, layers_kwargs=None,
         if drusen:
             visible = np.zeros(self.drusen.shape)
             visible[self.drusen] = 1.0
-            ax.imshow(self.drusen, alpha=visible, cmap="Reds")
+            ax.imshow(self.drusen[region], alpha=visible[region], cmap="Reds")
         for layer in layers:
             color = layers_color[layer]
             try: