kpriv/freebsd-ports
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
freebsd-ports/graphics/openfx-misc/files/patch-CImg_nlmeans.h
Olivier Cochard e03ec2f3ff New port: graphics/openfx-misc 
Miscellaneous OFX / OpenFX / Open Effects plugins. These plugins were
primarily developped for Natron, but may be used with other OpenFX hosts.
2018-05-01 23:51:34 +00:00

245 lines
12 KiB
C++
Raw Blame History

--- CImg/nlmeans.h.orig 2018-04-30 23:16:26 UTC
+++ CImg/nlmeans.h
@@ -0,0 +1,242 @@
+/*
+ #
+ # 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 */
Reference in a new issue View git blame Copy permalink