Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Newer
Older
100644 127 lines (104 sloc) 4.407 kB
9769e01 @chriskelvinlee gaussian too difficult to implement in gpu for now, using norm smoothing
authored
1 from pylab import *
2 import os
3 import sys
4 import numpy as np
5 import time
6
7 total_start_time = time.time()
8 setup_start_time = time.time()
9
10 #Get the input filename from the command line
11 try:
12 file_name = sys.argv[1]; MaxRad = float(sys.argv[2]); Threshold = float(sys.argv[3])
13 except:
14 print "Usage:",sys.argv[0], "infile maxrad threshold"; sys.exit(1)
15
16 original_image_rgb = imread(file_name)
17
18 # Image is black and white so R=B=G
19 IMG = array( original_image_rgb[:,:,0])
20
21 # Get image data
22 Lx = int32( IMG.shape[0])
23 Ly = int32( IMG.shape[1])
24
25 print "Processing file %s (%d x %d image)" % (file_name, Lx, Ly)
26
27 # Allocate memory
28 # size of the box needed to reach the threshold value or maxrad value
29 BOX = np.zeros((Lx, Ly), dtype=np.float64)
30 # normalized array
31 NORM = np.zeros((Lx, Ly), dtype=np.float64)
32 # output array
33 OUT = np.zeros((Lx, Ly), dtype=np.float64)
34
35 setup_stop_time = time.time()
36 kernel_start_time = time.time()
37
38 # Convolve the image with the gaussian kernel
39 for xx in range(Lx):
40 for yy in range(Ly):
41 qq = 1.0 # var to increase size of box
42 sum = 0.0 # value of the sum
43 ksum = 0.0 # value of the kernal sum
44 ss = qq # size of the box
45
46 # Continue until parameters are met
47 while (sum < Threshold) and (qq < MaxRad):
48 ss = qq
49 # ARE THE FOLLOWING TWO LINES NECESSARY?
50 # SUM = O WILL JUST NULL THE ABOVE THRESHOLD REQUIREMENT
51 sum = 0.0
52 ksum = 0.0
53
54 #check for boundary condition
55 if((xx + ss < Lx) and (yy + ss < Ly)):
56 #create a weighted gaussian sum
57 gx, gy = mgrid[-ss:ss+1, -ss:ss+1]
58 ww = exp(-(gx**2/float(ss)+gy**2/float(ss)))
59 #loop over the box to determine the values
60 for ii in xrange( int(-ss), int(ss+1) ):
61 for jj in xrange( int(-ss), int(ss+1) ):
62 sum += IMG[xx + ii][yy + jj] * ww[ii + ss][jj + ss]
63 ksum += ww[ii + ss][jj + ss]
64 else:
65 break
66 qq += 1
67
68 # save the size of the box determined to reach the threshold value
69 BOX[xx][yy] = ss
70
71 # Determine the normalization for each box
72 # Norm can't be determined from the above loop because it relies on the
73 # total ksum value, if placed above the incorrect ksum value will be
74 # divided.
75 if((xx + ss < Lx) and (yy + ss < Ly)):
76 #create a weighted gaussian sum
77 gx, gy = mgrid[-ss:ss+1, -ss:ss+1]
78 ww = exp(-(gx**2/float(ss)+gy**2/float(ss)))
79 for ii in xrange( int(-ss), int(ss+1) ):
80 for jj in xrange( int(-ss), int(ss+1) ):
81 if(ksum != 0):
82 NORM[xx+ii][yy+jj] += (ww[ii + ss][jj + ss] / ksum)
83 #---------------------------------------------------------------
84
85 # Normalize the image
86 for xx in range(Lx):
87 for yy in range(Ly):
88 if(NORM[xx][yy] != 0):
89 IMG[xx][yy] /= NORM[xx][yy]
90
91 #---------------------------------------------------------------
92
93 # Output file will be smoothed with the normalized IMG
94 for xx in range(Lx):
95 for yy in range(Ly):
96 ss = BOX[xx][yy]
97 sum = 0.0
98 ksum = 0.0
99
100 if((xx + ss < Lx) and (yy + ss < Ly)):
101 #create a weighted gaussian sum for the specific BOX size
102 gx, gy = mgrid[-ss:ss+1, -ss:ss+1]
103 ww = exp(-(gx**2/float(ss)+gy**2/float(ss)))
104 for ii in xrange( int(-ss), int(ss+1) ):
105 for jj in xrange( int(-ss), int(ss+1) ):
106 sum += (IMG[xx+ii][yy+jj] * ww[ii + ss][jj + ss])
107 ksum += ww[ii + ss][jj + ss]
108 #check for divide by zero
109 if(ksum != 0):
110 OUT[xx][yy] = sum / ksum
111 else:
112 OUT[xx][yy] = 0
113
114
115
116 kernel_stop_time = time.time()
117 total_stop_time = time.time()
118 #---------------------------------------------------------------
119
120 # Save the current image.
121 imsave('{}_serial_smoothed.png'.format(os.path.splitext(file_name)[0]), OUT, cmap=cm.gray, vmin=0, vmax=1)
122
123 # Print results & save
124 print "Total Time: %f" % (total_stop_time - total_start_time)
125 print "Setup Time: %f" % (setup_stop_time - setup_start_time)
126 print "Kernel Time: %f" % (kernel_stop_time - kernel_start_time)
Something went wrong with that request. Please try again.