gaussian_short_conv.c

Go to the documentation of this file.

#include "gaussian_short_conv.h"
#include <assert.h>
#include <math.h>

#ifndef M_PI_4

#define M_PI_4      0.78539816339744830961566084581987572
#endif

/* Define coefficients for computing short DCTs */
#ifndef M_1_SQRT2

#define M_1_SQRT2   0.70710678118654752440084436210484904
#endif
#ifndef M_SQRT3

#define M_SQRT3     1.73205080756887729352744634150587237
#endif
#ifndef M_COSPI_8

#define M_COSPI_8   0.92387953251128675612818318939678829
#endif
#ifndef M_COSPI3_8

#define M_COSPI3_8  0.38268343236508977172845998403039887
#endif


void gaussian_short_conv(num *dest, const num *src,
    long N, long stride, num sigma)
{
    double alpha = (2.0 * M_PI_4 * M_PI_4) * (sigma * sigma);
    
    assert(dest && src && 0 < N && N <= 4 && stride != 0 && sigma > 0);
    
    switch (N)
    {
    case 1:
        dest[0] = src[0];
        break;
    case 2:
        {
            double F0 = 0.5 * (src[0] + src[stride]);
            double F1 = 0.5 * (src[0] - src[stride]);
            F1 *= exp(-alpha);
            dest[0] = (num)(F0 + F1);
            dest[stride] = (num)(F0 - F1);
        }
        break;
    case 3:
        {
            double F0 = (src[0] + src[stride] + src[stride * 2]) / 3.0;
            double F1 = ((0.5*M_SQRT3) * (src[0] - src[stride * 2])) / 3.0;
            double F2 = (0.5*(src[0] + src[stride * 2]) - src[stride]) / 3.0;
            F1 *= exp(-alpha);
            F2 *= exp(-4.0 * alpha);
            dest[0] = (num)(F0 + M_SQRT3 * F1 + F2);
            dest[stride] = (num)(F0 - 2.0 * F2);
            dest[stride * 2] = (num)(F0 - M_SQRT3 * F1 + F2);
        }
        break;
    case 4:
        {
            double F0 = (src[0] + src[stride]
                + src[stride * 2] + src[stride * 3]) / 4.0;
            double F1 = (M_COSPI_8 * (src[0] - src[stride * 3])
                + M_COSPI3_8 * (src[stride] - src[stride * 2])) / 4.0;
            double F2 = (M_1_SQRT2 * (src[0] - src[stride]
                + src[stride * 2] - src[stride * 3])) / 4.0;
            double F3 = (M_COSPI3_8 * (src[0] - src[stride * 3])
                - M_COSPI_8 * (src[stride] - src[stride * 2])) / 4.0;
            F1 *= exp(-alpha);
            F2 *= exp(-4.0 * alpha);
            F3 *= exp(-9.0 * alpha);
            dest[0] = (num)(F0 + 2.0 *
                (M_COSPI_8 * F1 + M_1_SQRT2 * F2 + M_COSPI3_8 * F3));
            dest[stride] = (num)(F0 + 2.0 *
                (M_COSPI3_8 * F1 - M_1_SQRT2 * F2 - M_COSPI_8 * F3));
            dest[stride * 2] = (num)(F0 + 2.0 *
                (-M_COSPI3_8 * F1 - M_1_SQRT2 * F2 + M_COSPI_8 * F3));
            dest[stride * 3] = (num)(F0 + 2.0 *
                (-M_COSPI_8 * F1 + M_1_SQRT2 * F2 - M_COSPI3_8 * F3));
        }
        break;
    }
    
    return;
}