improc ass4: Added 1d Gaussian function to timer option.

improc/ass4/gauss.py

 #!/usr/bin/env python
-from numpy import zeros, arange, pi, e, ceil, meshgrid
+from numpy import zeros, arange, pi, e, ceil, meshgrid, array
 from matplotlib.pyplot import imread, imshow, plot, xlabel, ylabel, show, \
         subplot, xlim, savefig
 from mpl_toolkits.mplot3d import Axes3D
-from scipy.ndimage import convolve
+from scipy.ndimage import convolve, convolve1d
 from time import time
 from sys import argv, exit
 
 def exit_with_usage():
     """Print an error message with the program's usage and exit the program."""
-    print 'Usage: python %s timer REPEAT | diff SCALE' % argv[0]
+    print 'Usage: python %s timer METHOD [ REPEAT ] | diff SCALE' % argv[0]
     exit(1)
 
 def Gauss(s):
-    """Create a sampled Gaussian function of scale s."""
+    """Sample a two-dimensional Gaussian function of scale s."""
     size = int(ceil(3 * s))
-    r = 2 * size + 1
+    r = 2 * size
     W = zeros((r, r))
-    t = s ** 2.
+    t = float(s) ** 2
     a = 1 / (2 * pi * t)
 
     # Sample the Gaussian function
@@ -28,57 +28,92 @@ def Gauss(s):
     # Make sure that the sum of all kernel values is equal to one
     return W / W.sum()
 
-if len(argv) < 2:
-    exit_with_usage()
-
-F = imread('cameraman.png')
-
-if argv[1] == 'timer':
-    # Time for multiple scales
-    S = [1, 2, 3, 5, 7, 9, 11, 15, 19]
-    repeat = int(argv[2])
-    timings = []
-
-    for i, s in enumerate(S):
-        t = 0
+def Gauss1(s):
+    """Sample a one-dimensional Gaussian function of scale s."""
+    size = int(ceil(3 * s))
+    r = 2 * size
+    W = zeros((r,))
+    t = float(s) ** 2
+    a = 1 / (2 * pi * t)
 
-        for k in xrange(repeat):
-            W = Gauss(s)
-            start = time()
-            convolve(F, W, mode='nearest')
-            t += time() - start
+    # Sample the Gaussian function
+    W = array([a * e ** -((x - size) ** 2 / (2 * t)) for x in xrange(r)])
 
-        timings.append(t / repeat)
+    # Make sure that the sum of all kernel values is equal to one
+    return W / W.sum()
 
-    xlim(S[0], S[-1])
-    xlabel('s')
-    ylabel('time (s)')
-    plot(S, timings, 'o-')
-elif argv[1] == 'diff':
-    # Calculate and plot the convolution of the given scale
-    if len(argv) < 3:
+if __name__ == '__main__':
+    if len(argv) < 2:
         exit_with_usage()
 
-    s = float(argv[2])
-    W = Gauss(s)
-    G = convolve(F, W, mode='nearest')
-
-    # Original image
-    subplot(131)
-    imshow(F, cmap='gray')
-
-    # Gauss function (3D plot)
-    x = arange(W.shape[0])
-    X, Y = meshgrid(x, x)
-    ax = subplot(132, projection='3d')
-    stride = s / 4
-    ax.plot_surface(X, Y, W, rstride=stride, cstride=stride, cmap='jet')
-    ax.set_xlabel('x')
-    ax.set_ylabel('y')
-    ax.set_zlabel('g(x, y)')
-
-    # Convolution
-    subplot(133)
-    imshow(G, cmap='gray')
-
-show()
+    F = imread('cameraman.png')
+
+    #W1 = Gauss1(10)
+    #X = arange(W1.shape[0])
+    #G = convolve1d(F, W1, axis=0, mode='nearest')
+    #subplot(121)
+    #imshow(F, cmap='gray')
+    #subplot(122)
+    #imshow(G, cmap='gray')
+    #show()
+    #exit(0)
+
+    if argv[1] == 'timer':
+        if len(argv) < 3:
+            exit_with_usage()
+
+        method = argv[2]
+
+        # Time for multiple scales
+        S = [1, 2, 3, 5, 7, 9, 11, 15, 19]
+        repeat = int(argv[3]) if len(argv) > 3 else 1
+        n = 0
+        times = []
+
+        for i, s in enumerate(S):
+            t = 0
+
+            for k in xrange(repeat):
+                start = time()
+
+                if method == '1d':
+                    convolve1d(F, Gauss1(s), axis=n, mode='nearest')
+                elif method == '2d':
+                    convolve(F, Gauss(s), mode='nearest')
+
+                t += time() - start
+
+            times.append(t / repeat)
+
+        xlim(S[0], S[-1])
+        xlabel('s')
+        ylabel('time (s)')
+        plot(S, times, 'o-')
+    elif argv[1] == 'diff':
+        # Calculate and plot the convolution of the given scale
+        if len(argv) < 3:
+            exit_with_usage()
+
+        s = float(argv[2])
+        W = Gauss(s)
+        G = convolve(F, W, mode='nearest')
+
+        # Original image
+        subplot(131)
+        imshow(F, cmap='gray')
+
+        # Gauss function (3D plot)
+        x = arange(W.shape[0])
+        Y, X = meshgrid(x, x)
+        ax = subplot(132, projection='3d')
+        stride = s / 4
+        ax.plot_surface(X, Y, W, rstride=stride, cstride=stride, cmap='jet')
+        ax.set_xlabel('x')
+        ax.set_ylabel('y')
+        ax.set_zlabel('g(x, y)')
+
+        # Convolution
+        subplot(133)
+        imshow(G, cmap='gray')
+
+    show()