improc ass4: Fixed Gauss filter.

improc/ass4/gauss.py

@@ -7,51 +7,28 @@ from scipy.ndimage import convolve
 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]
+    exit(1)
+
 def Gauss(s):
     """Create a sampled Gaussian function of scale s."""
     size = int(ceil(3 * s))
-    W = zeros((2 * size + 1, 2 * size + 1))
+    r = 2 * int(ceil(3 * s)) + 1
+    W = zeros((r, r))
     t = s ** 2
     a = 1 / (2 * pi * t)
 
     # Sample the Gaussian function
-    for x in xrange(-size, size + 1):
-        for y in xrange(-size, size + 1):
-            W[x, y] = a * e ** -((x ** 2 + y ** 2) / (2 * t))
+    for x in xrange(r):
+        for y in xrange(r):
+            W[x, y] = a * e ** -(((x - size) ** 2 + (y - size) ** 2) / (2 * t))
 
     # Make sure that the sum of all kernel values is equal to one
     return W / W.sum()
 
-def plot_diff(F, s):
-    """Calculate and plot the convolution of scale s of an image F."""
-    W = Gauss(s)
-    G = convolve(F, W, mode='nearest')
-
-    # Original image
-    subplot(131)
-    imshow(F, cmap='gray')
-
-    # Gauss function
-    x = arange(W.shape[0])
-    X, Y = meshgrid(x, x)
-    ax = subplot(132, projection='3d')
-    ax.plot_surface(X, Y, W, rstride=1, cstride=1, cmap='jet')
-    ax.set_zlim3d(0, 1)
-    ax.set_xlabel('x')
-    ax.set_ylabel('y')
-    ax.set_zlabel('g(x, y)')
-
-    # Convolution
-    subplot(133)
-    imshow(G, cmap='gray')
-
-    show()
-
-def exit_usage():
-    print 'Usage: python %s timer [ FILE ] | diff [ SCALE ]' % argv[0]
-    exit(1)
-
-if len(argv) < 1:
+if len(argv) < 2:
     exit_with_usage()
 
 F = imread('cameraman.png')
@@ -59,7 +36,7 @@ F = imread('cameraman.png')
 if argv[1] == 'timer':
     # Time for multiple scales
     S = [1, 2, 3, 5, 7, 9, 11, 15, 19]
-    repeat = 1
+    repeat = int(argv[2])
     timings = []
 
     for i, s in enumerate(S):
@@ -77,13 +54,29 @@ if argv[1] == 'timer':
     xlabel('s')
     ylabel('time (s)')
     plot(S, timings, 'o-')
-
-    if len(argv) > 2:
-        savefig(argv[2])
-
-    show()
 elif argv[1] == 'diff':
-    if len(argv) < 2:
+    # Calculate and plot the convolution of the given scale
+    if len(argv) < 3:
         exit_with_usage()
 
-    plot_diff(F, int(argv[2]))
+    W = Gauss(int(argv[2]))
+    G = convolve(F, W, mode='nearest')
+
+    # Original image
+    subplot(131)
+    imshow(F, cmap='gray')
+
+    # Gauss function
+    x = arange(W.shape[0])
+    X, Y = meshgrid(x, x)
+    ax = subplot(132, projection='3d')
+    ax.plot_surface(X, Y, W, rstride=1, cstride=1, cmap='jet')
+    ax.set_xlabel('x')
+    ax.set_ylabel('y')
+    ax.set_zlabel('g(x, y)')
+
+    # Convolution
+    subplot(133)
+    imshow(G, cmap='gray')
+
+show()