Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Newer
Older
100644 186 lines (167 sloc) 5.099 kb
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
1 {{{id=4|
2 import Image
3 import numpy as np
4 from matplotlib import pyplot as plt
5 import scipy.stats
6 import scipy.signal
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
7 import warnings
8
9 warnings.simplefilter('ignore', np.ComplexWarning)
10 warnings.simplefilter('ignore', DeprecationWarning)
11 np.set_printoptions(precision=2)
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
12
13 _plotidx = 0
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
14 def plot_images(images):
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
15 global _plotidx
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
16 if not isinstance(images, list):
17 images = [images]
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
18 fig = plt.figure()
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
19 for idx, img in enumerate(images):
20 ax = fig.add_subplot(1, len(images), idx + 1)
21 i = ax.imshow(img, interpolation='nearest', cmap=plt.cm.gray)
22 if img.ndim == 2:
23 fig.colorbar(i)
24 fig.set_size_inches(15, 10)
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
25 _plotidx += 1
26 plt.savefig('out%d.png' % _plotidx, dpi=60, )
27 ///
28 }}}
29
30 {{{id=1|
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
31 photo = np.asarray(Image.open('/home/dagss/munich.png'))
32 photo = photo.astype(np.float32) / 256
33 small_photo = photo[543:615, 290:365]
34 plot_images([photo, small_photo])
35 ///
36 }}}
37
38 {{{id=16|
39 ksize = 4
40 x, y = np.mgrid[-ksize:ksize + 1, -ksize:ksize + 1]
41 print x
42 print y
43 ///
44 }}}
45
46 {{{id=17|
47 np.sqrt(x**2 + y**2)
48 ///
49 }}}
50
51 {{{id=18|
52 np.round(np.exp(-0.25*(x**2 + y**2)), 2)
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
53 ///
54 }}}
55
56 {{{id=2|
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
57 kernel = np.exp(-0.1*(x**2 + y**2)).astype(np.float32)
58 kernel /= np.sum(kernel)
59 plot_images([small_photo, kernel])
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
60 ///
61 }}}
62
63 {{{id=8|
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
64 def color_convolve(convolve_func, img, kernel):
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
65 smoothed_img = np.zeros((img.shape[0] + kernel.shape[0] - 1, img.shape[1] + kernel.shape[1] - 1, 3), np.float32)
66 for cidx in range(3):
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
67 smoothed_img[:, :, cidx] = convolve_func(img[:, :, cidx], kernel)
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
68 return smoothed_img
69 ///
70 }}}
71
72 {{{id=3|
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
73 time smooth_photo = color_convolve(scipy.signal.convolve2d, photo, kernel)
74 plot_images([photo, smooth_photo])
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
75 ///
76 }}}
77
78 {{{id=7|
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
79 def py_convolve2d(f, g):
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
80 # f is an image and is indexed by (v, w)
81 # g is a filter kernel and is indexed by (s, t),
82 # it needs odd dimensions
83 # h is the output image and is indexed by (x, y),
84 # it is not cropped
85 if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
86 raise ValueError("Only odd dimensions on filter supported")
87 # smid and tmid are number of pixels between the center pixel
88 # and the edge, ie for a 5x5 filter they will be 2.
89 #
90 # The output size is calculated by adding smid, tmid to each
91 # side of the dimensions of the input image.
92 vmax = f.shape[0]
93 wmax = f.shape[1]
94 smax = g.shape[0]
95 tmax = g.shape[1]
96 smid = smax // 2
97 tmid = tmax // 2
98 xmax = vmax + 2*smid
99 ymax = wmax + 2*tmid
100 # Allocate result image.
101 h = np.zeros([xmax, ymax], dtype=f.dtype)
102 # Do convolution
103 for x in range(xmax):
104 for y in range(ymax):
105 # Calculate pixel value for h at (x,y). Sum one component
106 # for each pixel (s, t) of the filter g.
107 s_from = max(smid - x, -smid)
108 s_to = min((xmax - x) - smid, smid + 1)
109 t_from = max(tmid - y, -tmid)
110 t_to = min((ymax - y) - tmid, tmid + 1)
111 value = 0
112 for s in range(s_from, s_to):
113 for t in range(t_from, t_to):
114 v = x - smid + s
115 w = y - tmid + t
116 value += g[smid - s, tmid - t] * f[v, w]
117 h[x, y] = value
118 return h
119 ///
120 }}}
121
122 {{{id=9|
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
123 time smoothed_small_photo = color_convolve(py_convolve2d, small_photo, kernel)
124 plot_images([small_photo, smoothed_small_photo])
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
125 ///
126 }}}
127
128 {{{id=10|
129 %cython
130 import numpy as np
131 cimport numpy as np
132 cimport cython
133
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
134 def cy_convolve2d(np.ndarray[np.float32_t, ndim=2] f, np.ndarray[np.float32_t, ndim=2] g):
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
135 cdef Py_ssize_t vmax, wmax, smax, tmax, smid, tmid, xmax, ymax, x, y, s, t, v, w, s_from, s_to, t_from, t_to
136 cdef np.ndarray[np.float32_t, ndim=2] h
137 cdef np.float64_t value
138
139 if g is None or f is None:
140 raise TypeError("f and g must not be None")
141 if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
142 raise ValueError("Only odd dimensions on filter supported")
143 vmax = f.shape[0]
144 wmax = f.shape[1]
145 smax = g.shape[0]
146 tmax = g.shape[1]
147 smid = smax // 2
148 tmid = tmax // 2
149 xmax = vmax + 2*smid
150 ymax = wmax + 2*tmid
151 h = np.zeros([xmax, ymax], dtype=np.float32)
152 for x in range(xmax):
153 for y in range(ymax):
154 s_from = max(smid - x, -smid)
155 s_to = min((xmax - x) - smid, smid + 1)
156 t_from = max(tmid - y, -tmid)
157 t_to = min((ymax - y) - tmid, tmid + 1)
158 value = 0
159 for s in range(s_from, s_to):
160 for t in range(t_from, t_to):
161 v = x - smid + s
162 w = y - tmid + t
163 value += g[smid - s, tmid - t] * f[v, w]
164 h[x, y] = value
165 return h
166 ///
167 }}}
168
169 {{{id=11|
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
170 time smoothed_small_photo = color_convolve(cy_convolve2d, small_photo, kernel)
171 plot_images([small_photo, smoothed_small_photo])
172 ///
173 }}}
174
175 {{{id=19|
176 time smoothed1 = color_convolve(scipy.signal.convolve2d, photo, kernel)
177 time smoothed2 = color_convolve(cy_convolve2d, photo, kernel)
178 plot_images([smoothed1, smoothed2, np.sum(smoothed1 - smoothed2, axis=2)])
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
179 ///
180 }}}
181
925f28a @dagss Improvements to cython numerics sage worksheet
dagss authored
182 {{{id=21|
097e26b @dagss First edition of Gaussian smoothing worksheet
dagss authored
183
184 ///
185 }}}
Something went wrong with that request. Please try again.