gaussian_conv_am.c

Go to the documentation of this file.

#include <assert.h>
#include <math.h>
#include "filter_util.h"
#include "gaussian_conv_am.h"

static num am_left_boundary(const num *data, long N, long stride,
    num nu, long num_terms)
{
    num h = 1, accum = data[0];
    long m;
    
    for (m = 1; m < num_terms; ++m)
    {
        h *= nu;
        accum += h * data[stride * extension(N, -m)];
    }
    
    return accum;
}

void am_gaussian_conv(num *dest, const num *src, long N, long stride,
    double sigma, int K, num tol, int use_adjusted_q)
{
    const long stride_N = stride * N;
    double q, lambda, dnu;
    num nu, scale;
    long i, num_terms;
    int pass;
    
    assert(dest && src && N > 0 && stride != 0 && sigma > 0
        && K > 0 && tol > 0);
    
    if (use_adjusted_q)  /* Use a regression on q for improved accuracy. */
        q = sigma * (1.0 + (0.3165 * K + 0.5695)
                / ((K + 0.7818) * (K + 0.7818)));
    else                /* Use q = sigma as in the original A-M method. */
        q = sigma;
    
    /* Precompute the filter coefficient nu. */
    lambda = (q * q) / (2.0 * K);
    dnu = (1.0 + 2.0*lambda - sqrt(1.0 + 4.0*lambda))/(2.0*lambda);
    nu = (num)dnu;
    /* For handling the left boundary, determine the number of terms needed to
       approximate the sum with accuracy tol. */
    num_terms = (long)ceil(log((1.0 - dnu)*tol) / log(dnu));
    /* Precompute the constant scale factor. */
    scale = (num)(pow(dnu / lambda, K));
    
    /* Copy src to dest and multiply by the constant scale factor. */
    for (i = 0; i < stride_N; i += stride)
        dest[i] = src[i] * scale;
    
    /* Perform K passes of filtering. */
    for (pass = 0; pass < K; ++pass)
    {
        /* Initialize the recursive filter on the left boundary. */
        dest[0] = am_left_boundary(dest, N, stride, nu, num_terms);
        
        /* This loop applies the causal filter, implementing the pseudocode
           
           For n = 1, ..., N - 1
               dest(n) = dest(n) + nu dest(n - 1)
           
           Variable i = stride * n is the offset to the nth sample.  */
        for (i = stride; i < stride_N; i += stride)
            dest[i] += nu * dest[i - stride];
        
        /* Handle the right boundary. */
        i -= stride;
        dest[i] /= (1 - nu);
        
        /* Similarly, this loop applies the anticausal filter as
           
           For n = N - 1, ..., 1
               dest(n - 1) = dest(n - 1) + nu dest(n) */
        for (; i > 0; i -= stride)
            dest[i - stride] += nu * dest[i];
    }
    
    return;
}

void am_gaussian_conv_image(num *dest, const num *src,
    int width, int height, int num_channels,
    num sigma, int K, num tol, int use_adjusted_q)
{
    const long num_pixels = ((long)width) * ((long)height);
    int x, y, channel;
    
    assert(dest && src && num_pixels > 0 && sigma > 0
        && K > 0 && tol > 0);
    
    /* Loop over the image channels. */
    for (channel = 0; channel < num_channels; ++channel)
    {
        num *dest_y = dest;
        const num *src_y = src;
        
        /* Filter each row of the channel. */
        for (y = 0; y < height; ++y)
        {
            am_gaussian_conv(dest_y, src_y, width, 1,
                sigma, K, tol, use_adjusted_q);
            dest_y += width;
            src_y += width;
        }
        
        /* Filter each column of the channel. */
        for (x = 0; x < width; ++x)
            am_gaussian_conv(dest + x, dest + x, height, width,
                sigma, K, tol, use_adjusted_q);
        
        dest += num_pixels;
        src += num_pixels;
    }
    
    return;
}