https://doi.org/10.5201/ipol.2016.172
Tip revision: edce3fbb7ce8b72fd316f77f13ff3f1d75091a52 authored by Software Heritage on 19 February 2016, 00:00:00 UTC
ipol: Deposit 1282 in collection ipol
ipol: Deposit 1282 in collection ipol
Tip revision: edce3fb
gaussian.h
// 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 this program. If not, see
// <http://www.opensource.org/licenses/bsd-license.html>.
//
// Copyright (C) 2012, Javier Sánchez Pérez <jsanchez@dis.ulpgc.es>
// Copyright (C) 2014, Nelson Monzón López <nmonzon@ctim.es>
// All rights reserved.
#ifndef GAUSSIAN_H
#define GAUSSIAN_H
#include <cmath>
#include <iostream>
/**
*
* Convolution with a Gaussian
*
*/
void
gaussian (float *I, //input/output image
const int xdim, //image width
const int ydim, //image height
const int zdim, // number of color channels in the image
const double sigma, //Gaussian sigma
const int bc = 1, //boundary condition
const int precision = 5 //defines the size of the window
)
{
int i, j, k;
const double den = 2 * sigma * sigma;
const int size = (int) (precision * sigma) + 1;
const int bdx = xdim + size;
const int bdy = ydim + size;
if (bc && size > xdim){
std::cerr << "GaussianSmooth: sigma too large for this bc\n" << std::endl;
throw 1;
}
// compute the coefficients of the 1D convolution kernel
double *B = new double[size];
for (int i = 0; i < size; i++)
B[i] = 1 / (sigma * sqrt (2.0 * 3.1415926)) * exp (-i * i / den);
double norm = 0;
// normalize the 1D convolution kernel
for (int i = 0; i < size; i++)
norm += B[i];
norm *= 2;
norm -= B[0];
for (int i = 0; i < size; i++)
B[i] /= norm;
double *R = new double[size + xdim + size];
double *T = new double[size + ydim + size];
//Loop for every channel
for(int index_multichannel = 0; index_multichannel < zdim; index_multichannel++){
// convolution of each line of the input image
for (k = 0; k < ydim; k++)
{
for (i = size; i < bdx; i++)
R[i] = I[(k * xdim + i - size) * zdim + index_multichannel];
switch (bc)
{
case 0: // Dirichlet boundary conditions
for (i = 0, j = bdx; i < size; i++, j++)
R[i] = R[j] = 0;
break;
case 1: // Reflecting boundary conditions
for (i = 0, j = bdx; i < size; i++, j++)
{
R[i] = I[(k * xdim + size - i ) * zdim + index_multichannel];
R[j] = I[(k * xdim + xdim - i - 1) * zdim + index_multichannel ];
}
break;
case 2: // Periodic boundary conditions
for (i = 0, j = bdx; i < size; i++, j++)
{
R[i] = I[(k * xdim + xdim - size + i) * zdim + index_multichannel];
R[j] = I[(k * xdim + i) * zdim + index_multichannel];
}
break;
}
for (i = size; i < bdx; i++)
{
double sum = B[0] * R[i];
for (int j = 1; j < size; j++)
sum += B[j] * (R[i - j] + R[i + j]);
I[(k * xdim + i - size) * zdim + index_multichannel] = sum;
}
}
// convolution of each column of the input image
for (k = 0; k < xdim; k++)
{
for (i = size; i < bdy; i++)
T[i] = I[((i - size) * xdim + k) * zdim + index_multichannel];
switch (bc)
{
case 0: // Dirichlet boundary conditions
for (i = 0, j = bdy; i < size; i++, j++)
T[i] = T[j] = 0;
break;
case 1: // Reflecting boundary conditions
for (i = 0, j = bdy; i < size; i++, j++)
{
T[i] = I[((size - i) * xdim + k) * zdim + index_multichannel];
T[j] = I[((ydim - i - 1) * xdim + k) * zdim + index_multichannel];
}
break;
case 2: // Periodic boundary conditions
for (i = 0, j = bdx; i < size; i++, j++)
{
T[i] = I[((ydim - size + i) * xdim + k) * zdim + index_multichannel];
T[j] = I[(i * xdim + k) * zdim + index_multichannel];
}
break;
}
for (i = size; i < bdy; i++)
{
double sum = B[0] * T[i];
for (j = 1; j < size; j++)
sum += B[j] * (T[i - j] + T[i + j]);
I[((i - size) * xdim + k) * zdim + index_multichannel] = sum;
}
}
}
delete[]B;
delete[]R;
delete[]T;
}
#endif