https://doi.org/10.5201/ipol.2012.g-tvdc
Tip revision: 13cac45f201f2e0e4a336450abb835be0ddd9359 authored by Software Heritage on 12 June 2012, 00:00:00 UTC
ipol: Deposit 1194 in collection ipol
ipol: Deposit 1194 in collection ipol
Tip revision: 13cac45
kernels.c
/**
* @file kernels.c
* @brief Convolution kernels
* @author Pascal Getreuer <getreuer@gmail.com>
*
*
* Copyright (c) 2011, Pascal Getreuer
* All rights reserved.
*
* This program is free software: you can use, modify and/or
* redistribute it under the terms of the simplified BSD License. You
* should have received a copy of this license along with this program.
* If not, see <http://www.opensource.org/licenses/bsd-license.html>.
*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "kernels.h"
#ifndef M_SQRT1_2
/** @brief The constant sqrt(1/2) */
#define M_SQRT1_2 0.70710678118654752440084436210484904
#endif
#ifndef M_PI
/** @brief The constant 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;
}
/**
* @brief Read a kernel from the command line option stirng
* @param Kernel pointer to kernel image
* @param String command line optionstring
* @return 1 on success, 0 on failure
*
* The syntax of string can be "disk:<radius>" or "gaussian:<sigma>".
* Otherwise, the string is assumed to be a file name.
*/
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;
}