mirror of
https://github.com/RetroDECK/ES-DE.git
synced 2024-11-24 07:05:39 +00:00
243 lines
12 KiB
C
243 lines
12 KiB
C
|
/*
|
||
|
|
#
|
||
|
|
# File : nlmeans.h
|
||
|
|
# ( C++ header file - CImg plug-in )
|
||
|
|
#
|
||
|
|
# Description : CImg plugin that implements the non-local mean filter.
|
||
|
|
# This file is a part of the CImg Library project.
|
||
|
|
# ( http://cimg.eu )
|
||
|
|
#
|
||
|
|
# [1] Buades, A.; Coll, B.; Morel, J.-M.: A non-local algorithm for image denoising
|
||
|
|
# IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2005. CVPR 2005.
|
||
|
|
# Volume 2, 20-25 June 2005 Page(s):60 - 65 vol. 2
|
||
|
|
#
|
||
|
|
# [2] Buades, A. Coll, B. and Morel, J.: A review of image denoising algorithms, with a new one.
|
||
|
|
# Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal 4 (2004) 490-530
|
||
|
|
#
|
||
|
|
# [3] Gasser, T. Sroka,L. Jennen Steinmetz,C. Residual variance and residual pattern nonlinear regression.
|
||
|
|
# Biometrika 73 (1986) 625-659
|
||
|
|
#
|
||
|
|
# Copyright : Jerome Boulanger
|
||
|
|
# ( http://www.irisa.fr/vista/Equipe/People/Jerome.Boulanger.html )
|
||
|
|
#
|
||
|
|
# License : CeCILL v2.0
|
||
|
|
# ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )
|
||
|
|
#
|
||
|
|
# This software is governed by the CeCILL license under French law and
|
||
|
|
# abiding by the rules of distribution of free software. You can use,
|
||
|
|
# modify and/ or redistribute the software under the terms of the CeCILL
|
||
|
|
# license as circulated by CEA, CNRS and INRIA at the following URL
|
||
|
|
# "http://www.cecill.info".
|
||
|
|
#
|
||
|
|
# As a counterpart to the access to the source code and rights to copy,
|
||
|
|
# modify and redistribute granted by the license, users are provided only
|
||
|
|
# with a limited warranty and the software's author, the holder of the
|
||
|
|
# economic rights, and the successive licensors have only limited
|
||
|
|
# liability.
|
||
|
|
#
|
||
|
|
# In this respect, the user's attention is drawn to the risks associated
|
||
|
|
# with loading, using, modifying and/or developing or reproducing the
|
||
|
|
# software by the user in light of its specific status of free software,
|
||
|
|
# that may mean that it is complicated to manipulate, and that also
|
||
|
|
# therefore means that it is reserved for developers and experienced
|
||
|
|
# professionals having in-depth computer knowledge. Users are therefore
|
||
|
|
# encouraged to load and test the software's suitability as regards their
|
||
|
|
# requirements in conditions enabling the security of their systems and/or
|
||
|
|
# data to be ensured and, more generally, to use and operate it in the
|
||
|
|
# same conditions as regards security.
|
||
|
|
#
|
||
|
|
# The fact that you are presently reading this means that you have had
|
||
|
|
# knowledge of the CeCILL license and that you accept its terms.
|
||
|
|
#
|
||
|
|
*/
|
||
|
|
|
||
|
|
#ifndef cimg_plugin_nlmeans
|
||
|
|
#define cimg_plugin_nlmeans
|
||
|
|
|
||
|
|
//! NL-Means denoising algorithm.
|
||
|
|
/**
|
||
|
|
This is the in-place version of get_nlmean().
|
||
|
|
**/
|
||
|
|
CImg<T>& nlmeans(int patch_size=1, double lambda=-1, double alpha=3, double sigma=-1, int sampling=1){
|
||
|
|
if (!is_empty()){
|
||
|
|
if (sigma<0) sigma = std::sqrt(variance_noise()); // noise variance estimation
|
||
|
|
const double np = (2*patch_size + 1)*(2*patch_size + 1)*spectrum()/(double)sampling;
|
||
|
|
if (lambda<0) {// Bandwidth estimation
|
||
|
|
if (np<100)
|
||
|
|
lambda = ((((((1.1785e-12*np - 5.1827e-10)*np + 9.5946e-08)*np -
|
||
|
|
9.7798e-06)*np + 6.0756e-04)*np - 0.0248)*np + 1.9203)*np + 7.9599;
|
||
|
|
else
|
||
|
|
lambda = (-7.2611e-04*np + 1.3213)*np + 15.2726;
|
||
|
|
}
|
||
|
|
#if cimg_debug>=1
|
||
|
|
std::fprintf(stderr,"Size of the patch : %dx%d \n",
|
||
|
|
2*patch_size + 1,2*patch_size + 1);
|
||
|
|
std::fprintf(stderr,"Size of window where similar patch are looked for : %dx%d \n",
|
||
|
|
(int)(alpha*(2*patch_size + 1)),(int)(alpha*(2*patch_size + 1)));
|
||
|
|
std::fprintf(stderr,"Bandwidth of the kernel : %fx%f^2 \n",
|
||
|
|
lambda,sigma);
|
||
|
|
std::fprintf(stderr,"Noise standard deviation estimated to : %f \n",
|
||
|
|
sigma);
|
||
|
|
#endif
|
||
|
|
|
||
|
|
CImg<T> dest(width(),height(),depth(),spectrum(),0);
|
||
|
|
double *uhat = new double[spectrum()];
|
||
|
|
const double h2 = -.5/(lambda*sigma*sigma); // [Kervrann] notations
|
||
|
|
if (depth()!=1){ // 3D case
|
||
|
|
const CImg<> P = (*this).get_blur(1); // inspired from Mahmoudi&Sapiro SPletter dec 05
|
||
|
|
const int n_simu = 64;
|
||
|
|
CImg<> tmp(n_simu,n_simu,n_simu);
|
||
|
|
const double sig = std::sqrt(tmp.fill(0.f).noise(sigma).blur(1).pow(2.).sum()/(n_simu*n_simu*n_simu));
|
||
|
|
const int
|
||
|
|
patch_size_z = 0,
|
||
|
|
pxi = (int)(alpha*patch_size),
|
||
|
|
pyi = (int)(alpha*patch_size),
|
||
|
|
pzi = 2; //Define the size of the neighborhood in z
|
||
|
|
for (int zi = 0; zi<depth(); ++zi) {
|
||
|
|
#if cimg_debug>=1
|
||
|
|
std::fprintf(stderr,"\rProcessing : %3d %%",(int)((float)zi/(float)depth()*100.));fflush(stdout);
|
||
|
|
#endif
|
||
|
|
for (int yi = 0; yi<height(); ++yi)
|
||
|
|
for (int xi = 0; xi<width(); ++xi) {
|
||
|
|
cimg_forC(*this,v) uhat[v] = 0;
|
||
|
|
float sw = 0, wmax = -1;
|
||
|
|
for (int zj = std::max(0,zi - pzi); zj<std::min(depth(),zi + pzi + 1); ++zj)
|
||
|
|
for (int yj = std::max(0,yi - pyi); yj<std::min(height(),yi + pyi + 1); ++yj)
|
||
|
|
for (int xj = std::max(0,xi - pxi); xj<std::min(width(),xi + pxi + 1); ++xj)
|
||
|
|
if (cimg::abs(P(xi,yi,zi) - P(xj,yj,zj))/sig<3) {
|
||
|
|
double d = 0;
|
||
|
|
int n = 0;
|
||
|
|
if (xi!=xj && yi!=yj && zi!=zj){
|
||
|
|
for (int kz = -patch_size_z; kz<patch_size_z + 1; kz+=sampling) {
|
||
|
|
int
|
||
|
|
zj_ = zj + kz,
|
||
|
|
zi_ = zi + kz;
|
||
|
|
if (zj_>=0 && zj_<depth() && zi_>=0 && zi_<depth())
|
||
|
|
for (int ky = -patch_size; ky<=patch_size; ky+=sampling) {
|
||
|
|
int
|
||
|
|
yj_ = yj + ky,
|
||
|
|
yi_ = yi + ky;
|
||
|
|
if (yj_>=0 && yj_<height() && yi_>=0 && yi_<height())
|
||
|
|
for (int kx = -patch_size; kx<=patch_size; kx+=sampling) {
|
||
|
|
int
|
||
|
|
xj_ = xj + kx,
|
||
|
|
xi_ = xi + kx;
|
||
|
|
if (xj_>=0 && xj_<width() && xi_>=0 && xi_<width())
|
||
|
|
cimg_forC(*this,v) {
|
||
|
|
double d1 = (*this)(xj_,yj_,zj_,v) - (*this)(xi_,yi_,zi_,v);
|
||
|
|
d+=d1*d1;
|
||
|
|
++n;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
float w = (float)std::exp(d*h2);
|
||
|
|
wmax = w>wmax?w:wmax;
|
||
|
|
cimg_forC(*this,v) uhat[v]+=w*(*this)(xj,yj,zj,v);
|
||
|
|
sw+=w;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
// add the central pixel
|
||
|
|
cimg_forC(*this,v) uhat[v]+=wmax*(*this)(xi,yi,zi,v);
|
||
|
|
sw+=wmax;
|
||
|
|
if (sw) cimg_forC(*this,v) dest(xi,yi,zi,v) = (T)(uhat[v]/=sw);
|
||
|
|
else cimg_forC(*this,v) dest(xi,yi,zi,v) = (*this)(xi,yi,zi,v);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else { // 2D case
|
||
|
|
const CImg<> P = (*this).get_blur(1); // inspired from Mahmoudi&Sapiro SPletter dec 05
|
||
|
|
const int n_simu = 512;
|
||
|
|
CImg<> tmp(n_simu,n_simu);
|
||
|
|
const double sig = std::sqrt(tmp.fill(0.f).noise(sigma).blur(1).pow(2.).sum()/(n_simu*n_simu));
|
||
|
|
const int
|
||
|
|
pxi = (int)(alpha*patch_size),
|
||
|
|
pyi = (int)(alpha*patch_size); //Define the size of the neighborhood
|
||
|
|
for (int yi = 0; yi<height(); ++yi) {
|
||
|
|
#if cimg_debug>=1
|
||
|
|
std::fprintf(stderr,"\rProcessing : %3d %%",(int)((float)yi/(float)height()*100.));fflush(stdout);
|
||
|
|
#endif
|
||
|
|
for (int xi = 0; xi<width(); ++xi) {
|
||
|
|
cimg_forC(*this,v) uhat[v] = 0;
|
||
|
|
float sw = 0, wmax = -1;
|
||
|
|
for (int yj = std::max(0,yi - pyi); yj<std::min(height(),yi + pyi + 1); ++yj)
|
||
|
|
for (int xj = std::max(0,xi - pxi); xj<std::min(width(),xi + pxi + 1); ++xj)
|
||
|
|
if (cimg::abs(P(xi,yi) - P(xj,yj))/sig<3.) {
|
||
|
|
double d = 0;
|
||
|
|
int n = 0;
|
||
|
|
if (!(xi==xj && yi==yj)) //{
|
||
|
|
for (int ky = -patch_size; ky<patch_size + 1; ky+=sampling) {
|
||
|
|
int
|
||
|
|
yj_ = yj + ky,
|
||
|
|
yi_ = yi + ky;
|
||
|
|
if (yj_>=0 && yj_<height() && yi_>=0 && yi_<height())
|
||
|
|
for (int kx = -patch_size; kx<patch_size + 1; kx+=sampling) {
|
||
|
|
int
|
||
|
|
xj_ = xj + kx,
|
||
|
|
xi_ = xi + kx;
|
||
|
|
if (xj_>=0 && xj_<width() && xi_>=0 && xi_<width())
|
||
|
|
cimg_forC(*this,v) {
|
||
|
|
double d1 = (*this)(xj_,yj_,v) - (*this)(xi_,yi_,v);
|
||
|
|
d+=d1*d1;
|
||
|
|
n++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
//}
|
||
|
|
float w = (float)std::exp(d*h2);
|
||
|
|
cimg_forC(*this,v) uhat[v]+=w*(*this)(xj,yj,v);
|
||
|
|
wmax = w>wmax?w:wmax; // Store the maximum of the weights
|
||
|
|
sw+=w; // Compute the sum of the weights
|
||
|
|
}
|
||
|
|
}
|
||
|
|
// add the central pixel with the maximum weight
|
||
|
|
cimg_forC(*this,v) uhat[v]+=wmax*(*this)(xi,yi,v);
|
||
|
|
sw+=wmax;
|
||
|
|
|
||
|
|
// Compute the estimate for the current pixel
|
||
|
|
if (sw) cimg_forC(*this,v) dest(xi,yi,v) = (T)(uhat[v]/=sw);
|
||
|
|
else cimg_forC(*this,v) dest(xi,yi,v) = (*this)(xi,yi,v);
|
||
|
|
}
|
||
|
|
} // main loop
|
||
|
|
} // 2d
|
||
|
|
delete [] uhat;
|
||
|
|
dest.move_to(*this);
|
||
|
|
#if cimg_debug>=1
|
||
|
|
std::fprintf(stderr,"\n"); // make a new line
|
||
|
|
#endif
|
||
|
|
} // is empty
|
||
|
|
return *this;
|
||
|
|
}
|
||
|
|
|
||
|
|
//! Get the result of the NL-Means denoising algorithm.
|
||
|
|
/**
|
||
|
|
\param patch_size = radius of the patch (1=3x3 by default)
|
||
|
|
\param lambda = bandwidth ( -1 by default : automatic selection)
|
||
|
|
\param alpha = size of the region where similar patch are searched (3 x patch_size = 9x9 by default)
|
||
|
|
\param sigma = noise standard deviation (-1 = estimation)
|
||
|
|
\param sampling = sampling of the patch (1 = uses all point, 2 = uses one point on 4, etc)
|
||
|
|
If the image has three dimensions then the patch is only in 2D and the neighborhood extent in time is only 5.
|
||
|
|
If the image has several channel (color images), the distance between the two patch is computed using
|
||
|
|
all the channels.
|
||
|
|
The greater the patch is the best is the result.
|
||
|
|
Lambda parameter is function of the size of the patch size. The automatic Lambda parameter is taken
|
||
|
|
in the Chi2 table at a significiance level of 0.01. This diffear from the original paper [1].
|
||
|
|
The weighted average becomes then:
|
||
|
|
\f$$ \hat{f}(x,y) = \sum_{x',y'} \frac{1}{Z} exp(\frac{P(x,y)-P(x',y')}{2 \lambda \sigma^2}) f(x',y') $$\f
|
||
|
|
where \f$ P(x,y) $\f denotes the patch in (x,y) location.
|
||
|
|
|
||
|
|
An a priori is also used to increase the speed of the algorithm in the spirit of Sapiro et al. SPletter dec 05
|
||
|
|
|
||
|
|
This very basic version of the Non-Local Means algorithm provides an output image which contains
|
||
|
|
some residual noise with a relatively small variance (\f$\sigma<5$\f).
|
||
|
|
|
||
|
|
[1] A non-local algorithm for image denoising
|
||
|
|
Buades, A.; Coll, B.; Morel, J.-M.;
|
||
|
|
Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on
|
||
|
|
Volume 2, 20-25 June 2005 Page(s):60 - 65 vol. 2
|
||
|
|
**/
|
||
|
|
CImg<T> get_nlmeans( int patch_size=1, double lambda=-1, double alpha=3 ,double sigma=-1, int sampling=1) const {
|
||
|
|
return CImg<T>(*this).nlmeans(patch_size,lambda,alpha,sigma,sampling);
|
||
|
|
}
|
||
|
|
|
||
|
|
#endif /* cimg_plugin_nlmeans */
|