Posted to tcl by kbk at Sun Sep 20 02:48:30 GMT 2009view raw
- /*
- * Implementation of a Gaussian IIR filter after
- * Ian T. Young and Lucas J. van Vliet. "Recursive implementation of
- * the Gaussian filter." _Signal Processing_ 44 (1995), pp. 139-151.
- */
- /* Initialize the filter coefficients for the IIR filter. */
- static inline void
- iir_coefficients(float sigma, /* Radius of the filter */
- float* coefs) /* Output coefficients, normalized */
- {
- float q;
- float b0, b1, b2, b3;
- /* Correct the filter width (Eq. 11a in the paper). */
- if (sigma > 2.5f) {
- q = 0.98711f * sigma - 0.96330f;
- } else {
- q = 3.97156f - 4.14554f * sqrt(1.0f - 0.26891f * sigma);
- }
- /* Determine the filter coefficients */
- b0 = 1.57825f + q * (2.44413f + q * (1.42810f + q * 0.422205f));
- b1 = q * (2.44413f + q * (2.85619f + q * 1.26661f));
- b2 = q * q * (-1.42810f + q * -1.26661f);
- b3 = q * q * q * 0.422205f;
- /* Determine the normalization constant (Eq. 10) */
- coefs[0] = 1.0f - (b1 + b2 + b3) / b0;
- /* Scale the filter coefficients */
- coefs[1] = b1 / b0;
- coefs[2] = b2 / b0;
- coefs[3] = b3 / b0;
- }
- /* Apply the IIR filter to the input value and the last three output values */
- static inline float
- iir_kernel(float in, /* Input value */
- float out_1, /* Last output value */
- float out_2, /* Penultimate output value */
- float out_3, /* Antepenultimate output value */
- float* coefs) /* Filter coefficients */
- {
- return coefs[0] * in
- + coefs[1] * out_1 + coefs[2] * out_2 + coefs[3] * out_3;
- }
- /* Apply the IIR filter across a row of pixels, first left then right.
- * Filter "in place": replace the row with a filtered copy.
- * The 'z' array is expected to be at least 'pad' elements in length
- * and is used to handle edge effects. */
- static inline void
- iir_row (unsigned width, /* Width of the image */
- unsigned pad, /* Number of pixels of pad at the edges */
- float* coefs, /* Filter coefficients */
- float* row, /* One row of pixels to filter */
- float* z) /* Temporary array for starting and stopping */
- {
- int i;
- /* Start the filter running on the row in reverse to "prime the pump" */
- z[0] = iir_kernel(row[pad-1], 0.0f, 0.0f, 0.0f, coefs);
- z[1] = iir_kernel(row[pad-2], z[0], 0.0f, 0.0f, coefs);
- z[2] = iir_kernel(row[pad-3], z[1], z[0], 0.0f, coefs);
- for (i = 3; i < pad; ++i) {
- z[i] = iir_kernel(row[pad-i-1], z[i-1], z[i-2], z[i-3], coefs);
- }
- /* Turn around and start the row. */
- row[0] = iir_kernel(row[0], z[pad-1], z[pad-2], z[pad-3], coefs);
- row[1] = iir_kernel(row[1], row[0], z[pad-1], z[pad-2], coefs);
- row[2] = iir_kernel(row[2], row[1], row[0], z[pad-1], coefs);
- for (i = 3; i+pad < width; ++i) {
- row[i] = iir_kernel(row[i], row[i-1], row[i-2], row[i-3], coefs);
- }
- /* Stash the last 'pad' columns off in the 'z' array before
- * overwriting them */
- for (; i < width; ++i) {
- z[i+pad-width] = row[i];
- row[i] = iir_kernel(row[i], row[i-1], row[i-2], row[i-3], coefs);
- }
- /* Continue for an overrun of 'pad' columns */
- z[pad-1] = iir_kernel(z[pad-1], row[width-1], row[width-2], row[width-3],
- coefs);
- z[pad-2] = iir_kernel(z[pad-2], z[pad-1], row[width-1], row[width-2],
- coefs);
- z[pad-3] = iir_kernel(z[pad-3], z[pad-2], z[pad-1], row[width-1], coefs);
- for (i = pad-4; i >= 0; --i) {
- z[i] = iir_kernel(z[i], z[i+1], z[i+2], z[i+3], coefs);
- }
- /* Now prime the filter in the opposite direction. */
- for (i = 3; i < pad; ++i) {
- z[i] = iir_kernel(z[i], z[i-1], z[i-2], z[i-3], coefs);
- }
- /* Turn around and start filling the row */
- row[width-1] = iir_kernel(row[width-1], z[pad-1], z[pad-2], z[pad-3],
- coefs);
- row[width-2] = iir_kernel(row[width-2], row[width-1], z[pad-1], z[pad-2],
- coefs);
- row[width-3] = iir_kernel(row[width-3], row[width-2], row[width-1],
- z[pad-1], coefs);
- /* Run the filter from right to left */
- for (i = width-4; i >= 0; --i) {
- row[i] = iir_kernel(row[i], row[i+1], row[i+2], row[i+3], coefs);
- }
- }
- /*
- * Apply the IIR blur filter across the columns on one row of an image
- */
- static inline void
- iir_v_row(unsigned width, /* Width of the image */
- float* coefs, /* Filter coefficients */
- float* in_row, /* Input row */
- float* out_row_1, /* Previous output row */
- float* out_row_2, /* Penultimate output row */
- float* out_row_3, /* Antepenultimate output row */
- float* out_row) /* OUTPUT: current output row. May overwrite
- * in_row */
- {
- unsigned i;
- for (i = 0; i < width; ++i) {
- out_row[i] = iir_kernel(in_row[i], out_row_1[i], out_row_2[i],
- out_row_3[i], coefs);
- }
- }
- /*
- * Apply the IIR blur filter to an entire image
- */
- void
- iir_blur_image(unsigned height, /* Height of the image */
- unsigned width, /* Width of the image */
- float* image, /* Pixels of the image, row-major order */
- float sigma, /* Sigma of the Gaussian filter */
- unsigned pad) /* Half-length of the support */
- {
- int i;
- float coefs[4];
- float* y = malloc(pad * sizeof(float));
- float* z = malloc(pad * width * sizeof(float));
- float* zero = malloc(width * sizeof(float));
- memset(zero, 0, width * sizeof(float));
- /* Compute the coefficients */
- iir_coefficients(sigma, coefs);
- /* Start a filter running on the top fringe of the image in reverse
- * to 'prime the pump.' */
- iir_v_row(width, coefs, image + (pad-1)*width, zero, zero, zero, z);
- iir_v_row(width, coefs, image + (pad-2)*width, z, zero, zero, z+width);
- iir_v_row(width, coefs, image + (pad-3)*width, z+width, z, zero,
- z+2*width);
- for (i = 3; i < pad; ++i) {
- iir_v_row(width, coefs, image+(pad-i-1)*width, z+(i-1)*width,
- z+(i-2)*width, z+(i-3)*width, z+i*width);
- }
- /* Turn around and run the filter forward to the bottom fringe */
- iir_v_row(width, coefs, image, z+(pad-1)*width, z+(pad-2)*width,
- z+(pad-3)*width, image);
- iir_v_row(width, coefs, image+width, image, z+(pad-1)*width,
- z+(pad-2)*width, image+width);
- iir_v_row(width, coefs, image+2*width, image+width, image,
- z+(pad-1)*width, image+2*width);
- for (i = 3; i + pad < height; ++i) {
- iir_v_row(width, coefs, image+i*width, image+(i-1)*width,
- image+(i-2)*width, image+(i-3)*width, image+i*width);
- }
- /* Safely store the last 'pad' rows in the 'z' array before
- * overwriting them */
- for (; i < height; ++i) {
- memcpy(z+(i+pad-height)*width, image+i*width, width*sizeof(float));
- iir_v_row(width, coefs, image+i*width, image+(i-1)*width,
- image+(i-2)*width, image+(i-3)*width, image+i*width);
- }
- /* Continue for an overrun of 'pad' rows */
- iir_v_row(width, coefs, z+(pad-1)*width, image+(height-1)*width,
- image+(height-2)*width,
- image+(height-3)*width, z+(pad-1)*width);
- iir_v_row(width, coefs, z+(pad-2)*width, z+(pad-1)*width,
- image+(height-1)*width, image+(height-2)*width,
- z+(pad-2)*width);
- iir_v_row(width, coefs, z+(pad-3)*width, z+(pad-2)*width,
- z+(pad-1)*width, image+(height-1)*width, z+(pad-3)*width);
- for (i = pad-4; i >= 0; --i) {
- iir_v_row(width, coefs, z+i*width, z+(i+1)*width, z+(i+2)*width,
- z+(i+3)*width, z+i*width);
- }
- /* Now start the filter in the upward direction */
- for (i = 3; i < pad; ++i) {
- iir_v_row(width, coefs, z+i*width, z+(i-1)*width, z+(i-2)*width,
- z+(i-3)*width, z+i*width);
- }
- iir_v_row(width, coefs, image+(height-1)*width, z+(pad-1)*width,
- z+(pad-2)*width, z+(pad-3)*width, image+(height-1)*width);
- iir_v_row(width, coefs, image+(height-2)*width, image+(height-1)*width,
- z+(pad-1)*width, z+(pad-2)*width, image+(height-2)*width);
- iir_v_row(width, coefs, image+(height-3)*width, image+(height-2)*width,
- image+(height-1)*width, z+(pad-1)*width, image+(height-3)*width);
- /* Run the filter upward to the top of the image */
- for (i=height-4; i >= 0; --i) {
- iir_v_row(width, coefs, image+i*width, image+(i+1)*width,
- image+(i+2)*width, image+(i+3)*width, image+i*width);
- iir_row(width, pad, coefs, image+(i+3)*width, y);
- }
- iir_row(width, pad, coefs, image+2*width, y);
- iir_row(width, pad, coefs, image+width, y);
- iir_row(width, pad, coefs, image, y);
- free(zero);
- free(z);
- free(y);
- }