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.  
  81. row[0] = iir_kernel(row[0], z[pad-1], z[pad-2], z[pad-3], coefs);
  82. row[1] = iir_kernel(row[1], row[0], z[pad-1], z[pad-2], coefs);
  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) {
  92. z[i+pad-width] = row[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.  
  98. z[pad-1] = iir_kernel(z[pad-1], row[width-1], row[width-2], row[width-3],
  99. coefs);
  100. z[pad-2] = iir_kernel(z[pad-2], z[pad-1], row[width-1], row[width-2],
  101. coefs);
  102. z[pad-3] = iir_kernel(z[pad-3], z[pad-2], z[pad-1], row[width-1], 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.  
  115. row[width-1] = iir_kernel(row[width-1], z[pad-1], z[pad-2], z[pad-3],
  116. coefs);
  117. row[width-2] = iir_kernel(row[width-2], row[width-1], z[pad-1], z[pad-2],
  118. coefs);
  119. row[width-3] = iir_kernel(row[width-3], row[width-2], row[width-1],
  120. z[pad-1], coefs);
  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) {
  181. iir_v_row(width, coefs, image+(pad-i-1)*width, z+(i-1)*width,
  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.  
  187. iir_v_row(width, coefs, image, z+(pad-1)*width, z+(pad-2)*width,
  188. z+(pad-3)*width, image);
  189. iir_v_row(width, coefs, image+width, image, z+(pad-1)*width,
  190. z+(pad-2)*width, image+width);
  191. iir_v_row(width, coefs, image+2*width, image+width, image,
  192. z+(pad-1)*width, image+2*width);
  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) {
  202. memcpy(z+(i+pad-height)*width, image+i*width, width*sizeof(float));
  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.  
  209. iir_v_row(width, coefs, z+(pad-1)*width, image+(height-1)*width,
  210. image+(height-2)*width,
  211. image+(height-3)*width, z+(pad-1)*width);
  212. iir_v_row(width, coefs, z+(pad-2)*width, z+(pad-1)*width,
  213. image+(height-1)*width, image+(height-2)*width,
  214. z+(pad-2)*width);
  215. iir_v_row(width, coefs, z+(pad-3)*width, z+(pad-2)*width,
  216. z+(pad-1)*width, image+(height-1)*width, z+(pad-3)*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. }
  228. iir_v_row(width, coefs, image+(height-1)*width, z+(pad-1)*width,
  229. z+(pad-2)*width, z+(pad-3)*width, image+(height-1)*width);
  230. iir_v_row(width, coefs, image+(height-2)*width, image+(height-1)*width,
  231. z+(pad-1)*width, z+(pad-2)*width, image+(height-2)*width);
  232. iir_v_row(width, coefs, image+(height-3)*width, image+(height-2)*width,
  233. image+(height-1)*width, z+(pad-1)*width, image+(height-3)*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. }