improc ass4: Used linear independency of Gauss function in gD.

improc/ass4/canny.py

@@ -19,10 +19,10 @@ 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)
 
     # Find intensity gradient
-    mask = Gauss1(1, 1)
+    mask = Gauss1(1.4, 1)
     Gx = convolve1d(F, mask, axis=1, mode='nearest')
     Gy = convolve1d(F, mask, axis=0, mode='nearest')
     G = zeros(F.shape)

improc/ass4/gauss.py

 #!/usr/bin/env python
-from numpy import zeros, arange, pi, e, ceil, meshgrid, array
+from numpy import zeros, arange, pi, e, ceil, meshgrid, array, dot
 from matplotlib.pyplot import imread, imshow, plot, xlabel, ylabel, show, \
         subplot, xlim, savefig
 from mpl_toolkits.mplot3d import Axes3D
@@ -29,32 +29,35 @@ def Gauss(s):
     # Make sure that the sum of all kernel values is equal to one
     return W / W.sum()
 
-def gauss(x, s):
+def f_gauss(x, s):
     return e ** -(x ** 2 / (2 * s ** 2)) / (2 * pi * s ** 2)
 
-def gauss_der_1(x, s):
+def f_gauss_der_1(x, s):
     return -x * e ** -(x ** 2 / (2 * s ** 2)) / (2 * pi * s ** 4)
 
-def gauss_der_2(x, s):
+def f_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]
+    f = [f_gauss, f_gauss_der_1, f_gauss_der_2][order]
     s = float(s)
-    size = int(ceil(3 * s))
-    r = 2 * size + 1
-    W = zeros(r)
+    r = int(ceil(3 * s))
+    size = 2 * r + 1
+    W = zeros(size)
     #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([f(x - size, s) for x in xrange(r)])
+    W = array([f(x - r, s) for x in xrange(size)])
 
     # Make sure that the sum of all kernel values is equal to one
-    return W / W.sum()
+    if not order:
+        W /= W.sum()
+
+    return W
 
 def plot_mask(W, ax):
     """"""
@@ -64,30 +67,17 @@ def plot_mask(W, ax):
     #ax.plot_surface(X, Y, W, rstride=stride, cstride=stride, cmap='jet')
     ax.plot_surface(X, Y, W, rstride=1, cstride=1, linewidth=0, \
             antialiased=True, cmap='jet')
-    ax.set_xlabel('x')
-    ax.set_ylabel('y')
+    ax.set_xlabel('y')
+    ax.set_ylabel('x')
     ax.set_zlabel('g(x, y)')
 
 def gD(F, s, iorder, jorder):
     """Create the Gaussian derivative convolution of image F."""
-    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))
-    W = zeros((r, r))
-    Fx = funcs[iorder]
-    Fy = funcs[jorder]
-
-    for x in iterator:
-        for y in iterator:
-            W[x, y] = Fx(x - size, s) * Fy(y - size, s)
+    Fy = Gauss1(s, iorder)
+    Fx = Fy if jorder == iorder else Gauss1(s, jorder)
+    W = dot(array([Fy]).T, array([Fx]))
 
-    return W, convolve(F, W, mode='nearest')
+    return convolve(F, W, mode='nearest')
 
 if __name__ == '__main__':
     if len(argv) < 2:
@@ -100,7 +90,14 @@ if __name__ == '__main__':
             exit_with_usage()
 
         s = float(argv[2])
-        W, G = gD(F, s, int(argv[3]), int(argv[4]))
+        iorder = int(argv[3])
+        jorder = int(argv[4])
+
+        Fy = Gauss1(s, iorder)
+        Fx = Fy if jorder == iorder else Gauss1(s, jorder)
+        W = dot(array([Fy]).T, array([Fx]))
+        G = gD(F, s, iorder, jorder)
+
         subplot(131)
         imshow(F, cmap='gray')
         plot_mask(W, subplot(132, projection='3d'))