/*
* Copyright (C) 2014 Jacques Froment <Jacques.Froment@univ-ubs.fr>
*
* --------------------------------------------------------------------
* Patent warning
* --------------------------------------------------------------------
* This file implements algorithms possibly linked to the patents
*
* # A. Buades, T. Coll and J.M. Morel, Image data processing method by
* reducing image noise, and camera integrating means for implementing
* said method, EP Patent 1,749,278 (Feb. 7, 2007).
*
* This file is made available for the exclusive aim of serving as
* scientific tool to verify the soundness and completeness of the
* algorithm description. Compilation, execution and redistribution
* of this file may violate patents rights in certain countries.
* The situation being different for every country and changing
* over time, it is your responsibility to determine which patent
* rights restrictions apply to you before you compile, use,
* modify, or redistribute this file. A patent lawyer is qualified
* to make this determination.
* If and only if they don't conflict with any patent terms, you
* can benefit from the following license terms attached to this
* file.
* --------------------------------------------------------------------
*
* 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>.
*/
/**
* @file nlmp.c
* @brief Pixelwise Non-Local Means (NLM-P and NLM-Pa), basic and SIL algorithms
* @author Jacques Froment <Jacques.Froment@univ-ubs.fr>
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#ifdef _OPENMP
#include <omp.h>
#endif
/* Usual min and max functions defined as macros */
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))
/* The following is related to the LUT (Look-Up Table) defined to get a fast
exp() by tabulating exp(-x) for x in [0, MaxExpLUT] with step 1/SampExpLUT:
*/
/* Maximum value of x when tabulating exp(-x), x>=0.
The LUT returns 0 if x is greater than this threshold
*/
#define MaxExpLUT 30.0
/* Sampling frequency for x when tabulating exp(-x), x>=0.
The sampling period (interval between two adjacent samples) is 1/SampExpLUT.
*/
#define SampExpLUT 1000
/* Minimum distance to an integer to not to be considered as an integer */
#define EpsExpLUT 1e-10
/**
* @brief Return the value of exp(-x) as tabulated by the given LUT, with
* linear interpolation between samples.
*
* @param
* LUT Precomputed LUT for exp(-x)
* x x in exp(-x), must be >= 0
*/
double explut(float *LUT, double x)
{
int ix; /* Lowest index for x in the LUT */
double xs; /* x*SampExpLUT : if it is an integer, no need for interpolation */
double xsmix; /* xs - ix as a double */
float y0; /* Value of the LUT for index ix */
float y1; /* Value of the LUT for index ix+1 */
if (x > MaxExpLUT) return(0.0);
xs=x*SampExpLUT;
ix=(int)xs;
xsmix=xs-ix;
y0=LUT[ix];
if (xsmix<EpsExpLUT) return(y0);
/* xs is not an integer, perform linear interpolation */
y1=LUT[ix+1];
return((y1-y0)*xsmix+y0);
}
/**
* @brief Main NLM-P and NLM-Pa function
*
* @param
* V Input noisy image
* sigma Noise standard deviation
* ds Patch side length is d=2*ds+1
* Ds Search window side length is D=2*Ds+1
* N1,N2 Size of the input image: N2xN1 matrix with
* N2=number of rows (dx2) and N1=number of columns (dx1)
* Nc Number of channels in the input image (1 or 3)
* a Parameter of the Gaussian Euclidean norm (NLP-Pa) or 0 for the
* plain Euclidean norm (NLM-P)
* h Filtering parameter in the Gaussian weights NLM-P and NLM-Pa
* BasicAlgo 1 to run the basic algorithm, 0 to run fast SIL algorithm
*
*/
float *nlmp(float *V, int ds, int Ds, int N1, int N2, int Nc, float a, float h,
char BasicAlgo)
{
/* Declaration of shared local variables.
* By shared variables we mean variables that are shared between threads
* (if inside a parallel section).
* Inside a parallel section, this is safe only for read-only variables.
*/
float *Vd; /* Output denoised image (return value) */
float *Vsym; /* Symmetrized noisy image */
float *ExpLUT; /* LUT (Look-Up Table) for a fast exp(): tabulate exp(-x) */
int lExpLUT; /* Length of ExpLUT (number of samples) */
double *K; /* Kernel used to compute patchs similarity */
int adr; /* Pixel address in V or in Vd */
int symadr; /* Pixel address in Vsym */
int l,lsym; /* Shift due to channel c in V/Vd and in Vsym */
int *padr; /* Patch pixels addresses */
int N; /* Number of pixels of the input image */
int NNc; /* Size of the input image (N * Nc) */
int border; /* Border size Ds+ds to be added to the input image */
int N2sym,N1sym;/* Size of the symmetrized image */
int NS; /* N for the symmetrized image */
int d; /* Patch side length */
int d2; /* Patch side size */
int D; /* Search window side length */
int D2; /* Search window side size */
int Dd; /* D*d */
int D2d; /* D2*d */
int i,j; /* Various indexes */
int x1,x2,c; /* Pixel (x1,x2) and color channel c */
double s; /* Multi-purpose variable */
/* Initialization */
N=N2*N1;NNc=Nc*N; /* Image size (one plane;all planes) */
border=Ds+ds; /* Border size to be added to perform symmetrization */
d=2*ds+1; d2=d*d; /* Patch side length and size */
D=2*Ds+1; D2=D*D; /* Search window side length and size */
Dd=D*d; D2d=D2*d;
if ((N1<=border)||(N2<=border))
{
fprintf(stderr,
"Image of size (%d,%d) too small regarding border size %d\n",
N1,N2,border);
fprintf(stderr,"Please decrease ds and/or Ds parameter.\n");
exit(EXIT_FAILURE);
}
/* Memory allocation of output image */
Vd =(float *)malloc(NNc*sizeof(float));
if (!Vd)
{
fprintf(stderr,"Not enough memory (Vd)\n");
exit(EXIT_FAILURE);
}
/* Memory allocation of symmetrized noisy image Vsym */
N1sym=N1+2*border; N2sym=N2+2*border; NS=N2sym*N1sym;
Vsym =(float *)malloc(Nc*NS*sizeof(float));
if (!Vsym)
{
fprintf(stderr,"Not enough memory (Vsym)\n");
exit(EXIT_FAILURE);
}
/* Memory allocation of the LUT (Look-Up Table) for fast exp() access */
lExpLUT=1+ceil((double)MaxExpLUT*SampExpLUT);
ExpLUT=(float *)malloc(lExpLUT*sizeof(float));
if (!ExpLUT)
{
fprintf(stderr,"Not enough memory (ExpLUT)\n");
exit(EXIT_FAILURE);
}
/* Fill the LUT */
ExpLUT[0]=1.0;
for (i=1; i<lExpLUT; i++) ExpLUT[i]=exp(-(double)i/SampExpLUT);
/* Compute the symmetrized image Vsym */
for (j=0;j<N2sym;j++)
{
if (j<border )
x2=border-j;
else
if (j>N2+border-1)
x2=2*N2+border-j-2;
else
x2=j-border;
for (i=0;i<N1sym;i++)
{
if (i<border)
x1=border-i;
else
if (i>N1+border-1)
x1=2*N1+border-i-2;
else
x1=i-border;
/* Pixel address in V and in Vsym */
adr=x2*N1+x1;
symadr=j*N1sym+i;
for (l=lsym=c=0; c<Nc; c++, l+=N, lsym+=NS)
Vsym[lsym+symadr]=V[l+adr];
}
}
/* Basic Dist2 computation is the fastest way when ds=1 */
if (ds<=1) BasicAlgo=1;
/* Allocation of the kernel and padr tabs */
if (BasicAlgo) K=(double*)malloc(d2*sizeof(double));
else K=(double*)malloc(d*sizeof(double));
if (!K)
{
fprintf(stderr,"Not enough memory (K)\n");
exit(EXIT_FAILURE);
}
padr = (int*)malloc(d2*sizeof(int));
if (!padr)
{
fprintf(stderr,"Not enough memory (padr)\n");
exit(EXIT_FAILURE);
}
/* Computation of the kernel K */
a*=2*a;
