1 /* -*- C++ -*-
2 *
3 * This file is part of RawTherapee.
4 *
5 * RawTherapee is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * RawTherapee is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with RawTherapee. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #include "improcfun.h"
20 #include "imagesource.h"
21 #include "mytime.h"
22 #include "rt_algo.h"
23 #include "ipdenoise.h"
24 #ifdef _OPENMP
25 #include <omp.h>
26 #endif
27
28 namespace rtengine {
29
30 extern const Settings *settings;
31 using namespace procparams;
32
33 namespace {
34
adjust_params(procparams::DenoiseParams & dnparams,double scale)35 void adjust_params(procparams::DenoiseParams &dnparams, double scale)
36 {
37 if (scale <= 1.0) {
38 return;
39 }
40
41 double scale_factor = 1.0 / scale;
42 double noise_factor_c = std::pow(scale_factor, 0.46);
43 double noise_factor_l = std::pow(scale_factor, 0.62);
44 dnparams.luminance *= noise_factor_l;
45 dnparams.luminanceDetail *= (1.0 + std::pow(1.0 - scale_factor, 2.2));
46 dnparams.chrominance *= noise_factor_c;
47 dnparams.chrominanceRedGreen *= noise_factor_c;
48 dnparams.chrominanceBlueYellow *= noise_factor_c;
49 }
50
51
calcautodn_info(const ProcParams * params,float & chaut,float & delta,int Nb,int levaut,float maxmax,float lumema,float chromina,int mode,int lissage,float redyel,float skinc,float nsknc)52 void calcautodn_info(const ProcParams *params, float &chaut, float &delta, int Nb, int levaut, float maxmax, float lumema, float chromina, int mode, int lissage, float redyel, float skinc, float nsknc)
53 {
54
55 float reducdelta = 1.f;
56
57 if (params->denoise.aggressive) {
58 reducdelta = static_cast<float>(settings->nrhigh);
59 }
60
61 chaut = (chaut * Nb - maxmax) / (Nb - 1); //suppress maximum for chaut calcul
62
63 if ((redyel > 5000.f || skinc > 1000.f) && nsknc < 0.4f && chromina > 3000.f) {
64 chaut *= 0.45f; //reduct action in red zone, except skin for high / med chroma
65 } else if ((redyel > 12000.f || skinc > 1200.f) && nsknc < 0.3f && chromina > 3000.f) {
66 chaut *= 0.3f;
67 }
68
69 if (mode == 0 || mode == 2) { //Preview or Auto multizone
70 if (chromina > 10000.f) {
71 chaut *= 0.7f; //decrease action for high chroma (visible noise)
72 } else if (chromina > 6000.f) {
73 chaut *= 0.9f;
74 } else if (chromina < 3000.f) {
75 chaut *= 1.2f; //increase action in low chroma==> 1.2 /==>2.0 ==> curve CC
76 } else if (chromina < 2000.f) {
77 chaut *= 1.5f; //increase action in low chroma==> 1.5 / ==>2.7
78 }
79
80 if (lumema < 2500.f) {
81 chaut *= 1.3f; //increase action for low light
82 } else if (lumema < 5000.f) {
83 chaut *= 1.2f;
84 } else if (lumema > 20000.f) {
85 chaut *= 0.9f; //decrease for high light
86 }
87 } else if (mode == 1) {//auto ==> less coefficient because interaction
88 if (chromina > 10000.f) {
89 chaut *= 0.8f; //decrease action for high chroma (visible noise)
90 } else if (chromina > 6000.f) {
91 chaut *= 0.9f;
92 } else if (chromina < 3000.f) {
93 chaut *= 1.5f; //increase action in low chroma
94 } else if (chromina < 2000.f) {
95 chaut *= 2.2f; //increase action in low chroma
96 }
97
98 if (lumema < 2500.f) {
99 chaut *= 1.2f; //increase action for low light
100 } else if (lumema < 5000.f) {
101 chaut *= 1.1f;
102 } else if (lumema > 20000.f) {
103 chaut *= 0.9f; //decrease for high light
104 }
105 }
106
107 if (levaut == 0) { //Low denoise
108 if (chaut > 300.f) {
109 chaut = 0.714286f * chaut + 85.71428f;
110 }
111 }
112
113 delta = maxmax - chaut;
114 delta *= reducdelta;
115
116 if (lissage == 1 || lissage == 2) {
117 if (chaut < 200.f && delta < 200.f) {
118 delta *= 0.95f;
119 } else if (chaut < 200.f && delta < 400.f) {
120 delta *= 0.5f;
121 } else if (chaut < 200.f && delta >= 400.f) {
122 delta = 200.f;
123 } else if (chaut < 400.f && delta < 400.f) {
124 delta *= 0.4f;
125 } else if (chaut < 400.f && delta >= 400.f) {
126 delta = 120.f;
127 } else if (chaut < 550.f) {
128 delta *= 0.15f;
129 } else if (chaut < 650.f) {
130 delta *= 0.1f;
131 } else { /*if (chaut >= 650.f)*/
132 delta *= 0.07f;
133 }
134
135 if (mode == 0 || mode == 2) { //Preview or Auto multizone
136 if (chromina < 6000.f) {
137 delta *= 1.4f; //increase maxi
138 }
139
140 if (lumema < 5000.f) {
141 delta *= 1.4f;
142 }
143 } else if (mode == 1) { //Auto
144 if (chromina < 6000.f) {
145 delta *= 1.2f; //increase maxi
146 }
147
148 if (lumema < 5000.f) {
149 delta *= 1.2f;
150 }
151 }
152 }
153
154 if (lissage == 0) {
155 if (chaut < 200.f && delta < 200.f) {
156 delta *= 0.95f;
157 } else if (chaut < 200.f && delta < 400.f) {
158 delta *= 0.7f;
159 } else if (chaut < 200.f && delta >= 400.f) {
160 delta = 280.f;
161 } else if (chaut < 400.f && delta < 400.f) {
162 delta *= 0.6f;
163 } else if (chaut < 400.f && delta >= 400.f) {
164 delta = 200.f;
165 } else if (chaut < 550.f) {
166 delta *= 0.3f;
167 } else if (chaut < 650.f) {
168 delta *= 0.2f;
169 } else { /*if (chaut >= 650.f)*/
170 delta *= 0.15f;
171 }
172
173 if (mode == 0 || mode == 2) { //Preview or Auto multizone
174 if (chromina < 6000.f) {
175 delta *= 1.4f; //increase maxi
176 }
177
178 if (lumema < 5000.f) {
179 delta *= 1.4f;
180 }
181 } else if (mode == 1) { //Auto
182 if (chromina < 6000.f) {
183 delta *= 1.2f; //increase maxi
184 }
185
186 if (lumema < 5000.f) {
187 delta *= 1.2f;
188 }
189 }
190 }
191
192 }
193
194
RGB_denoise_infoGamCurve(const procparams::DenoiseParams & dnparams,bool isRAW,LUTf & gamcurve,float & gam,float & gamthresh,float & gamslope)195 void RGB_denoise_infoGamCurve(const procparams::DenoiseParams & dnparams, bool isRAW, LUTf &gamcurve, float &gam, float &gamthresh, float &gamslope)
196 {
197 gam = dnparams.gamma;
198 gamthresh = 0.001f;
199
200 if (!isRAW) {//reduce gamma under 1 for Lab mode ==> TIF and JPG
201 if (gam < 1.9f) {
202 gam = 1.f - (1.9f - gam) / 3.f; //minimum gamma 0.7
203 } else if (gam >= 1.9f && gam <= 3.f) {
204 gam = (1.4f / 1.1f) * gam - 1.41818f;
205 }
206 }
207
208 gamslope = exp(log(static_cast<double>(gamthresh)) / gam) / gamthresh;
209 Color::gammaf2lut(gamcurve, gam, gamthresh, gamslope, 65535.f, 32768.f);
210 }
211
212
RGB_denoise_info(ImProcData & im,Imagefloat * src,Imagefloat * provicalc,const bool isRAW,LUTf & gamcurve,float gam,float gamthresh,float gamslope,const procparams::DenoiseParams & dnparams,const double expcomp,float & chaut,int & Nb,float & redaut,float & blueaut,float & maxredaut,float & maxblueaut,float & minredaut,float & minblueaut,float & chromina,float & sigma,float & lumema,float & sigma_L,float & redyel,float & skinc,float & nsknc)213 void RGB_denoise_info(ImProcData &im, Imagefloat * src, Imagefloat * provicalc, const bool isRAW, LUTf &gamcurve, float gam, float gamthresh, float gamslope, const procparams::DenoiseParams & dnparams, const double expcomp, float &chaut, int &Nb, float &redaut, float &blueaut, float &maxredaut, float &maxblueaut, float &minredaut, float &minblueaut, float &chromina, float &sigma, float &lumema, float &sigma_L, float &redyel, float &skinc, float &nsknc)
214 {
215 const ProcParams *params = im.params;
216 double scale = im.scale;
217 bool multiThread = im.multiThread;
218
219 if (dnparams.chrominanceMethod != procparams::DenoiseParams::ChrominanceMethod::AUTOMATIC) {
220 //nothing to do
221 return;
222 }
223
224 int hei, wid;
225 float** lumcalc;
226 float** acalc;
227 float** bcalc;
228 hei = provicalc->getHeight();
229 wid = provicalc->getWidth();
230 TMatrix wprofi = ICCStore::getInstance()->workingSpaceMatrix(params->icm.workingProfile);
231
232 const float wpi[3][3] = {
233 {static_cast<float>(wprofi[0][0]), static_cast<float>(wprofi[0][1]), static_cast<float>(wprofi[0][2])},
234 {static_cast<float>(wprofi[1][0]), static_cast<float>(wprofi[1][1]), static_cast<float>(wprofi[1][2])},
235 {static_cast<float>(wprofi[2][0]), static_cast<float>(wprofi[2][1]), static_cast<float>(wprofi[2][2])}
236 };
237
238 lumcalc = new float*[hei];
239
240 for (int i = 0; i < hei; ++i) {
241 lumcalc[i] = new float[wid];
242 }
243
244 acalc = new float*[hei];
245
246 for (int i = 0; i < hei; ++i) {
247 acalc[i] = new float[wid];
248 }
249
250 bcalc = new float*[hei];
251
252 for (int i = 0; i < hei; ++i) {
253 bcalc[i] = new float[wid];
254 }
255
256 #ifdef _OPENMP
257 #pragma omp parallel for if (multiThread)
258 #endif
259
260 for (int ii = 0; ii < hei; ++ii) {
261 for (int jj = 0; jj < wid; ++jj) {
262 float LLum, AAum, BBum;
263 float RL = provicalc->r(ii, jj);
264 float GL = provicalc->g(ii, jj);
265 float BL = provicalc->b(ii, jj);
266 // determine luminance for noisecurve
267 float XL, YL, ZL;
268 Color::rgbxyz(RL, GL, BL, XL, YL, ZL, wpi);
269 Color::XYZ2Lab(XL, YL, ZL, LLum, AAum, BBum);
270 lumcalc[ii][jj] = LLum;
271 acalc[ii][jj] = AAum;
272 bcalc[ii][jj] = BBum;
273 }
274 }
275
276 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
277
278 const int imheight = src->getHeight(), imwidth = src->getWidth();
279 const float gain = pow(2.0f, float(expcomp));
280
281 int tilesize;
282 int overlap;
283
284 if (settings->leveldnti == 0) {
285 tilesize = 1024;
286 overlap = 128;
287 }
288
289 if (settings->leveldnti == 1) {
290 tilesize = 768;
291 overlap = 96;
292 }
293
294 int numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip;
295
296 //always no Tiles
297 int kall = 0;
298 rtengine::denoise::Tile_calc(tilesize, overlap, kall, imwidth, imheight, numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip);
299
300 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
301
302 TMatrix wprof = ICCStore::getInstance()->workingSpaceMatrix(params->icm.workingProfile);
303 const float wp[3][3] = {
304 {static_cast<float>(wprof[0][0]), static_cast<float>(wprof[0][1]), static_cast<float>(wprof[0][2])},
305 {static_cast<float>(wprof[1][0]), static_cast<float>(wprof[1][1]), static_cast<float>(wprof[1][2])},
306 {static_cast<float>(wprof[2][0]), static_cast<float>(wprof[2][1]), static_cast<float>(wprof[2][2])}
307 };
308
309 float chau = 0.f;
310 float chred = 0.f;
311 float chblue = 0.f;
312 float maxchred = 0.f;
313 float maxchblue = 0.f;
314 float minchred = 100000000.f;
315 float minchblue = 100000000.f;
316 int nb = 0;
317 int comptlevel = 0;
318
319 for (int tiletop = 0; tiletop < imheight; tiletop += tileHskip) {
320 for (int tileleft = 0; tileleft < imwidth; tileleft += tileWskip) {
321
322 int tileright = MIN(imwidth, tileleft + tilewidth);
323 int tilebottom = MIN(imheight, tiletop + tileheight);
324 int width = tileright - tileleft;
325 int height = tilebottom - tiletop;
326 LabImage * labdn = new LabImage(width, height);
327 float** noisevarlum = new float*[(height + 1) / 2];
328
329 for (int i = 0; i < (height + 1) / 2; ++i) {
330 noisevarlum[i] = new float[(width + 1) / 2];
331 }
332
333 float** noisevarchrom = new float*[(height + 1) / 2];
334
335 for (int i = 0; i < (height + 1) / 2; ++i) {
336 noisevarchrom[i] = new float[(width + 1) / 2];
337 }
338
339 float** noisevarhue = new float*[(height + 1) / 2];
340
341 for (int i = 0; i < (height + 1) / 2; ++i) {
342 noisevarhue[i] = new float[(width + 1) / 2];
343 }
344
345 float realred, realblue;
346 float interm_med = static_cast<float>(dnparams.chrominance) / 10.0;
347 float intermred, intermblue;
348
349 if (dnparams.chrominanceRedGreen > 0.) {
350 intermred = (dnparams.chrominanceRedGreen / 10.);
351 } else {
352 intermred = static_cast<float>(dnparams.chrominanceRedGreen) / 7.0; //increase slower than linear for more sensit
353 }
354
355 if (dnparams.chrominanceBlueYellow > 0.) {
356 intermblue = (dnparams.chrominanceBlueYellow / 10.);
357 } else {
358 intermblue = static_cast<float>(dnparams.chrominanceBlueYellow) / 7.0; //increase slower than linear for more sensit
359 }
360
361 realred = interm_med + intermred;
362
363 if (realred < 0.f) {
364 realred = 0.001f;
365 }
366
367 realblue = interm_med + intermblue;
368
369 if (realblue < 0.f) {
370 realblue = 0.001f;
371 }
372
373 //fill tile from image; convert RGB to "luma/chroma"
374
375 if (isRAW) {//image is raw; use channel differences for chroma channels
376 #ifdef _OPENMP
377 #pragma omp parallel for if (multiThread)
378 #endif
379
380 for (int i = tiletop; i < tilebottom; i += 2) {
381 int i1 = i - tiletop;
382 #ifdef __SSE2__
383 __m128 aNv, bNv;
384 __m128 c100v = _mm_set1_ps(100.f);
385 int j;
386
387 for (j = tileleft; j < tileright - 7; j += 8) {
388 int j1 = j - tileleft;
389 aNv = LVFU(acalc[i >> 1][j >> 1]);
390 bNv = LVFU(bcalc[i >> 1][j >> 1]);
391 _mm_storeu_ps(&noisevarhue[i1 >> 1][j1 >> 1], xatan2f(bNv, aNv));
392 _mm_storeu_ps(&noisevarchrom[i1 >> 1][j1 >> 1], vmaxf(vsqrtf(SQRV(aNv) + SQRV(bNv)),c100v));
393 }
394
395 for (; j < tileright; j += 2) {
396 int j1 = j - tileleft;
397 float aN = acalc[i >> 1][j >> 1];
398 float bN = bcalc[i >> 1][j >> 1];
399 float cN = sqrtf(SQR(aN) + SQR(bN));
400 noisevarhue[i1 >> 1][j1 >> 1] = xatan2f(bN, aN);
401
402 if (cN < 100.f) {
403 cN = 100.f; //avoid divided by zero
404 }
405
406 noisevarchrom[i1 >> 1][j1 >> 1] = cN;
407 }
408
409 #else
410
411 for (int j = tileleft; j < tileright; j += 2) {
412 int j1 = j - tileleft;
413 float aN = acalc[i >> 1][j >> 1];
414 float bN = bcalc[i >> 1][j >> 1];
415 float cN = sqrtf(SQR(aN) + SQR(bN));
416 float hN = xatan2f(bN, aN);
417
418 if (cN < 100.f) {
419 cN = 100.f; //avoid divided by zero
420 }
421
422 noisevarchrom[i1 >> 1][j1 >> 1] = cN;
423 noisevarhue[i1 >> 1][j1 >> 1] = hN;
424 }
425
426 #endif
427 }
428
429 #ifdef _OPENMP
430 #pragma omp parallel for if (multiThread)
431 #endif
432
433 for (int i = tiletop; i < tilebottom; i += 2) {
434 int i1 = i - tiletop;
435
436 for (int j = tileleft; j < tileright; j += 2) {
437 int j1 = j - tileleft;
438 float Llum = lumcalc[i >> 1][j >> 1];
439 Llum = Llum < 2.f ? 2.f : Llum; //avoid divided by zero ?
440 Llum = Llum > 32768.f ? 32768.f : Llum; // not strictly necessary
441 noisevarlum[i1 >> 1][j1 >> 1] = Llum;
442 }
443 }
444
445 for (int i = tiletop/*, i1=0*/; i < tilebottom; ++i/*, ++i1*/) {
446 int i1 = i - tiletop;
447
448 for (int j = tileleft/*, j1=0*/; j < tileright; ++j/*, ++j1*/) {
449 int j1 = j - tileleft;
450
451 float X = gain * src->r(i, j);
452 float Y = gain * src->g(i, j);
453 float Z = gain * src->b(i, j);
454
455 X = X < 65535.f ? gamcurve[X] : (Color::gammaf(X / 65535.f, gam, gamthresh, gamslope) * 32768.f);
456 Y = Y < 65535.f ? gamcurve[Y] : (Color::gammaf(Y / 65535.f, gam, gamthresh, gamslope) * 32768.f);
457 Z = Z < 65535.f ? gamcurve[Z] : (Color::gammaf(Z / 65535.f, gam, gamthresh, gamslope) * 32768.f);
458
459 // labdn->a[i1][j1] = (X - Y);
460 // labdn->b[i1][j1] = (Y - Z);
461 float l, u, v;
462 Color::rgb2yuv(X, Y, Z, l, u, v, wp);
463 labdn->a[i1][j1] = v;
464 labdn->b[i1][j1] = u;
465 }
466 }
467 } else {//image is not raw; use Lab parametrization
468 for (int i = tiletop/*, i1=0*/; i < tilebottom; ++i/*, ++i1*/) {
469 int i1 = i - tiletop;
470
471 for (int j = tileleft/*, j1=0*/; j < tileright; ++j/*, ++j1*/) {
472 int j1 = j - tileleft;
473 //float L, a, b;
474 float rLum = src->r(i, j) ; //for luminance denoise curve
475 float gLum = src->g(i, j) ;
476 float bLum = src->b(i, j) ;
477
478 //use gamma sRGB, not good if TIF (JPG) Output profil not with gamma sRGB (eg : gamma =1.0, or 1.8...)
479 //very difficult to solve !
480 // solution ==> save TIF with gamma sRGB and re open
481 float rtmp = Color::igammatab_srgb[ src->r(i, j) ];
482 float gtmp = Color::igammatab_srgb[ src->g(i, j) ];
483 float btmp = Color::igammatab_srgb[ src->b(i, j) ];
484 //modification Jacques feb 2013
485 // gamma slider different from raw
486 rtmp = rtmp < 65535.f ? gamcurve[rtmp] : (Color::gammanf(rtmp / 65535.f, gam) * 32768.f);
487 gtmp = gtmp < 65535.f ? gamcurve[gtmp] : (Color::gammanf(gtmp / 65535.f, gam) * 32768.f);
488 btmp = btmp < 65535.f ? gamcurve[btmp] : (Color::gammanf(btmp / 65535.f, gam) * 32768.f);
489
490 // float X, Y, Z;
491 // Color::rgbxyz(rtmp, gtmp, btmp, X, Y, Z, wp);
492
493 // //convert Lab
494 // Color::XYZ2Lab(X, Y, Z, L, a, b);
495 float Y, u, v;
496 Color::rgb2yuv(rtmp, gtmp, btmp, Y, u, v, wp);
497
498 if (((i1 | j1) & 1) == 0) {
499 float Llum, alum, blum;
500 float XL, YL, ZL;
501 Color::rgbxyz(rLum, gLum, bLum, XL, YL, ZL, wp);
502 Color::XYZ2Lab(XL, YL, ZL, Llum, alum, blum);
503 float kN = Llum;
504
505 if (kN < 2.f) {
506 kN = 2.f;
507 }
508
509 if (kN > 32768.f) {
510 kN = 32768.f;
511 }
512
513 noisevarlum[i1 >> 1][j1 >> 1] = kN;
514 float aN = alum;
515 float bN = blum;
516 float hN = xatan2f(bN, aN);
517 float cN = sqrt(SQR(aN) + SQR(bN));
518
519 if (cN < 100.f) {
520 cN = 100.f; //avoid divided by zero
521 }
522
523 noisevarchrom[i1 >> 1][j1 >> 1] = cN;
524 noisevarhue[i1 >> 1][j1 >> 1] = hN;
525 }
526
527 labdn->a[i1][j1] = v;
528 labdn->b[i1][j1] = u;
529 }
530 }
531 }
532
533 int datalen = labdn->W * labdn->H;
534
535 //now perform basic wavelet denoise
536 //last two arguments of wavelet decomposition are max number of wavelet decomposition levels;
537 //and whether to subsample the image after wavelet filtering. Subsampling is coded as
538 //binary 1 or 0 for each level, eg subsampling = 0 means no subsampling, 1 means subsample
539 //the first level only, 7 means subsample the first three levels, etc.
540
541 wavelet_decomposition* adecomp;
542 wavelet_decomposition* bdecomp;
543
544 int schoice = 0;//shrink method
545
546 if (dnparams.aggressive) {
547 schoice = 2;
548 }
549
550 const int levwav = max(2, int(5 - std::ceil(std::log(scale))));
551 #ifdef _OPENMP
552 #pragma omp parallel sections if (multiThread)
553 #endif
554 {
555 #ifdef _OPENMP
556 #pragma omp section
557 #endif
558 {
559 adecomp = new wavelet_decomposition(labdn->data + datalen, labdn->W, labdn->H, levwav, 1);
560 }
561 #ifdef _OPENMP
562 #pragma omp section
563 #endif
564 {
565 bdecomp = new wavelet_decomposition(labdn->data + 2 * datalen, labdn->W, labdn->H, levwav, 1);
566 }
567 }
568
569 if (comptlevel == 0) {
570 denoise::WaveletDenoiseAll_info(
571 levwav,
572 *adecomp,
573 *bdecomp,
574 noisevarlum,
575 noisevarchrom,
576 noisevarhue,
577 chaut,
578 Nb,
579 redaut,
580 blueaut,
581 maxredaut,
582 maxblueaut,
583 minredaut,
584 minblueaut,
585 schoice,
586 chromina,
587 sigma,
588 lumema,
589 sigma_L,
590 redyel,
591 skinc,
592 nsknc,
593 maxchred,
594 maxchblue,
595 minchred,
596 minchblue,
597 nb,
598 chau,
599 chred,
600 chblue
601 ); // Enhance mode
602 }
603
604 comptlevel += 1;
605 delete adecomp;
606 delete bdecomp;
607 delete labdn;
608
609 for (int i = 0; i < (height + 1) / 2; ++i) {
610 delete[] noisevarlum[i];
611 }
612
613 delete[] noisevarlum;
614
615 for (int i = 0; i < (height + 1) / 2; ++i) {
616 delete[] noisevarchrom[i];
617 }
618
619 delete[] noisevarchrom;
620
621 for (int i = 0; i < (height + 1) / 2; ++i) {
622 delete[] noisevarhue[i];
623 }
624
625 delete[] noisevarhue;
626
627 }//end of tile row
628 }//end of tile loop
629
630 for (int i = 0; i < hei; ++i) {
631 delete[] lumcalc[i];
632 }
633
634 delete[] lumcalc;
635
636 for (int i = 0; i < hei; ++i) {
637 delete[] acalc[i];
638 }
639
640 delete[] acalc;
641
642 for (int i = 0; i < hei; ++i) {
643 delete[] bcalc[i];
644 }
645
646 delete[] bcalc;
647
648 #undef TS
649 //#undef fTS
650 #undef offset
651 #undef epsilon
652
653 } // End of main RGB_denoise
654
655 } // namespace
656
657
denoiseComputeParams(ImageSource * imgsrc,const ColorTemp & currWB,DenoiseInfoStore & store,procparams::DenoiseParams & dnparams)658 void ImProcFunctions::denoiseComputeParams(ImageSource *imgsrc, const ColorTemp &currWB, DenoiseInfoStore &store, procparams::DenoiseParams &dnparams)
659 {
660 float autoNR = settings->nrauto;
661 float autoNRmax = settings->nrautomax;
662
663 int tilesize;
664 int overlap;
665
666 if (settings->leveldnti == 0) {
667 tilesize = 1024;
668 overlap = 128;
669 }
670
671 if (settings->leveldnti == 1) {
672 tilesize = 768;
673 overlap = 96;
674 }
675
676 int numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip;
677 int kall = 2;
678
679 int widIm, heiIm;
680 int tr = getCoarseBitMask(params->coarse);
681 imgsrc->getFullSize(widIm, heiIm, tr);
682
683 denoise::Tile_calc(tilesize, overlap, kall, widIm, heiIm, numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip);
684 kall = 0;
685
686 float min_b[9];
687 float min_r[9];
688 float lumL[9];
689 float chromC[9];
690 float ry[9];
691 float sk[9];
692 float pcsk[9];
693 std::vector<int> centerTile_X(numtiles_W);
694 std::vector<int> centerTile_Y(numtiles_H);
695
696 for (int cX = 0; cX < numtiles_W; cX++) {
697 centerTile_X[cX] = tileWskip / 2 + tileWskip * cX;
698 }
699
700 for (int cY = 0; cY < numtiles_H; cY++) {
701 centerTile_Y[cY] = tileHskip / 2 + tileHskip * cY;
702 }
703
704 assert(settings->leveldnautsimpl == 0);
705
706 if (!store.valid && dnparams.chrominanceMethod == procparams::DenoiseParams::ChrominanceMethod::AUTOMATIC) {
707 MyTime t1aue, t2aue;
708 t1aue.set();
709
710 int crW = 100; // settings->leveldnv == 0
711 int crH = 100; // settings->leveldnv == 0
712
713 if (settings->leveldnv == 1) {
714 crW = 250;
715 crH = 250;
716 }
717
718 // if(settings->leveldnv ==2) {crW=int(tileWskip/2);crH=int((tileWskip/2));}//adapted to scale of preview
719 if (settings->leveldnv == 2) {
720 crW = int (tileWskip / 2);
721 crH = int (tileHskip / 2);
722 }
723
724 if (settings->leveldnv == 3) {
725 crW = tileWskip - 10;
726 crH = tileHskip - 10;
727 }
728
729 float lowdenoise = 1.f;
730 int levaut = settings->leveldnaut;
731
732 if (levaut == 1) { //Standard
733 lowdenoise = 0.7f;
734 }
735
736 LUTf gamcurve(65536, 0);
737 float gam, gamthresh, gamslope;
738 RGB_denoise_infoGamCurve(dnparams, imgsrc->isRAW(), gamcurve, gam, gamthresh, gamslope);
739 int Nb[9];
740
741 #ifdef _OPENMP
742 #pragma omp parallel if (multiThread)
743 #endif
744 {
745 Imagefloat *origCropPart = new Imagefloat(crW, crH); //allocate memory
746 Imagefloat *provicalc = new Imagefloat((crW + 1) / 2, (crH + 1) / 2); //for denoise curves
747
748 int coordW[3];//coordinate of part of image to measure noise
749 int coordH[3];
750 int begW = 50;
751 int begH = 50;
752 coordW[0] = begW;
753 coordW[1] = widIm / 2 - crW / 2;
754 coordW[2] = widIm - crW - begW;
755 coordH[0] = begH;
756 coordH[1] = heiIm / 2 - crH / 2;
757 coordH[2] = heiIm - crH - begH;
758
759 #ifdef _OPENMP
760 # pragma omp for schedule(dynamic) collapse(2) nowait
761 #endif
762
763 for (int wcr = 0; wcr <= 2; wcr++) {
764 for (int hcr = 0; hcr <= 2; hcr++) {
765 PreviewProps ppP(coordW[wcr], coordH[hcr], crW, crH, 1);
766 imgsrc->getImage(currWB, tr, origCropPart, ppP, params->exposure, params->raw);
767
768 // we only need image reduced to 1/4 here
769 for (int ii = 0; ii < crH; ii += 2) {
770 for (int jj = 0; jj < crW; jj += 2) {
771 provicalc->r(ii >> 1, jj >> 1) = origCropPart->r(ii, jj);
772 provicalc->g(ii >> 1, jj >> 1) = origCropPart->g(ii, jj);
773 provicalc->b(ii >> 1, jj >> 1) = origCropPart->b(ii, jj);
774 }
775 }
776
777 imgsrc->convertColorSpace(provicalc, params->icm, currWB); //for denoise luminance curve
778
779 float pondcorrec = 1.0f;
780 float chaut = 0.f, redaut = 0.f, blueaut = 0.f, maxredaut = 0.f, maxblueaut = 0.f, minredaut = 0.f, minblueaut = 0.f, chromina = 0.f, sigma = 0.f, lumema = 0.f, sigma_L = 0.f, redyel = 0.f, skinc = 0.f, nsknc = 0.f;
781 int nb = 0;
782 ImProcData im(params, scale, multiThread);
783 RGB_denoise_info(im, origCropPart, provicalc, imgsrc->isRAW(), gamcurve, gam, gamthresh, gamslope, dnparams, std::log(5.f)/std::log(2.f)/*imgsrc->getDirPyrDenoiseExpComp()*/, chaut, nb, redaut, blueaut, maxredaut, maxblueaut, minredaut, minblueaut, chromina, sigma, lumema, sigma_L, redyel, skinc, nsknc);
784
785 //printf("DCROP skip=%d cha=%f red=%f bl=%f redM=%f bluM=%f chrom=%f sigm=%f lum=%f\n",skip, chaut,redaut,blueaut, maxredaut, maxblueaut, chromina, sigma, lumema);
786 Nb[hcr * 3 + wcr] = nb;
787 store.ch_M[hcr * 3 + wcr] = pondcorrec * chaut;
788 store.max_r[hcr * 3 + wcr] = pondcorrec * maxredaut;
789 store.max_b[hcr * 3 + wcr] = pondcorrec * maxblueaut;
790 min_r[hcr * 3 + wcr] = pondcorrec * minredaut;
791 min_b[hcr * 3 + wcr] = pondcorrec * minblueaut;
792 lumL[hcr * 3 + wcr] = lumema;
793 chromC[hcr * 3 + wcr] = chromina;
794 ry[hcr * 3 + wcr] = redyel;
795 sk[hcr * 3 + wcr] = skinc;
796 pcsk[hcr * 3 + wcr] = nsknc;
797
798 }
799 }
800
801 delete provicalc;
802 delete origCropPart;
803 }
804 float chM = 0.f;
805 float MaxR = 0.f;
806 float MaxB = 0.f;
807 float MinR = 100000000000.f;
808 float MinB = 100000000000.f;
809 float maxr = 0.f;
810 float maxb = 0.f;
811 float Max_R[9] = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f};
812 float Max_B[9] = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f};
813 float Min_R[9];
814 float Min_B[9];
815 float MaxRMoy = 0.f;
816 float MaxBMoy = 0.f;
817 float MinRMoy = 0.f;
818 float MinBMoy = 0.f;
819
820 float multip = 1.f;
821
822 if (!imgsrc->isRAW()) {
823 multip = 2.f; //take into account gamma for TIF / JPG approximate value...not good for gamma=1
824 }
825
826 float adjustr = 1.f;
827
828 // if (params->icm.workingProfile == "ProPhoto") {
829 // adjustr = 1.f; //
830 // } else if (params->icm.workingProfile == "Adobe RGB") {
831 // adjustr = 1.f / 1.3f;
832 // } else if (params->icm.workingProfile == "sRGB") {
833 // adjustr = 1.f / 1.3f;
834 // } else if (params->icm.workingProfile == "WideGamut") {
835 // adjustr = 1.f / 1.1f;
836 // } else if (params->icm.workingProfile == "Beta RGB") {
837 // adjustr = 1.f / 1.2f;
838 // } else if (params->icm.workingProfile == "BestRGB") {
839 // adjustr = 1.f / 1.2f;
840 // } else if (params->icm.workingProfile == "BruceRGB") {
841 // adjustr = 1.f / 1.2f;
842 // }
843
844 float delta[9];
845 int mode = 1;
846 int lissage = settings->leveldnliss;
847
848 for (int k = 0; k < 9; k++) {
849 float maxmax = max(store.max_r[k], store.max_b[k]);
850 calcautodn_info(params, store.ch_M[k], delta[k], Nb[k], levaut, maxmax, lumL[k], chromC[k], mode, lissage, ry[k], sk[k], pcsk[k]);
851 // printf("ch_M=%f delta=%f\n",ch_M[k], delta[k]);
852 }
853
854 for (int k = 0; k < 9; k++) {
855 if (store.max_r[k] > store.max_b[k]) {
856 Max_R[k] = (delta[k]) / ((autoNRmax * multip * adjustr * lowdenoise) / 2.f);
857 Min_B[k] = - (store.ch_M[k] - min_b[k]) / (autoNRmax * multip * adjustr * lowdenoise);
858 Max_B[k] = 0.f;
859 Min_R[k] = 0.f;
860 } else {
861 Max_B[k] = (delta[k]) / ((autoNRmax * multip * adjustr * lowdenoise) / 2.f);
862 Min_R[k] = - (store.ch_M[k] - min_r[k]) / (autoNRmax * multip * adjustr * lowdenoise);
863 Min_B[k] = 0.f;
864 Max_R[k] = 0.f;
865 }
866 }
867
868 for (int k = 0; k < 9; k++) {
869 // printf("ch_M= %f Max_R=%f Max_B=%f min_r=%f min_b=%f\n",ch_M[k],Max_R[k], Max_B[k],Min_R[k], Min_B[k]);
870 chM += store.ch_M[k];
871 MaxBMoy += Max_B[k];
872 MaxRMoy += Max_R[k];
873 MinRMoy += Min_R[k];
874 MinBMoy += Min_B[k];
875
876 if (Max_R[k] > MaxR) {
877 MaxR = Max_R[k];
878 }
879
880 if (Max_B[k] > MaxB) {
881 MaxB = Max_B[k];
882 }
883
884 if (Min_R[k] < MinR) {
885 MinR = Min_R[k];
886 }
887
888 if (Min_B[k] < MinB) {
889 MinB = Min_B[k];
890 }
891 }
892
893 chM /= 9;
894 MaxBMoy /= 9;
895 MaxRMoy /= 9;
896 MinBMoy /= 9;
897 MinRMoy /= 9;
898
899 if (MaxR > MaxB) {
900 maxr = MaxRMoy + (MaxR - MaxRMoy) * 0.66f; //#std Dev
901 //maxb=MinB;
902 maxb = MinBMoy + (MinB - MinBMoy) * 0.66f;
903 } else {
904 maxb = MaxBMoy + (MaxB - MaxBMoy) * 0.66f;
905 maxr = MinRMoy + (MinR - MinRMoy) * 0.66f;
906 }
907
908 // printf("DCROP skip=%d cha=%f red=%f bl=%f \n",skip, chM,maxr,maxb);
909 dnparams.chrominance = chM / (autoNR * multip * adjustr);
910 dnparams.chrominanceRedGreen = maxr;
911 dnparams.chrominanceBlueYellow = maxb;
912
913 dnparams.chrominance *= dnparams.chrominanceAutoFactor;
914 dnparams.chrominanceRedGreen *= dnparams.chrominanceAutoFactor;
915 dnparams.chrominanceBlueYellow *= dnparams.chrominanceAutoFactor;
916
917 store.valid = true;
918
919 if (settings->verbose) {
920 t2aue.set();
921 printf("Info denoise auto performed in %d usec:\n", t2aue.etime(t1aue));
922 }
923 //end evaluate noise
924 }
925 }
926
927
denoise(ImageSource * imgsrc,const ColorTemp & currWB,Imagefloat * img,const DenoiseInfoStore & store,const procparams::DenoiseParams & dnparams)928 void ImProcFunctions::denoise(ImageSource *imgsrc, const ColorTemp &currWB, Imagefloat *img, const DenoiseInfoStore &store, const procparams::DenoiseParams &dnparams)
929 {
930 if (!dnparams.enabled) {
931 return;
932 }
933
934 if (plistener) {
935 plistener->setProgressStr("PROGRESSBAR_DENOISING");
936 plistener->setProgress(0);
937 }
938
939 procparams::DenoiseParams denoiseParams = dnparams;
940 NoiseCurve noiseLCurve;
941 NoiseCurve noiseCCurve;
942
943 const int W = img->getWidth();
944 const int H = img->getHeight();
945
946 Imagefloat *calclum = nullptr;
947 {
948 const int fw = W;
949 const int fh = H;
950 // we only need image reduced to 1/4 here
951 calclum = new Imagefloat((fw + 1) / 2, (fh + 1) / 2); //for luminance denoise curve
952 #ifdef _OPENMP
953 # pragma omp parallel for if (multiThread)
954 #endif
955
956 for (int ii = 0; ii < fh; ii += 2) {
957 for (int jj = 0; jj < fw; jj += 2) {
958 calclum->r(ii >> 1, jj >> 1) = img->r(ii, jj);
959 calclum->g(ii >> 1, jj >> 1) = img->g(ii, jj);
960 calclum->b(ii >> 1, jj >> 1) = img->b(ii, jj);
961 }
962 }
963 imgsrc->convertColorSpace(calclum, params->icm, currWB);
964 }
965
966 float nresi, highresi;
967 DenoiseInfoStore &dnstore = const_cast<DenoiseInfoStore &>(store);
968
969 adjust_params(denoiseParams, scale);
970
971 noiseCCurve.Set({
972 FCT_MinMaxCPoints,
973 0.05,
974 0.50,
975 0.35,
976 0.35,
977 0.35,
978 0.05,
979 0.35,
980 0.35
981 });
982
983 if (plistener) {
984 plistener->setProgress(0.1);
985 }
986
987 ImProcData im(params, scale, multiThread);
988 denoise::RGB_denoise(im, 0, img, img, calclum, dnstore.ch_M, dnstore.max_r, dnstore.max_b, imgsrc->isRAW(), denoiseParams, 0, noiseLCurve, noiseCCurve, nresi, highresi);
989
990 if (plistener) {
991 plistener->setProgress(0.8);
992 }
993
994 if (denoiseParams.smoothingEnabled) {
995 denoise::denoiseGuidedSmoothing(im, img);
996 if (denoiseParams.nlStrength) {
997 img->setMode(Imagefloat::Mode::YUV, multiThread);
998 array2D<float> tmp(img->getWidth(), img->getHeight(), img->g.ptrs, ARRAY2D_BYREFERENCE);
999 denoise::NLMeans(tmp, 65535.f, denoiseParams.nlStrength, denoiseParams.nlDetail, scale, multiThread);
1000 img->setMode(Imagefloat::Mode::RGB, multiThread);
1001 }
1002 }
1003
1004 if (plistener) {
1005 plistener->setProgress(1);
1006 }
1007 }
1008
1009
1010 } // namespace rtengine
1011