improc ass4: Implemented Canny Edge Detector.

improc/ass4/canny.py

 #!/usr/bin/env python
-from matplotlib.pyplot import imread, imshow, show
-from numpy import arctan
+from matplotlib.pyplot import imread, imshow, subplot, show
+from numpy import arctan2, zeros, append, pi#, argmax
+from numpy.linalg import norm
 from gauss import gD
 
-def canny(F, s):
-    # Noise reduction by a Gauss filter
-    G = gD(F, s, 2, 2)
+def in_image(p, F):
+    """Check if given pixel coordinates p are within the bound of image F."""
+    return p[0] >= 0 and p[1] >= 0 and p[0] < F.shape[0] and p[1] < F.shape[1]
+
+#def zero_crossing(a, b, F):
+#    """Cech if there is a zero crossing point between F[a] and F[b]."""
+#    return in_image(a, F) and in_image(b, F) \
+#            and ((F[a] < 0 and F[b] > 0) or (F[a] > 0 and F[b] < 0))
+
+def canny(F, s, Tl=None, Th=None):
+    """Apply the Canny Edge Detection algorithm with Gauss scale s to an
+    image F. Optionally specify a low and high threshold (Tl and Th) for
+    hysteresis thresholding."""
+    # Noise reduction by a Gaussian filter
+    #F = gD(F, s, 0, 0)[1]
 
     # Find intensity gradient
+    #F = gD(F, s, 2, 2)[1]
+    Gx = gD(F, s, 1, 0)[1]
+    Gy = gD(F, s, 0, 1)[1]
+    G = zeros(F.shape)
+    A = zeros(F.shape, dtype=int)
+
+    for x in xrange(F.shape[0]):
+        for y in xrange(F.shape[1]):
+            p = (x, y)
+
+            # Gradient norm and angle
+            G[p] = norm(append(Gx[p], Gy[p]))
+            A[p] = int(round(arctan2(Gy[p], Gx[p]) * 4 / pi + 1)) % 4
+
+            #p = (x, y)
+            #compare = [(x, y - 1), (x, y + 1), (x + 1, y - 1), \
+            #           (x - 1, y + 1), (x - 1, y), (x + 1, y), \
+            #           (x - 1, y - 1), (x + 1, y + 1)]
+            #norms = zeros(8)
+
+            #for i, c in enumerate(compare):
+            #    if zero_crossing(p, c, F):
+            #        norms[i] = abs(F[p]) + abs(F[c])
 
-    return G
+            #m = argmax(norms)
+            #G[p] = norms[m]
+            #A[p] = m >> 1
 
     # Non-maximum suppression
+    E = zeros(F.shape)
+
+    for x in xrange(F.shape[0]):
+        for y in xrange(F.shape[1]):
+            g = G[x, y]
+            a = A[x, y]
+            compare = [((x, y - 1), (x, y + 1)), ((x - 1, y - 1), \
+                       (x + 1, y + 1)), ((x - 1, y), (x + 1, y)), \
+                       ((x + 1, y - 1), (x - 1, y + 1))]
+            na, nb = compare[a]
+
+            if (not in_image(na, G) or g > G[na]) \
+                    and (not in_image(nb, G) or g > G[nb]):
+                E[x, y] = g
+
+    # Only execute hysteresis thresholding if the thresholds are specified
+    if Tl is None or Th is None:
+        return E
+
+    # Hysteresis thresholding
+    Tl *= (E.max() - E.min()) / 255
+    Th *= (E.max() - E.min()) / 255
+    T = zeros(F.shape, dtype=bool)
+
+    # Clear image borders
+    for x in xrange(F.shape[0]):
+        E[x, 0] = E[x, F.shape[1] - 1] = 0
 
+    for y in xrange(1, F.shape[1] - 1):
+        E[0, y] = E[F.shape[0] - 1, y] = 0
 
-    # Tracing edges
+    def follow_nb(x, y):
+        """Follow the neighbouring pixels of an edge pixel in E recursively."""
+        if T[x, y]:
+            return
 
+        T[x, y] = True
+
+        for nx in xrange(-1, 2):
+            for ny in xrange(-1, 2):
+                if (nx != 0 or ny != 0) and E[nx, ny] > Tl:
+                    follow_nb(nx, ny)
+
+    # Follow edges that have a starting value above Th
+    for x in xrange(F.shape[0]):
+        for y in xrange(F.shape[1]):
+            if E[x, y] > Th:
+                follow_nb(x, y)
+
+    return E, T
 
 if __name__ == '__main__':
     from sys import argv, exit
 
-    if len(argv) < 2:
-        print 'Usage: python %s SCALE' % argv[0]
+    if len(argv) < 2 or len(argv) == 3:
+        print 'Usage: python %s SCALE [ TL TH ]' % argv[0]
         exit(1)
 
-    s = float(argv[1])
     F = imread('cameraman.png')
+    #F = imread('flags.png')
+
+    s = float(argv[1])
+
+    if len(argv) > 3:
+        # Execute with tracing edges
+        E, T = canny(F, s, float(argv[2]), float(argv[3]))
+
+        subplot(131)
+        imshow(F, cmap='gray')
+        subplot(132)
+        imshow(E, cmap='gray')
+        subplot(133)
+        imshow(T, cmap='gray')
+    else:
+        # Execute until nn-maximum suppression
         E = canny(F, s)
+
+        subplot(121)
+        imshow(F, cmap='gray')
+        subplot(122)
         imshow(E, cmap='gray')
+
     show()

improc/ass4/gauss.py

@@ -57,7 +57,7 @@ def gD(F, s, iorder, jorder):
     funcs = [lambda x: e ** -(x ** 2 / (2 * s ** 2)) / (2 * pi * s ** 2), \
              lambda x: -x * e ** -(x ** 2 / (2 * s ** 2)) \
                        / (2 * pi * s ** 4), \
-             lambda x: -(x ** 2 - s ** 2) * e ** -(x ** 2 / (2 * s ** 2)) \
+             lambda x: (x ** 2 - s ** 2) * e ** -(x ** 2 / (2 * s ** 2)) \
                        / (2 * pi * s ** 6)]
     size = int(ceil(3 * s))
     r = 2 * size + 1