kernels.c

Go to the documentation of this file.

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "kernels.h"

#ifndef M_SQRT1_2

#define M_SQRT1_2   0.70710678118654752440084436210484904
#endif

#ifndef M_PI

#define M_PI        3.14159265358979323846264338327950288
#endif


static int KernelRescale(image *Kernel)
{
    const int NumEl = Kernel->Width * Kernel->Height;
    num *Data = Kernel->Data;
    num Sum = 0;
    int n;

    for(n = 0; n < NumEl; n++)
        Sum += Data[n];

    if(Sum == 0)
    {
        fprintf(stderr, "Kernel must have nonzero sum.\n");
        return 0;
    }

    for(n = 0; n < NumEl; n++)
        Data[n] /= Sum;

    return 1;
}


static int GaussianKernel(image *Kernel, num Sigma)
{
    const int R = (int)ceil(4*Sigma);
    const int W = 2*R + 1;
    const double ExpDenom = 2.0*(double)Sigma*(double)Sigma;
    num *Data;
    double Sum;
    int x, y;
    

    if(!Kernel || Sigma < 0.0 || !(AllocImageObj(Kernel, W, W, 1)))
        return 0;
    
    Data = Kernel->Data;
    
    if(Sigma == 0.0)
        Data[0] = 1;
    else
    {
        for(y = -R, Sum = 0; y <= R; y++)
            for(x = -R; x <= R; x++)
            {
                Data[(R + x) + W*(R + y)] = (num)exp(-(x*x + y*y)/ExpDenom);
                Sum += Data[(R + x) + W*(R + y)];
            }

        for(y = -R; y <= R; y++)
            for(x = -R; x <= R; x++)
                Data[(R + x) + W*(R + y)] =
                    (num)(Data[(R + x) + W*(R + y)]/Sum);
    }
    
    return 1;
}


static int DiskKernel(image *Kernel, num Radius)
{
    const int R = (int)ceil(Radius - 0.5);
    const int W = 2*R + 1;
    const int Res = 8;
    const double RadiusSqr = (double)Radius*(double)Radius;
    const double RadiusInnerSqr =
        ((double)Radius - M_SQRT1_2)*((double)Radius - M_SQRT1_2);
    const double RadiusOuterSqr =
        ((double)Radius + M_SQRT1_2)*((double)Radius + M_SQRT1_2);
    num *Data;
    double Sum, xl, yl, Start = -0.5 + 0.5/Res, Step = 1.0/Res;
    int c, x, y, m, n;
    

    if(!Kernel || Radius < 0.0 || !(AllocImageObj(Kernel, W, W, 1)))
        return 0;
    
    Data = Kernel->Data;
    
    if(Radius <= 0.5)
        Data[0] = 1;
    else
    {
        for(y = -R, Sum = 0; y <= R; y++)
            for(x = -R; x <= R; x++)
            {
                if(x*x + y*y <= RadiusInnerSqr)
                    c = Res*Res;
                else if(x*x + y*y > RadiusOuterSqr)
                    c = 0;
                else
                    for(n = 0, yl = y + Start, c = 0; n < Res; n++, yl += Step)
                        for(m = 0, xl = x + Start; m < Res; m++, xl += Step)
                            if(xl*xl + yl*yl <= RadiusSqr)
                                c++;
                            
                Data[(R + x) + W*(R + y)] = (num)c;
                Sum += c;
            }

        for(y = -R; y <= R; y++)
            for(x = -R; x <= R; x++)
                Data[(R + x) + W*(R + y)] =
                    (num)(Data[(R + x) + W*(R + y)]/Sum);
    }

    return 1;
}

    
static int MakeNamedKernel(image *Kernel, const char *String)
{
    const char *ColonPtr;
    num KernelParam;
    int Length;
    char KernelName[32];
    
    if(!Kernel || !(ColonPtr = strchr(String, ':'))
        || (Length = (int)(ColonPtr - String)) > 9)
        return 0;
                
    strncpy(KernelName, String, Length);
    KernelName[Length] = '\0';
    KernelParam = (num)atof(ColonPtr + 1);
    
    if(!strcmp(KernelName, "Gaussian") || !strcmp(KernelName, "gaussian"))
        return GaussianKernel(Kernel, KernelParam);
    else if(!strcmp(KernelName, "Disk") || !strcmp(KernelName, "disk"))
        return DiskKernel(Kernel, KernelParam);
    else
        return 0;
}


int ReadKernel(image *Kernel, const char *String)
{
    if(Kernel->Data)
    {
        FreeImageObj(*Kernel);
        *Kernel = NullImage;
    }
    
    if(MakeNamedKernel(Kernel, String)
        || ReadMatrixFromFile(Kernel, String, KernelRescale))
        return 1;
    else
        return 0;
}