improc ass4: Added usage of 1d convolution in Canny Edge Detector instead of 2d.

improc/ass4/canny.py

@@ -2,7 +2,8 @@
 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
+from scipy.ndimage import convolve1d
+from gauss import Gauss1, gD
 
 def in_image(p, F):
     """Check if given pixel coordinates p are within the bound of image F."""
@@ -18,20 +19,20 @@ def canny(F, s, Tl=None, Th=None):
     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]
+    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]
+    mask = Gauss1(1, 1)
+    Gx = convolve1d(F, mask, axis=1, mode='nearest')
+    Gy = convolve1d(F, mask, axis=0, mode='nearest')
     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)
+    for y in xrange(F.shape[0]):
+        for x in xrange(F.shape[1]):
+            p = (y, x)
 
-            # Gradient norm and angle
+            # Gradient norm and rounded angle
             G[p] = norm(append(Gx[p], Gy[p]))
             A[p] = int(round(arctan2(Gy[p], Gx[p]) * 4 / pi + 1)) % 4
 

improc/ass4/gauss.py

@@ -29,16 +29,29 @@ def Gauss(s):
     # Make sure that the sum of all kernel values is equal to one
     return W / W.sum()
 
-def Gauss1(s):
+def gauss(x, s):
+    return e ** -(x ** 2 / (2 * s ** 2)) / (2 * pi * s ** 2)
+
+def gauss_der_1(x, s):
+    return -x * e ** -(x ** 2 / (2 * s ** 2)) / (2 * pi * s ** 4)
+
+def gauss_der_2(x, s):
+    return (x ** 2 - s ** 2) * e ** -(x ** 2 / (2 * s ** 2)) \
+            / (2 * pi * s ** 6)
+
+def Gauss1(s, order=0):
     """Sample a one-dimensional Gaussian function of scale s."""
+    f = [gauss, gauss_der_1, gauss_der_2][order]
+    s = float(s)
     size = int(ceil(3 * s))
     r = 2 * size + 1
-    W = zeros((r,))
-    t = float(s) ** 2
-    a = 1 / (2 * pi * t)
+    W = zeros(r)
+    #t = float(s) ** 2
+    #a = 1 / (2 * pi * t)
 
     # Sample the Gaussian function
-    W = array([a * e ** -((x - size) ** 2 / (2 * t)) for x in xrange(r)])
+    #W = array([a * e ** -((x - size) ** 2 / (2 * t)) for x in xrange(r)])
+    W = array([f(x - size, s) for x in xrange(r)])
 
     # Make sure that the sum of all kernel values is equal to one
     return W / W.sum()
@@ -54,11 +67,12 @@ def plot_mask(W, ax, stride):
 
 def gD(F, s, iorder, jorder):
     """Create the Gaussian derivative convolution of image F."""
-    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)) \
-                       / (2 * pi * s ** 6)]
+    funcs = [gauss, gauss_der_1, gauss_der_2]
+    #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)) \
+    #                   / (2 * pi * s ** 6)]
     size = int(ceil(3 * s))
     r = 2 * size + 1
     iterator = map(float, range(r))
@@ -68,7 +82,7 @@ def gD(F, s, iorder, jorder):
 
     for x in iterator:
         for y in iterator:
-            W[x, y] = Fx(x - size) * Fy(y - size)
+            W[x, y] = Fx(x - size, s) * Fy(y - size, s)
 
     return W, convolve(F, W, mode='nearest')