Started to add drusen computation from layer segmentation.

eyepy/core/octbase.py

@@ -9,6 +9,10 @@
 
 class Oct(ABC):
 
+    @abstractmethod
+    def __init__(self):
+        self._drusen = None
+
     @abstractmethod
     def __getitem__(self, key):
         raise NotImplementedError()
@@ -47,6 +51,12 @@ def segmentation_raw(self):
     def meta(self):
         raise NotImplementedError()
 
+    @property
+    def drusen(self):
+        if self._drusen is None:
+            self._drusen = np.stack([x.drusen for x in self._bscans], axis=1)
+        return self._drusen
+
     def plot(self, ax=None, bscan_positions=None, line_kwargs=None):
         """ Plot Slo with localization of corresponding B-Scans"""
 
@@ -116,6 +126,166 @@ def segmentation(self):
     def segmentation_raw(self):
         raise NotImplementedError()
 
+    @property
+    @abstractmethod
+    def drusen(self):
+        """ Compute drusen from the RPE layer segmentation.
+
+        First estimate the normal RPE by fitting a polynomial to the RPE. Then
+        compute drusen as the area between the RPE and the normal RPE
+        """
+        rpe_height = self.segmentation["RPE"]
+        bm_height = self.segmentation["BM"]
+        normal_rpe_height = normal_RPE_estimation()
+
+        # def normal_RPE_estimation(layer, degree=3, it=3, s_ratio=1, \
+        #                          farDiff=5, ignoreFarPoints=True, returnImg=False, \
+        #                          useBM=False, useWarping=True, xloc=[], yloc=[], polyFitType="Reguilarized"):
+        """
+        Given the RPE layer, estimate the normal RPE. By first warping the RPE
+        layer with respect to the BM layer. Then fit a third degree polynomial
+        on the RPE layer, and warp the resulting curve back.
+        """
+        # if (useWarping):
+
+        #    yn, xn = warp_BM(layer)
+        #    return yn, xn
+
+        # def warp_BM(seg_img, returnWarpedImg=False):
+        #    """
+        #    Warp seg_img with respect to the BM layer in the segmentation image.
+        #    Return the location of the warped RPE layer in the end. If
+        #    returnWarpedImg is set to True, return the whole warped seg_img.
+        #    """
+        h, w = self.scan.shape
+        # yr, xr = get_RPE_location(seg_img)
+        # yb, xb = get_BM_location(seg_img)
+        rmask = np.zeros((h, w), dtype='int')
+        bmask = np.zeros((h, w), dtype='int')
+
+        rmask[yr, xr] = 255
+        bmask[yb, xb] = 255
+
+        vis_img = np.copy(seg_img)
+
+        wvector = np.empty((w), dtype='int')
+        wvector.fill(h - (h / 2))
+
+        zero_x = []
+        zero_part = False
+        last_nonzero_diff = 0
+        for i in range(w):
+            # get BM height(s) in A-Scan
+            bcol = np.where(bmask[:, i] > 0)[0]
+            # Set to half the image height - largest BM height in this column
+            # (we have only one height per column) This is how much the BM has
+            # to get shifted to align to the horizontal center line
+            wvector[i] = wvector[i] - np.max(bcol) if len(bcol) > 0 else 0
+
+            # In case there is no BM given in this column
+            if (len(bcol) == 0):
+                # Keep track of the gaps in the BM
+                zero_part = True
+                zero_x.append(i)
+            if (len(bcol) > 0 and zero_part):
+                # Set the gaps found sofar to the current BM diff
+                diff = wvector[i]
+                zero_part = False
+                wvector[zero_x] = diff
+                zero_x = []
+            if (len(bcol) > 0):
+                last_nonzero_diff = wvector[i]
+            if (i == w - 1 and zero_part):
+                # If we are done and there are open gaps fill it with last BM diff value
+                wvector[zero_x] = last_nonzero_diff
+
+        shifted = np.zeros(vis_img.shape)
+        nrmask = np.zeros((h, w), dtype='int')
+        nbmask = np.zeros((h, w), dtype='int')
+        # Shift BM, RPE and and combination of both
+        for i in range(w):
+            nrmask[:, i] = np.roll(rmask[:, i], wvector[i])
+            nbmask[:, i] = np.roll(bmask[:, i], wvector[i])
+            shifted[:, i] = np.roll(vis_img[:, i], wvector[i])
+        # now shift the RPE location vector as well
+        shifted_yr = []
+        for x, y in enumerate(rpe_height):
+            shifted_yr.append(y[i] + wvector[x])
+
+        # yn, xn = normal_RPE_estimation(rmask, it=5, useWarping=False, \
+        #                               xloc=xr, yloc=shifted_yr)
+        xloc = range(len(rpe_height))
+        yloc = shifted_yr
+        y = yloc
+        x = xloc
+        tmpx = np.copy(x)
+        tmpy = np.copy(y)
+        origy = np.copy(y)
+        origx = np.copy(x)
+        finalx = np.copy(tmpx)
+        finaly = tmpy
+
+        degree = 3
+        it = 5
+        s_ratio = 1
+        farDiff = 5
+        ignoreFarPoints = True
+        polyFitType = "Reguilarized"
+
+        for i in range(it):
+            if (s_ratio > 1):
+                s_rate = len(tmpx) / s_ratio
+                rand = np.random.rand(s_rate) * len(tmpx)
+                rand = rand.astype('int')
+
+                sx = tmpx[rand]
+                sy = tmpy[rand]
+                if (polyFitType == 'None'):
+                    z = np.polyfit(sx, sy, deg=degree)
+                else:
+                    z = self.compute_reguilarized_fit(sx, sy, deg=degree)
+            else:
+                if (polyFitType == 'None'):
+                    z = np.polyfit(tmpx, tmpy, deg=degree)
+                else:
+                    z = self.compute_reguilarized_fit(tmpx, tmpy, deg=degree)
+            p = np.poly1d(z)
+
+            new_y = p(finalx).astype('int')
+            if (ignoreFarPoints):
+                tmpx = []
+                tmpy = []
+                for i in range(0, len(origx)):
+                    diff = new_y[i] - origy[i]
+                    if diff < farDiff:
+                        tmpx.append(origx[i])
+                        tmpy.append(origy[i])
+            else:
+                tmpy = np.maximum(new_y, tmpy)
+            finaly = new_y
+
+        yn, xn = finaly, finalx
+
+        for i in range(len(xn)):
+            yn[i] = yn[i] - wvector[xn[i]]
+
+        return yn, xn
+
+        # Create drusen map
+        drusen = np.zeros(self.scan.shape)
+        # Exclude normal RPE and RPE from the drusen area.
+        drusen[normal_rpe_height + 1:rpe_height] = 1
+        return drusen
+
+    def compute_reguilarized_fit(self, x, y, deg):
+        resMat = np.zeros((deg + 1, deg + 1))
+        for d in range(deg + 1):
+            z = np.polyfit(x, y, deg=d)
+            for i in range(len(z)):
+                resMat[d, -1 - i] = z[-1 - i]
+        weightedAvg = np.average(resMat, axis=0, weights=[1., 1., 0.1 * 2, 0.1 ** 4])
+        return weightedAvg
+
     def plot(self, ax=None, layers=None, layers_kwargs=None, layers_color=None, layers_only=False):
         """ Plot B-Scan with segmented Layers """
         if ax is None:
@@ -137,7 +307,7 @@ def plot(self, ax=None, layers=None, layers_kwargs=None, layers_color=None, laye
             layers_color = config.layers_color
         else:
             layers_color = {**config.layers_color, **layers_color}
-        
+
         if not layers_only:
             ax.imshow(self.scan, cmap="gray")
         for layer in layers:

eyepy/io/heyex/heyex.py

@@ -62,6 +62,7 @@ def __init__(self, bscans, slo, meta):
         slo :
         meta :
         """
+        super().__init__()
         self._bscans = bscans
         self._sloreader = slo
         self._meta = meta