gaussian_conv_sii.c

Go to the documentation of this file.

#include "gaussian_conv_sii.h"
#include <assert.h>
#include <stdio.h>
#include "filter_util.h"

#ifndef M_PI

#define M_PI        3.14159265358979323846264338327950288
#endif

void sii_precomp(sii_coeffs *c, double sigma, int K)
{
    /* Elboher and Werman's optimal radii and weights. */
    const double sigma0 = 100.0 / M_PI;
    static const short radii0[SII_MAX_K - SII_MIN_K + 1][SII_MAX_K] =
        {{76, 46, 23, 0, 0},
         {82, 56, 37, 19, 0},
         {85, 61, 44, 30, 16}};
    static const float weights0[SII_MAX_K - SII_MIN_K + 1][SII_MAX_K] =
        {{0.1618f, 0.5502f, 0.9495f, 0, 0},
         {0.0976f, 0.3376f, 0.6700f, 0.9649f, 0},
         {0.0739f, 0.2534f, 0.5031f, 0.7596f, 0.9738f}};

    const int i = K - SII_MIN_K;
    double sum;
    int k;
    
    assert(c && sigma > 0 && SII_VALID_K(K));
    c->K = K;
    
    for (k = 0, sum = 0; k < K; ++k)
    {
        c->radii[k] = (long)(radii0[i][k] * (sigma / sigma0) + 0.5);
        sum += weights0[i][k] * (2 * c->radii[k] + 1);
    }
    
    for (k = 0; k < K; ++k)
        c->weights[k] = (num)(weights0[i][k] / sum);
    
    return;
}

long sii_buffer_size(sii_coeffs c, long N)
{
    long pad = c.radii[0] + 1;
    return N + 2 * pad;
}

void sii_gaussian_conv(sii_coeffs c, num *dest, num *buffer,
    const num *src, long N, long stride)
{
    num accum;
    long pad, n;
    int k;
    
    assert(dest && buffer && src && dest != buffer && src != buffer
        && N > 0 && stride != 0);
    
    pad = c.radii[0] + 1;
    buffer += pad;
    
    /* Compute cumulative sum of src over n = -pad,..., N + pad - 1. */
    for (n = -pad, accum = 0; n < N + pad; ++n)
    {
        accum += src[stride * extension(N, n)];
        buffer[n] = accum;
    }
    
    /* Compute stacked box filters. */
    for (n = 0; n < N; ++n, dest += stride)
    {
        accum = c.weights[0] * (buffer[n + c.radii[0]]
            - buffer[n - c.radii[0] - 1]);
        
        for (k = 1; k < c.K; ++k)
            accum += c.weights[k] * (buffer[n + c.radii[k]]
                - buffer[n - c.radii[k] - 1]);
        
        *dest = accum;
    }
    
    return;
}

void sii_gaussian_conv_image(sii_coeffs c, num *dest, num *buffer,
    const num *src, int width, int height, int num_channels)
{
    long num_pixels = ((long)width) * ((long)height);
    int x, y, channel;
    
    assert(dest && buffer && 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 row of the channel. */
        for (y = 0; y < height; ++y)
        {
            sii_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)
            sii_gaussian_conv(c,
                dest + x, buffer, dest + x, height, width);
        
        dest += num_pixels;
        src += num_pixels;
    }
    
    return;
}