Gaussian IIR filter

/* 
 * 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);
}