gaussian_conv_dct.c

Go to the documentation of this file.

#include "gaussian_conv_dct.h"
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#ifndef M_PI

#define M_PI        3.14159265358979323846264338327950288
#endif

int dct_precomp(dct_coeffs *c, num *dest, const num *src,
    long N, long stride, double sigma)
{
    const fftw_r2r_kind dct_2 = FFTW_REDFT10, dct_3 = FFTW_REDFT01;
    double temp;
    int length = N;
    
    assert(c && dest && src && N > 0 && stride != 0 && sigma > 0);
    c->forward_plan = c->inverse_plan = NULL;
    
    if (!(c->forward_plan = FFT(plan_many_r2r)(1, &length, 1, (num *)src,
        NULL, stride, 0, dest, NULL, stride, 0, &dct_2,
        FFTW_ESTIMATE | ((src != dest) ? FFTW_PRESERVE_INPUT : 0)))
        || !(c->inverse_plan = FFT(plan_many_r2r)(1, &length, 1, dest,
        NULL, stride, 0, dest, NULL, stride, 0, &dct_3, FFTW_ESTIMATE)))
    {
        dct_free(c);
        return 0;
    }
    
    c->dest = dest;
    c->src = src;
    c->conv_type = DCT_GAUSSIAN_1D;
    temp = (sigma * M_PI) / N;
    c->dims.one.alpha = (num)(temp * temp / 2);
    c->dims.one.N = N;
    c->dims.one.stride = stride;
    return 1;
}

int dct_precomp_image(dct_coeffs *c, num *dest, const num *src,
    int width, int height, int num_channels, double sigma)
{
    const fftw_r2r_kind dct_2[] = {FFTW_REDFT10, FFTW_REDFT10};
    const fftw_r2r_kind dct_3[] = {FFTW_REDFT01, FFTW_REDFT01};
    const int dist = width * height;
    double temp;
    int size[2];
    
    assert(c && dest && src && width > 0 && height > 0 && sigma > 0);
    size[1] = width;
    size[0] = height;
    c->forward_plan = c->inverse_plan = NULL;
    
    if (!(c->forward_plan = FFT(plan_many_r2r)(2, size, num_channels,
        (num *)src, NULL, 1, dist, dest, NULL, 1, dist, dct_2,
        FFTW_ESTIMATE | ((src != dest) ? FFTW_PRESERVE_INPUT : 0)))
        || !(c->inverse_plan = FFT(plan_many_r2r)(2, size, num_channels,
        dest, NULL, 1, dist, dest, NULL, 1, dist, dct_3, FFTW_ESTIMATE)))
    {
        dct_free(c);
        return 0;
    }
    
    c->dest = dest;
    c->src = src;
    c->conv_type = DCT_GAUSSIAN_IMAGE;
    temp = (sigma * M_PI) / width;
    c->dims.image.alpha_x = (num)(temp * temp / 2);
    temp = (sigma * M_PI) / height;
    c->dims.image.alpha_y = (num)(temp * temp / 2);
    c->dims.image.width = width;
    c->dims.image.height = height;
    c->dims.image.num_channels = num_channels;
    return 1;
}

void dct_gaussian_conv(dct_coeffs c)
{
    assert(c.forward_plan && c.inverse_plan);
    
    /* Perform the forward DCT-II transform. */
    FFT(execute)(c.forward_plan);
    
    /* Perform spectral domain multiplication with the DCT-I transform of the
       Gaussian, which is exp(-alpha n^2) for the nth mode. */
    if (c.conv_type == DCT_GAUSSIAN_1D)
    {   /* Multiplication in one dimension. */
        num denom = 2 * c.dims.one.N;
        long n;
        
        for (n = 0; n < c.dims.one.N; ++n)
        {
            *c.dest *= ((num)exp(-c.dims.one.alpha * n * n)) / denom;
            c.dest += c.dims.one.stride;
        }
    }
    else
    {   /* Multiplication in two dimensions. */
        num denom = 4 * ((num)c.dims.image.width)
            * ((num)c.dims.image.height);
        int x, y, channel;
        
        for (channel = 0; channel < c.dims.image.num_channels; ++channel)
            for (y = 0; y < c.dims.image.height; ++y)
            {
                for (x = 0; x < c.dims.image.width; ++x)
                    c.dest[x] *= ((num)exp(
                        -c.dims.image.alpha_x * x * x
                        -c.dims.image.alpha_y * y * y)) / denom;
                
                c.dest += c.dims.image.width;
            }
    }
    
    /* Perform the inverse DCT-II transform. */
    FFT(execute)(c.inverse_plan);
    return;
}

void dct_free(dct_coeffs *c)
{
    assert(c);
    
    if (c->inverse_plan)
        FFT(destroy_plan)(c->inverse_plan);
    if (c->forward_plan)
        FFT(destroy_plan)(c->forward_plan);
    
    FFT(cleanup)();
    c->forward_plan = c->inverse_plan = NULL;
    return;
}