gaussian_conv_ebox.c

Go to the documentation of this file.

#include "gaussian_conv_ebox.h"
#include <assert.h>
#include <math.h>
#include <string.h>
#include "filter_util.h"

void ebox_precomp(ebox_coeffs *c, double sigma, int K)
{
    double alpha;
    
    assert(c && sigma > 0 && K > 0);
    
    /* Set the parameters according to Gwosdek et al. */
    c->r = (long)(0.5 * sqrt((12.0 * sigma * sigma) / K + 1.0) - 0.5);
    alpha = (2 * c->r + 1) * (c->r * (c->r + 1) - 3.0 * sigma * sigma / K)
        / (6.0 * (sigma * sigma / K - (c->r + 1) * (c->r + 1)));
    c->c_1 = alpha / (2.0 * (alpha + c->r) + 1);
    c->c_2 = (1.0 - alpha) / (2.0 * (alpha + c->r) + 1);
    c->K = K;
    return;
}

static void ebox_filter(num *dest, long dest_stride,
    const num *src, long src_stride, long N, long r, num c_1, num c_2)
{
    long n;
    num accum = 0;
    
    assert(dest && src && dest != src && N > 0 && r >= 0);
    
    for (n = -r; n <= r; ++n)
        accum += src[src_stride * extension(N, n)];
    
    dest[0] = accum = c_1 * (src[src_stride * extension(N, r + 1)]
        + src[src_stride * extension(N, -r - 1)])
        + (c_1 + c_2) * accum;
    
    for (n = 1; n < N; ++n)
    {
        accum += c_1 * (src[src_stride * extension(N, n + r + 1)]
            - src[src_stride * extension(N, n - r - 2)])
            + c_2 * (src[src_stride * extension(N, n + r)]
            - src[src_stride * extension(N, n - r - 1)]);
        dest[dest_stride * n] = accum;
    }
    
    return;
}

void ebox_gaussian_conv(ebox_coeffs c, num *dest_data, num *buffer_data,
    const num *src, long N, long stride)
{
    struct
    {
        num *data;
        long stride;
    } dest, buffer, cur, next;
    int step;
    
    assert(dest_data && buffer_data && src
        && dest_data != buffer_data && N > 0);
    
    dest.data = dest_data;
    dest.stride = stride;
    buffer.data = buffer_data;
    buffer.stride = (buffer_data == src) ? stride : 1;
    
    next = (buffer_data == src || (dest_data != src && c.K % 2 == 1))
        ? dest : buffer;
    /* Perform the first filtering pass. */
    ebox_filter(next.data, next.stride, src, stride,
        N, c.r, c.c_1, c.c_2);
    
    /* Perform (K - 1) filtering passes, alternating the roles of the dest
       and buffer arrays. */
    for (step = 1; step < c.K; ++step)
    {
        cur = next;
        next = (cur.data == buffer_data) ? dest : buffer;
        ebox_filter(next.data, next.stride, cur.data, cur.stride,
            N, c.r, c.c_1, c.c_2);
    }
    
    /* If necessary, copy the result to the destination array. */
    if (next.data != dest_data)
    {
        if (stride == 1)
            memcpy(dest_data, buffer_data, sizeof(num) * N);
        else
        {
            long n, i;
            
            for (n = i = 0; n < N; ++n, i += stride)
                dest_data[i] = buffer_data[n];
        }
    }
    
    return;
}

void ebox_gaussian_conv_image(ebox_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(dest && buffer && src && dest != buffer && 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 row of the channel. */
        for (y = 0; y < height; ++y)
        {
            ebox_gaussian_conv(c, dest_y, buffer, src_y, width, 1);
            dest_y += width;
            src_y += width;
        }
        
        /* Filter each column of the channel. */
        for (x = 0; x < width; ++x)
            ebox_gaussian_conv(c, dest + x, buffer, dest + x, height, width);
        
        dest += num_pixels;
        src += num_pixels;
    }
    
    return;
}