add canny harris

Perception/README.md

@@ -0,0 +1,12 @@
+Overview
+--------
+
+Canny Edge Detection
+--------------------
+
+![canny](Vision/figure/canny.png)
+
+Harris Corner Detection
+-----------------------
+
+![harris](Vision/figure/corner_detection.png)

Perception/Vision/cannny.py

@@ -0,0 +1,111 @@
+"""
+J. Canny. A computational approach to edge detection. DOI:10.1109/TPAMI.1986.4767851
+"""
+
+import numpy as np  
+import matplotlib.pyplot as plt
+from skimage import data, filters, color
+
+
+def gaussian_blur(image, sigma):
+    """Use Gaussian filter to reduce noise
+    """
+    return filters.gaussian(image, sigma)
+
+def calc_gradient(image):
+    """Calculate image gradient
+    """
+    dx = filters.sobel(image,axis=1)
+    dy = filters.sobel(image,axis=0)
+    g_mag = np.hypot(dx, dy)
+    g_direct = np.arctan2(dy, dx)
+    return g_mag, g_direct
+
+def non_maximal_suppresion(g_mag, g_direct):
+    """Compare the pixel gradient magnitude with its neighbor
+    in gradient direction, only the local maximal are preserved
+    """
+    nms = np.zeros_like(g_mag)
+    for y in range(1,g_mag.shape[0]-1):
+        for x in range(1,g_mag.shape[1]-1):
+            if g_mag[y,x]>0:
+                theta = g_direct[y, x]
+                dx, dy = np.cos(theta), np.sin(theta)
+                dx = int(round(dx))
+                dy = int(round(dy))
+                q = g_mag[y+dy, x+dx]
+                r = g_mag[y-dy, x-dx]
+                if g_mag[y,x]>q and g_mag[y,x]>r:
+                    nms[y,x] = g_mag[y,x]
+    return nms
+
+def hysteresis_threshold(image,low,high):
+    """Preserve the image area above high and area above low which are connected to high area.
+
+    Step:
+        1. add pixels above high to a queue
+        2. breadth first search
+    """
+
+    OPEN = []
+    CLOSE = []
+    y_inds, x_inds = np.nonzero(image>high)
+    for y_ind, x_ind in zip(y_inds, x_inds):
+        OPEN.append([y_ind, x_ind])
+
+    acts = [[dy, dx] for dx in range(-1,2) for dy in range(-1,2)]
+    acts.remove([0,0])
+    while len(OPEN)>0:
+        y_ind, x_ind = OPEN.pop(0)
+        if (y_ind, x_ind) not in CLOSE:
+            CLOSE.append([y_ind, x_ind])
+        for dy, dx in acts:
+            v = [y_ind+dy, x_ind+dx]
+            if 0<=v[0]<image.shape[0] and 0<=v[1]<image.shape[1]:
+                if image[v[0],v[1]] > low and v not in CLOSE:
+                    OPEN.append(v)
+    ht = np.zeros_like(image)
+    for y_ind, x_ind in CLOSE:
+        ht[y_ind, x_ind] = image[y_ind, x_ind]
+    return ht
+
+def apply_canny(image, sigma, low, high):
+    if image.shape[-1] == 3: # rgb
+        image = color.rgb2gray(image)
+    smoothed = gaussian_blur(image, sigma=sigma)
+    g_mag, g_direct = calc_gradient(smoothed)
+    nms = non_maximal_suppresion(g_mag, g_direct)
+    ht = hysteresis_threshold(nms, low, high)
+    return ht
+
+def main():
+    sigma = 1
+    low = 0.05
+    high = 0.2
+    image = data.coins()
+    smoothed = gaussian_blur(image, sigma=sigma)
+    g_mag, g_direct = calc_gradient(smoothed)
+    nms = non_maximal_suppresion(g_mag, g_direct)
+    ht = hysteresis_threshold(nms, low, high)
+
+    fig, ax = plt.subplots(nrows=2, ncols=2)
+
+    ax[0, 0].imshow(image, cmap='gray')
+    ax[0, 0].set_title('Original image')
+
+    ax[0, 1].imshow(g_mag, cmap='magma')
+    ax[0, 1].set_title('Gradient')
+
+    ax[1, 0].imshow(nms, cmap='magma')
+    ax[1, 0].set_title('Non maximal suppresion')
+
+    ax[1, 1].imshow(ht>0, cmap='magma')
+    ax[1, 1].set_title('Hysteresis threshold')
+    for a in ax.ravel():
+        a.axis('off')
+
+    plt.tight_layout()
+    plt.show()
+
+if __name__ == "__main__":
+    main()

Perception/Vision/sh_corner.py

@@ -0,0 +1,156 @@
+
+
+"""
+Ref:
+https://github.com/scikit-image/scikit-image/tree/master/skimage/feature
+"""
+
+import numpy as np  
+import matplotlib.pyplot as plt
+from skimage import data, filters, color
+from itertools import combinations_with_replacement
+
+def _compute_derivatives(image, mode='constant', cval=0):
+    """Compute derivatives in axis directions using the Sobel operator.
+    Parameters
+    ----------
+    image : ndarray
+        Input image.
+    mode : {'constant', 'reflect', 'wrap', 'nearest', 'mirror'}, optional
+        How to handle values outside the image borders.
+    cval : float, optional
+        Used in conjunction with mode 'constant', the value outside
+        the image boundaries.
+    Returns
+    -------
+    derivatives : list of ndarray
+        Derivatives in each axis direction.
+    """
+
+    derivatives = [filters.sobel(image, axis=i, mode=mode, cval=cval)
+                   for i in range(image.ndim)]
+
+    return derivatives
+
+
+def structure_tensor(image, sigma=1, mode='constant', cval=0):
+    """Compute structure tensor using sum of squared differences.
+    The (2-dimensional) structure tensor A is defined as::
+        A = [Arr Arc]
+            [Arc Acc]
+    which is approximated by the weighted sum of squared differences in a local
+    window around each pixel in the image. This formula can be extended to a
+    larger number of dimensions (see [1]_).
+    Parameters
+    ----------
+    image : ndarray
+        Input image.
+    sigma : float, optional
+        Standard deviation used for the Gaussian kernel, which is used as a
+        weighting function for the local summation of squared differences.
+    mode : {'constant', 'reflect', 'wrap', 'nearest', 'mirror'}, optional
+        How to handle values outside the image borders.
+    cval : float, optional
+        Used in conjunction with mode 'constant', the value outside
+        the image boundaries.
+    """
+    derivatives = _compute_derivatives(image, mode=mode, cval=cval)
+
+    # structure tensor
+    A_elems = [filters.gaussian(der0 * der1, sigma, mode=mode, cval=cval)
+               for der0, der1 in combinations_with_replacement(derivatives, 2)]
+    return A_elems
+
+def gaussian_blur(image, sigma):
+    """Use Gaussian filter to reduce noise
+    """
+    return filters.gaussian(image, sigma)
+
+
+def corner_harris(image, sigma=1, k=0.05):
+    """Calculate Harris response of the image
+    """
+    Arr, Arc, Acc = structure_tensor(image, sigma=sigma, mode='nearest')
+    # determinant
+    detA = Arr * Acc - Arc ** 2
+    # trace
+    traceA = Arr + Acc
+
+    response = detA - k * traceA ** 2
+    return response
+
+def corner_shi_tomasi(image, sigma=1):
+    """Compute Shi-Tomasi (Kanade-Tomasi) corner measure response image.
+    """
+
+    Arr, Arc, Acc = structure_tensor(image, sigma, mode='nearest')
+
+    # minimum eigenvalue of A
+    response = ((Arr + Acc) - np.sqrt((Arr - Acc) ** 2 + 4 * Arc ** 2)) / 2
+
+    return response
+
+def corner_kitchen_rosenfeld(image, mode='constant', cval=0):
+    """Compute Kitchen and Rosenfeld corner measure response image.
+    The corner measure is calculated as follows::
+        (imxx * imy**2 + imyy * imx**2 - 2 * imxy * imx * imy)
+            / (imx**2 + imy**2)
+    Where imx and imy are the first and imxx, imxy, imyy the second
+    derivatives.
+    """
+
+    imy, imx = _compute_derivatives(image, mode=mode, cval=cval)
+    imxy, imxx = _compute_derivatives(imx, mode=mode, cval=cval)
+    imyy, imyx = _compute_derivatives(imy, mode=mode, cval=cval)
+
+    numerator = (imxx * imy ** 2 + imyy * imx ** 2 - 2 * imxy * imx * imy)
+    denominator = (imx ** 2 + imy ** 2)
+
+    response = np.zeros_like(image, dtype=np.double)
+
+    mask = denominator != 0
+    response[mask] = numerator[mask] / denominator[mask]
+
+    return response
+
+
+def main():
+    from matplotlib import pyplot as plt
+
+    from skimage import data
+    from skimage.transform import warp, AffineTransform
+    from skimage.draw import ellipse
+
+    # Sheared checkerboard
+    tform = AffineTransform(scale=(1.3, 1.1), rotation=1, shear=0.7,
+                            translation=(110, 30))
+    image = warp(data.checkerboard()[:90, :90], tform.inverse,
+                output_shape=(200, 310))
+    # Ellipse
+    rr, cc = ellipse(160, 175, 10, 100)
+    image[rr, cc] = 1
+    # Two squares
+    image[30:80, 200:250] = 1
+    image[80:130, 250:300] = 1
+
+
+    fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(8,6))
+    ax[0,0].imshow(image, cmap='gray')
+    ax[0,0].set_title('Original image')
+
+    ax[0,1].imshow(corner_kitchen_rosenfeld(image), cmap='magma')
+    ax[0,1].set_title('Kitchen')
+
+    ax[1,0].imshow(corner_harris(image), cmap='magma')
+    ax[1,0].set_title('Harris')
+
+    ax[1,1].imshow(corner_shi_tomasi(image), cmap='magma')
+    ax[1,1].set_title('Shi-Tomasi')
+    for a in ax.ravel():
+        a.axis('off')
+
+    plt.tight_layout()
+    plt.show()
+
+if __name__ == "__main__":
+    main()