Posted to tcl by kbk at Sun Sep 20 02:48:30 GMT 2009view raw

1. /*
2. * Implementation of a Gaussian IIR filter after
3. * Ian T. Young and Lucas J. van Vliet. "Recursive implementation of
4. * the Gaussian filter." _Signal Processing_ 44 (1995), pp. 139-151.
5. */
6.
7. /* Initialize the filter coefficients for the IIR filter. */
8.
9. static inline void
10. iir_coefficients(float sigma, /* Radius of the filter */
11. float* coefs) /* Output coefficients, normalized */
12. {
13. float q;
14. float b0, b1, b2, b3;
15.
16. /* Correct the filter width (Eq. 11a in the paper). */
17.
18. if (sigma > 2.5f) {
19. q = 0.98711f * sigma - 0.96330f;
20. } else {
21. q = 3.97156f - 4.14554f * sqrt(1.0f - 0.26891f * sigma);
22. }
23.
24. /* Determine the filter coefficients */
25.
26. b0 = 1.57825f + q * (2.44413f + q * (1.42810f + q * 0.422205f));
27. b1 = q * (2.44413f + q * (2.85619f + q * 1.26661f));
28. b2 = q * q * (-1.42810f + q * -1.26661f);
29. b3 = q * q * q * 0.422205f;
30.
31. /* Determine the normalization constant (Eq. 10) */
32.
33. coefs[0] = 1.0f - (b1 + b2 + b3) / b0;
34.
35. /* Scale the filter coefficients */
36.
37. coefs[1] = b1 / b0;
38. coefs[2] = b2 / b0;
39. coefs[3] = b3 / b0;
40.
41. }
42.
43. /* Apply the IIR filter to the input value and the last three output values */
44.
45. static inline float
46. iir_kernel(float in, /* Input value */
47. float out_1, /* Last output value */
48. float out_2, /* Penultimate output value */
49. float out_3, /* Antepenultimate output value */
50. float* coefs) /* Filter coefficients */
51. {
52. return coefs[0] * in
53. + coefs[1] * out_1 + coefs[2] * out_2 + coefs[3] * out_3;
54. }
55.
56. /* Apply the IIR filter across a row of pixels, first left then right.
57. * Filter "in place": replace the row with a filtered copy.
58. * The 'z' array is expected to be at least 'pad' elements in length
59. * and is used to handle edge effects. */
60.
61. static inline void
62. iir_row (unsigned width, /* Width of the image */
63. unsigned pad, /* Number of pixels of pad at the edges */
64. float* coefs, /* Filter coefficients */
65. float* row, /* One row of pixels to filter */
66. float* z) /* Temporary array for starting and stopping */
67. {
68. int i;
69.
70. /* Start the filter running on the row in reverse to "prime the pump" */
71.
72. z[0] = iir_kernel(row[pad-1], 0.0f, 0.0f, 0.0f, coefs);
73. z[1] = iir_kernel(row[pad-2], z[0], 0.0f, 0.0f, coefs);
74. z[2] = iir_kernel(row[pad-3], z[1], z[0], 0.0f, coefs);
75. for (i = 3; i < pad; ++i) {
76. z[i] = iir_kernel(row[pad-i-1], z[i-1], z[i-2], z[i-3], coefs);
77. }
78.
79. /* Turn around and start the row. */
80.
83. row[2] = iir_kernel(row[2], row[1], row[0], z[pad-1], coefs);
84. for (i = 3; i+pad < width; ++i) {
85. row[i] = iir_kernel(row[i], row[i-1], row[i-2], row[i-3], coefs);
86. }
87.
88. /* Stash the last 'pad' columns off in the 'z' array before
89. * overwriting them */
90.
91. for (; i < width; ++i) {
93. row[i] = iir_kernel(row[i], row[i-1], row[i-2], row[i-3], coefs);
94. }
95.
96. /* Continue for an overrun of 'pad' columns */
97.
99. coefs);
101. coefs);
103. for (i = pad-4; i >= 0; --i) {
104. z[i] = iir_kernel(z[i], z[i+1], z[i+2], z[i+3], coefs);
105. }
106.
107. /* Now prime the filter in the opposite direction. */
108.
109. for (i = 3; i < pad; ++i) {
110. z[i] = iir_kernel(z[i], z[i-1], z[i-2], z[i-3], coefs);
111. }
112.
113. /* Turn around and start filling the row */
114.
116. coefs);
118. coefs);
119. row[width-3] = iir_kernel(row[width-3], row[width-2], row[width-1],
121.
122. /* Run the filter from right to left */
123.
124. for (i = width-4; i >= 0; --i) {
125. row[i] = iir_kernel(row[i], row[i+1], row[i+2], row[i+3], coefs);
126. }
127.
128. }
129.
130. /*
131. * Apply the IIR blur filter across the columns on one row of an image
132. */
133.
134. static inline void
135. iir_v_row(unsigned width, /* Width of the image */
136. float* coefs, /* Filter coefficients */
137. float* in_row, /* Input row */
138. float* out_row_1, /* Previous output row */
139. float* out_row_2, /* Penultimate output row */
140. float* out_row_3, /* Antepenultimate output row */
141. float* out_row) /* OUTPUT: current output row. May overwrite
142. * in_row */
143. {
144. unsigned i;
145. for (i = 0; i < width; ++i) {
146. out_row[i] = iir_kernel(in_row[i], out_row_1[i], out_row_2[i],
147. out_row_3[i], coefs);
148. }
149. }
150.
151. /*
152. * Apply the IIR blur filter to an entire image
153. */
154.
155. void
156. iir_blur_image(unsigned height, /* Height of the image */
157. unsigned width, /* Width of the image */
158. float* image, /* Pixels of the image, row-major order */
159. float sigma, /* Sigma of the Gaussian filter */
160. unsigned pad) /* Half-length of the support */
161. {
162. int i;
163. float coefs[4];
164. float* y = malloc(pad * sizeof(float));
165. float* z = malloc(pad * width * sizeof(float));
166. float* zero = malloc(width * sizeof(float));
167. memset(zero, 0, width * sizeof(float));
168.
169. /* Compute the coefficients */
170.
171. iir_coefficients(sigma, coefs);
172.
173. /* Start a filter running on the top fringe of the image in reverse
174. * to 'prime the pump.' */
175.
176. iir_v_row(width, coefs, image + (pad-1)*width, zero, zero, zero, z);
177. iir_v_row(width, coefs, image + (pad-2)*width, z, zero, zero, z+width);
178. iir_v_row(width, coefs, image + (pad-3)*width, z+width, z, zero,
179. z+2*width);
180. for (i = 3; i < pad; ++i) {
182. z+(i-2)*width, z+(i-3)*width, z+i*width);
183. }
184.
185. /* Turn around and run the filter forward to the bottom fringe */
186.
189. iir_v_row(width, coefs, image+width, image, z+(pad-1)*width,
191. iir_v_row(width, coefs, image+2*width, image+width, image,
193. for (i = 3; i + pad < height; ++i) {
194. iir_v_row(width, coefs, image+i*width, image+(i-1)*width,
195. image+(i-2)*width, image+(i-3)*width, image+i*width);
196. }
197.
198. /* Safely store the last 'pad' rows in the 'z' array before
199. * overwriting them */
200.
201. for (; i < height; ++i) {
203. iir_v_row(width, coefs, image+i*width, image+(i-1)*width,
204. image+(i-2)*width, image+(i-3)*width, image+i*width);
205. }
206.
207. /* Continue for an overrun of 'pad' rows */
208.
210. image+(height-2)*width,
213. image+(height-1)*width, image+(height-2)*width,
217. for (i = pad-4; i >= 0; --i) {
218. iir_v_row(width, coefs, z+i*width, z+(i+1)*width, z+(i+2)*width,
219. z+(i+3)*width, z+i*width);
220. }
221.
222. /* Now start the filter in the upward direction */
223.
224. for (i = 3; i < pad; ++i) {
225. iir_v_row(width, coefs, z+i*width, z+(i-1)*width, z+(i-2)*width,
226. z+(i-3)*width, z+i*width);
227. }
230. iir_v_row(width, coefs, image+(height-2)*width, image+(height-1)*width,
232. iir_v_row(width, coefs, image+(height-3)*width, image+(height-2)*width,
234.
235. /* Run the filter upward to the top of the image */
236.
237. for (i=height-4; i >= 0; --i) {
238. iir_v_row(width, coefs, image+i*width, image+(i+1)*width,
239. image+(i+2)*width, image+(i+3)*width, image+i*width);
240. iir_row(width, pad, coefs, image+(i+3)*width, y);
241. }
242. iir_row(width, pad, coefs, image+2*width, y);
243. iir_row(width, pad, coefs, image+width, y);
244. iir_row(width, pad, coefs, image, y);
245.
246. free(zero);
247. free(z);
248. free(y);
249. }