Drusen code is clean, tests missing

eyepy/core/drusen.py

@@ -0,0 +1,91 @@
+import numpy as np
+from scipy.interpolate import interp1d
+
+
+def drusen(rpe_height, bm_height, scan_shape, degree=3, iterations=5,
+           outlier_threshold=5, poly_fit_type="regularized"):
+    """ 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
+    """
+    # Estimate normal RPE
+    normal_rpe_height = normal_rpe(rpe_height, bm_height, scan_shape, degree,
+                                   iterations, outlier_threshold,
+                                   poly_fit_type)
+    # Create drusen map
+    drusen_map = np.zeros(scan_shape)
+    # Exclude normal RPE and RPE from the drusen area.
+    drusen_map[normal_rpe_height + 1:rpe_height] = 1
+    return drusen_map
+
+
+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
+    image space.
+    """
+    h, w = scan_shape
+
+    # interpolate NANs in BM
+    nans = np.isnan(bm_height)
+    f = interp1d(np.arange(w)[~nans], bm_height[~nans], kind="nearest")
+    bm_height[nans] = f(np.argwhere(nans))
+
+    # compute shift needed to align the BM to the horizontal center line
+    shift = np.empty(w, dtype='int')
+    shift.fill(h - (h / 2))
+    shift = shift - bm_height
+
+    # now shift the RPE location vector as well
+    shifted_rpe_height = rpe_height + shift
+
+    # These variables change while outliers are removed and polynomial fitting
+    tmpx = range(len(rpe_height))
+    tmpy = np.copy(shifted_rpe_height)
+
+    it = 0
+    while True:
+        if poly_fit_type == 'regularized':
+            coeffs = compute_regularized_fit(tmpx, tmpy, deg=degree)
+        else:
+            coeffs = np.polyfit(tmpx, tmpy, deg=degree)
+
+        # Evaluate the polynomial for all x values
+        norm_rpe = np.polyval(coeffs, range(w)).astype('int')
+
+        # Outlier removal
+        if outlier_threshold:
+            # Compute norm rpe with these values in next iteration
+            inlier = norm_rpe - shifted_rpe_height < outlier_threshold
+            tmpx = np.argwhere(inlier)
+            tmpy = shifted_rpe_height[inlier]
+        else:
+            # TODO: Is this correct
+            # Element wise maximum of the polynomial fitted rpe and the rpe
+            tmpy = np.maximum(norm_rpe, tmpy)
+
+        it += 1
+        if it >= iterations:
+            break
+
+    # Shift back into original image space
+    norm_rpe = norm_rpe - shift
+
+    return norm_rpe
+
+
+def compute_regularized_fit(x, y, deg):
+    result_matrix = 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)):
+            result_matrix[d, -1 - i] = z[-1 - i]
+    # The highest degree has the lowest weight
+    weighted_average = np.average(result_matrix, axis=0,
+                                  weights=[1., 1., 0.1 * 2, 0.1 ** 4])
+    return weighted_average

eyepy/core/octbase.py

@@ -5,6 +5,7 @@
 from matplotlib import pyplot as plt
 
 from eyepy.core import config
+from eyepy.core.drusen import drusen
 
 
 class Oct(ABC):
@@ -129,162 +130,9 @@ def segmentation_raw(self):
     @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
+        """Return drusen computed from the RPE and BM layer segmentation"""
+        return drusen(self.segmentation["RPE"], self.segmentation["BM"],
+                      self.scan.shape)
 
     def plot(self, ax=None, layers=None, layers_kwargs=None, layers_color=None, layers_only=False):
         """ Plot B-Scan with segmented Layers """

setup.py

@@ -16,7 +16,8 @@
     "numpy==1.17.3",
     "matplotlib==3.1.1",
     "scikit-image==0.16.2",
-    "untangle==1.1.1"
+    "untangle==1.1.1",
+    "scipy"
 ]
 
 setup_requirements = ["pytest-runner"]