improc ass4: Changed gD() function to use the Gauss separability property.

improc/ass4/gauss.py

 #!/usr/bin/env python
-from numpy import zeros, arange, meshgrid, array, dot
+from numpy import zeros, arange, meshgrid, array, matrix
 from math import ceil, exp, pi, sqrt
 from matplotlib.pyplot import imread, imshow, plot, xlabel, ylabel, show, \
         subplot, xlim, savefig
@@ -52,38 +52,33 @@ def Gauss1(s, order=0):
     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 - r, s) for x in xrange(size)])
 
-    # Make sure that the sum of all kernel values is equal to one
     if not order:
+        # Make sure that the sum of all kernel values is equal to one
         W /= W.sum()
 
     return W
 
-def plot_mask(W, ax):
-    """"""
-    x = arange(W.shape[0])
-    Y, X = meshgrid(x, x)
-    #stride = s / 4
-    #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('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."""
     Fy = Gauss1(s, iorder)
     Fx = Fy if jorder == iorder else Gauss1(s, jorder)
-    W = dot(array([Fy]).T, array([Fx]))
+    G = convolve1d(F, Fy, axis=0, mode='nearest')
 
-    return convolve(F, W, mode='nearest')
+    return convolve1d(G, Fx, axis=1, mode='nearest')
+
+def plot_kernel(W, ax):
+    """Create a 3D plot of a kernel."""
+    x = arange(W.shape[0])
+    Y, X = meshgrid(x, x)
+    ax.plot_surface(X, Y, array(W), rstride=1, cstride=1, linewidth=0, \
+                    antialiased=True, cmap='jet')
+    ax.set_xlabel('y')
+    ax.set_ylabel('x')
+    ax.set_zlabel('g(x, y)')
 
 if __name__ == '__main__':
     if len(argv) < 2:
@@ -95,18 +90,21 @@ if __name__ == '__main__':
         if len(argv) < 5:
             exit_with_usage()
 
+        # Calculate the gaussian kernel using derivatives of the specified
+        # order in both directions
         s = float(argv[2])
         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]))
+        Fy = matrix([Gauss1(s, iorder)])
+        Fx = Fy if jorder == iorder else matrix([Gauss1(s, jorder)])
+        W = Fy.T * Fx
         G = gD(F, s, iorder, jorder)
 
+        # Show the original image, kernel and convoluted image respectively
         subplot(131)
         imshow(F, cmap='gray')
-        plot_mask(W, subplot(132, projection='3d'))
+        plot_kernel(W, subplot(132, projection='3d'))
         subplot(133)
         imshow(G, cmap='gray')
     elif argv[1] == 'timer':
@@ -120,14 +118,17 @@ if __name__ == '__main__':
         S = [1, 2, 3, 5, 7, 9, 11, 15, 19]
         times = []
 
-        for i, s in enumerate(S):
+        for s in S:
             t = 0
 
+            # Average a number of timings to eliminate noise
             for k in xrange(repeat):
                 start = time()
 
                 if method == '1d':
-                    convolve1d(F, Gauss1(s), axis=0, mode='nearest')
+                    W = Gauss1(s)
+                    G = convolve1d(F, W, axis=0, mode='nearest')
+                    convolve1d(G, W, axis=1, mode='nearest')
                 elif method == '2d':
                     convolve(F, Gauss(s), mode='nearest')
 
@@ -148,14 +149,10 @@ if __name__ == '__main__':
         W = Gauss(s)
         G = convolve(F, W, mode='nearest')
 
-        # Original image
+        # Show the original image, kernel and convoluted image respectively
         subplot(131)
         imshow(F, cmap='gray')
-
-        # Gauss function (3D plot)
-        plot_mask(W, subplot(132, projection='3d'))
-
-        # Convolution
+        plot_kernel(W, subplot(132, projection='3d'))
         subplot(133)
         imshow(G, cmap='gray')