gaussian_conv_fir.c

Go to the documentation of this file.

#include "gaussian_conv_fir.h"
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "filter_util.h"
#include "inverfc_acklam.h"

static num *make_g_trunc(num sigma, long r)
{
    num *g_trunc = NULL;
    
    if ((g_trunc = (num *)malloc(sizeof(num) * (r + 1))))
    {
        num accum = g_trunc[0] = 1.0;
        long m;
        
        for (m = 1; m <= r; ++m)
        {
            num temp = m / sigma;
            g_trunc[m] = exp(-0.5 * temp * temp);
            accum += 2.0 * g_trunc[m];
        }
        
        /* Normalize such that the filter has unit sum. */
        for (m = 0; m <= r; ++m)
            g_trunc[m] /= accum;
    }
    
    return g_trunc;
}

static void conv_sym(num *dest, const num *src, long N, long stride,
    const num *h, long r)
{
    long n;
    
    for (n = 0; n < N; ++n)
    {
        num accum = h[0] * src[stride * n];
        long m;
        
        /* Compute \sum_m h_m ( src(n - m) + src(n + m) ). */
        for (m = 1; m <= r; ++m)
            accum += h[m] * (src[stride * extension(N, n - m)]
                + src[stride * extension(N, n + m)]);
        
        dest[stride * n] = accum;
    }
    
    return;
}

int fir_precomp(fir_coeffs *c, double sigma, num tol)
{
    assert(c && sigma > 0.0 && 0.0 < tol && tol < 1.0);
    c->radius = (long)ceil(M_SQRT2 * sigma * inverfc(0.5 * tol));
    return (c->g_trunc = make_g_trunc(sigma, c->radius)) ? 1 : 0;
}

void fir_gaussian_conv(fir_coeffs c, num *dest, const num *src,
    long N, long stride)
{
    assert(c.g_trunc && dest && src && dest != src && N > 0 && stride != 0);
    conv_sym(dest, src, N, stride, c.g_trunc, c.radius);
    return;
}

void fir_gaussian_conv_image(fir_coeffs c, num *dest, num *buffer,
    const num *src, int width, int height, int num_channels)
{
    const long num_pixels = ((long)width) * ((long)height);
    int x, y, channel;
    
    assert(c.g_trunc && dest && buffer
        && src && dest != src && num_pixels > 0);
    
    /* Loop over the image channels. */
    for (channel = 0; channel < num_channels; ++channel)
    {
        num *dest_y = dest;
        const num *src_y = src;
        
        /* Filter each column of the channel. */
        for (x = 0; x < width; ++x)
            conv_sym(dest + x, src + x, height, width, c.g_trunc, c.radius);
        
        /* Filter each row of the channel. */
        for (y = 0; y < height; ++y)
        {
            conv_sym(buffer, dest_y, width, 1, c.g_trunc, c.radius);
            memcpy(dest_y, buffer, sizeof(num) * width);
            dest_y += width;
            src_y += width;
        }
        
        dest += num_pixels;
        src += num_pixels;
    }
    
    return;
}

void fir_free(fir_coeffs *c)
{
    if (c && c->g_trunc)
    {
        free(c->g_trunc);
        c->g_trunc = NULL;
    }
    
    return;
}