1 /* -*- C++ -*-
2 *
3 * This file is part of RawTherapee.
4 *
5 * Ported from G'MIC by Alberto Griggio <alberto.griggio@gmail.com>
6 *
7 * The original implementation in G'MIC was authored by Arto Huotari, and was
8 * released under the CeCILL free software license (see
9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.html)
10 *
11 * RawTherapee is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * RawTherapee is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with RawTherapee. If not, see <http://www.gnu.org/licenses/>.
23 */
24
25 #ifdef _OPENMP
26 #include <omp.h>
27 #endif
28
29 #include "improcfun.h"
30 #include "gauss.h"
31 #include "array2D.h"
32 #include "cplx_wavelet_dec.h"
33 #include "curves.h"
34 #include "labmasks.h"
35
36
37 namespace rtengine {
38
39 namespace {
40
41 class WavOpacityCurveWL {
42 private:
43 LUTf lutOpacityCurveWL; // 0xffff range
44 void Set(const Curve &pCurve);
45
46 public:
~WavOpacityCurveWL()47 virtual ~WavOpacityCurveWL() {};
48 WavOpacityCurveWL();
49
50 void Reset();
51 void Set(const Curve *pCurve);
52 void Set(const std::vector<double> &curvePoints);
operator [](float index) const53 float operator[](float index) const
54 {
55 return lutOpacityCurveWL ? lutOpacityCurveWL[index] : 0.f;
56 }
57
operator bool(void) const58 operator bool (void) const
59 {
60 return lutOpacityCurveWL;
61 }
62 };
63
WavOpacityCurveWL()64 WavOpacityCurveWL::WavOpacityCurveWL() {}
65
Reset()66 void WavOpacityCurveWL::Reset()
67 {
68 lutOpacityCurveWL.reset();
69 }
70
Set(const Curve & pCurve)71 void WavOpacityCurveWL::Set(const Curve &pCurve)
72 {
73 if (pCurve.isIdentity()) {
74 lutOpacityCurveWL.reset(); // raise this value if the quality suffers from this number of samples
75 return;
76 }
77
78 lutOpacityCurveWL(501); // raise this value if the quality suffers from this number of samples
79
80 for (int i = 0; i < 501; i++) {
81 lutOpacityCurveWL[i] = pCurve.getVal(double(i) / 500.);
82 }
83 }
84
Set(const std::vector<double> & curvePoints)85 void WavOpacityCurveWL::Set(const std::vector<double> &curvePoints)
86 {
87 if (!curvePoints.empty() && curvePoints[0] > FCT_Linear && curvePoints[0] < FCT_Unchanged) {
88 FlatCurve tcurve(curvePoints, false, CURVES_MIN_POLY_POINTS / 2);
89 tcurve.setIdentityValue(0.);
90 Set(tcurve);
91 } else {
92 Reset();
93 }
94 }
95
96
eval_avg(float * RESTRICT DataList,int datalen,float & averagePlus,float & averageNeg,float & max,float & min,bool multiThread)97 void eval_avg(float *RESTRICT DataList, int datalen, float &averagePlus, float &averageNeg, float &max, float &min, bool multiThread)
98 {
99 //find absolute mean
100 int countP = 0, countN = 0;
101 double averaP = 0.0, averaN = 0.0; // use double precision for large summations
102
103 float thres = 5.f;//different fom zero to take into account only data large enough
104 max = 0.f;
105 min = 0.f;
106
107 #ifdef _OPENMP
108 # pragma omp parallel if (multiThread)
109 #endif
110 {
111 float lmax = 0.f, lmin = 0.f;
112 #ifdef _OPENMP
113 #pragma omp for reduction(+:averaP,averaN,countP,countN) nowait
114 #endif
115
116 for(int i = 0; i < datalen; i++) {
117 if(DataList[i] >= thres) {
118 averaP += DataList[i];
119
120 if(DataList[i] > lmax) {
121 lmax = DataList[i];
122 }
123
124 countP++;
125 } else if(DataList[i] < -thres) {
126 averaN += DataList[i];
127
128 if(DataList[i] < lmin) {
129 lmin = DataList[i];
130 }
131
132 countN++;
133 }
134 }
135
136 #ifdef _OPENMP
137 #pragma omp critical
138 #endif
139 {
140 max = max > lmax ? max : lmax;
141 min = min < lmin ? min : lmin;
142 }
143 }
144
145 if(countP > 0) {
146 averagePlus = averaP / countP;
147 } else {
148 averagePlus = 0;
149 }
150
151 if(countN > 0) {
152 averageNeg = averaN / countN;
153 } else {
154 averageNeg = 0;
155 }
156 }
157
158
eval_sigma(float * RESTRICT DataList,int datalen,float averagePlus,float averageNeg,float & sigmaPlus,float & sigmaNeg,bool multiThread)159 void eval_sigma(float *RESTRICT DataList, int datalen, float averagePlus, float averageNeg, float &sigmaPlus, float &sigmaNeg, bool multiThread)
160 {
161 int countP = 0, countN = 0;
162 double variP = 0.0, variN = 0.0; // use double precision for large summations
163 float thres = 5.f;//different fom zero to take into account only data large enough
164
165 #ifdef _OPENMP
166 # pragma omp parallel for reduction(+:variP,variN,countP,countN) if (multiThread)
167 #endif
168 for(int i = 0; i < datalen; i++) {
169 if(DataList[i] >= thres) {
170 variP += SQR(DataList[i] - averagePlus);
171 countP++;
172 } else if(DataList[i] <= -thres) {
173 variN += SQR(DataList[i] - averageNeg);
174 countN++;
175 }
176 }
177
178 if(countP > 0) {
179 sigmaPlus = sqrt(variP / countP);
180 } else {
181 sigmaPlus = 0;
182 }
183
184 if(countN > 0) {
185 sigmaNeg = sqrt(variN / countN);
186 } else {
187 sigmaNeg = 0;
188 }
189 }
190
191
eval_level(float ** wl,int level,int W_L,int H_L,float * mean,float * meanN,float * sigma,float * sigmaN,float * MaxP,float * MaxN,bool multiThread)192 void eval_level(float **wl, int level, int W_L, int H_L, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN, bool multiThread)
193 {
194 float avLP[4], avLN[4];
195 float maxL[4], minL[4];
196 float sigP[4], sigN[4];
197 float AvL, AvN, SL, SN, maxLP, maxLN;
198
199 for (int dir = 1; dir < 4; dir++) {
200 eval_avg(wl[dir], W_L * H_L, avLP[dir], avLN[dir], maxL[dir], minL[dir], multiThread);
201 eval_sigma(wl[dir], W_L * H_L, avLP[dir], avLN[dir], sigP[dir], sigN[dir], multiThread);
202 }
203
204 AvL = 0.f;
205 AvN = 0.f;
206 SL = 0.f;
207 SN = 0.f;
208 maxLP = 0.f;
209 maxLN = 0.f;
210
211 for (int dir = 1; dir < 4; dir++) {
212 AvL += avLP[dir];
213 AvN += avLN[dir];
214 SL += sigP[dir];
215 SN += sigN[dir];
216 maxLP += maxL[dir];
217 maxLN += minL[dir];
218 }
219
220 AvL /= 3;
221 AvN /= 3;
222 SL /= 3;
223 SN /= 3;
224 maxLP /= 3;
225 maxLN /= 3;
226
227 mean[level] = AvL;
228 meanN[level] = AvN;
229 sigma[level] = SL;
230 sigmaN[level] = SN;
231 MaxP[level] = maxLP;
232 MaxN[level] = maxLN;
233 }
234
235
evaluate_params(wavelet_decomposition & wd,float * mean,float * meanN,float * sigma,float * sigmaN,float * MaxP,float * MaxN,bool multiThread)236 void evaluate_params(wavelet_decomposition &wd, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN, bool multiThread)
237 {
238 int maxlvl = wd.maxlevel();
239
240 for (int lvl = 0; lvl < maxlvl; lvl++) {
241 int Wlvl_L = wd.level_W(lvl);
242 int Hlvl_L = wd.level_H(lvl);
243
244 float **wl = wd.level_coeffs(lvl);
245
246 eval_level(wl, lvl, Wlvl_L, Hlvl_L, mean, meanN, sigma, sigmaN, MaxP, MaxN, multiThread);
247 }
248 }
249
250
local_contrast_wavelets(array2D<float> & Y,const LocalContrastParams::Region & params,double scale,bool multiThread)251 void local_contrast_wavelets(array2D<float> &Y, const LocalContrastParams::Region ¶ms, double scale, bool multiThread)
252 {
253 const int W = Y.width();
254 const int H = Y.height();
255
256 int wavelet_level = 7;
257 int dim = min(W, H);
258 while ((1 << wavelet_level) >= dim && wavelet_level > 1) {
259 --wavelet_level;
260 }
261 int skip = scale;
262 wavelet_decomposition wd(static_cast<float *>(Y), W, H, wavelet_level, 1, skip);
263
264 if (wd.memoryAllocationFailed) {
265 return;
266 }
267
268 const float contrast = params.contrast;
269 int maxlvl = wd.maxlevel();
270
271 if (contrast != 0) {
272 int W_L = wd.level_W(0);
273 int H_L = wd.level_H(0);
274 float *wl0 = wd.coeff0;
275
276 float maxh = 2.5f; //amplification contrast above mean
277 float maxl = 2.5f; //reduction contrast under mean
278 float multL = contrast * (maxl - 1.f) / 100.f + 1.f;
279 float multH = contrast * (maxh - 1.f) / 100.f + 1.f;
280 double avedbl = 0.0; // use double precision for large summations
281 float max0 = 0.f;
282 float min0 = FLT_MAX;
283
284 #ifdef _OPENMP
285 # pragma omp parallel for reduction(+:avedbl) if (multiThread)
286 #endif
287 for (int i = 0; i < W_L * H_L; i++) {
288 avedbl += wl0[i];
289 }
290
291 #ifdef _OPENMP
292 # pragma omp parallel if (multiThread)
293 #endif
294 {
295 float lminL = FLT_MAX;
296 float lmaxL = 0.f;
297
298 #ifdef _OPENMP
299 # pragma omp for
300 #endif
301 for (int i = 0; i < W_L * H_L; i++) {
302 lminL = min(lminL, wl0[i]);
303 lmaxL = max(lmaxL, wl0[i]);
304 }
305
306 #ifdef _OPENMP
307 # pragma omp critical
308 #endif
309 {
310 min0 = min(min0, lminL);
311 max0 = max(max0, lmaxL);
312 }
313 }
314
315 max0 /= 327.68f;
316 min0 /= 327.68f;
317 float ave = avedbl / double(W_L * H_L);
318 float av = ave / 327.68f;
319 float ah = (multH - 1.f) / (av - max0);
320 float bh = 1.f - max0 * ah;
321 float al = (multL - 1.f) / (av - min0);
322 float bl = 1.f - min0 * al;
323
324 if (max0 > 0.0) {
325 #ifdef _OPENMP
326 # pragma omp parallel for if (multiThread)
327 #endif
328 for (int i = 0; i < W_L * H_L; i++) {
329 if (wl0[i] < 32768.f) {
330 float prov;
331
332 if (wl0[i] > ave) {
333 float kh = ah * (wl0[i] / 327.68f) + bh;
334 prov = wl0[i];
335 wl0[i] = ave + kh * (wl0[i] - ave);
336 } else {
337 float kl = al * (wl0[i] / 327.68f) + bl;
338 prov = wl0[i];
339 wl0[i] = ave - kl * (ave - wl0[i]);
340 }
341
342 float diflc = wl0[i] - prov;
343 wl0[i] = prov + diflc;
344 }
345 }
346 }
347 }
348
349 float mean[10];
350 float meanN[10];
351 float sigma[10];
352 float sigmaN[10];
353 float MaxP[10];
354 float MaxN[10];
355 evaluate_params(wd, mean, meanN, sigma, sigmaN, MaxP, MaxN, multiThread);
356
357 WavOpacityCurveWL curve;
358 if (params.curve.empty() || params.curve[0] == FCT_Linear) {
359 LocalContrastParams::Region dflt;
360 curve.Set(dflt.curve);
361 } else {
362 curve.Set(params.curve);
363 }
364
365 for (int dir = 1; dir < 4; dir++) {
366 for (int level = 0; level < maxlvl; ++level) {
367 int W_L = wd.level_W(level);
368 int H_L = wd.level_H(level);
369 float **wl = wd.level_coeffs(level);
370
371 if (MaxP[level] > 0.f && mean[level] != 0.f && sigma[level] != 0.f) {
372 float insigma = 0.666f; //SD
373 float logmax = log(MaxP[level]); //log Max
374 float rapX = (mean[level] + sigma[level]) / MaxP[level]; //rapport between sD / max
375 float inx = log(insigma);
376 float iny = log(rapX);
377 float rap = inx / iny; //koef
378 float asig = 0.166f / sigma[level];
379 float bsig = 0.5f - asig * mean[level];
380 float amean = 0.5f / mean[level];
381
382 #ifdef _OPENMP
383 # pragma omp parallel for if (multiThread)
384 //schedule(dynamic, W_L * 16) if (multiThread)
385 #endif
386 for (int i = 0; i < W_L * H_L; i++) {
387 float absciss;
388 float &val = wl[dir][i];
389
390 if (std::isnan(val)) { // ALB -- TODO: this can happen in
391 // wavelet_decomposition, find out
392 // why
393 continue;
394 }
395
396 if (fabsf(val) >= (mean[level] + sigma[level])) { //for max
397 float valcour = xlogf(fabsf(val));
398 float valc = valcour - logmax;
399 float vald = valc * rap;
400 absciss = xexpf(vald);
401 } else if (fabsf(val) >= mean[level]) {
402 absciss = asig * fabsf(val) + bsig;
403 } else {
404 absciss = amean * fabsf(val);
405 }
406
407 float kc = curve[absciss * 500.f] - 0.5f;
408 float reduceeffect = kc <= 0.f ? 1.f : 1.5f;
409
410 float kinterm = 1.f + reduceeffect * kc;
411 kinterm = kinterm <= 0.f ? 0.01f : kinterm;
412
413 val *= kinterm;
414 }
415 }
416 }
417 }
418
419 wd.reconstruct(static_cast<float *>(Y), 1.f);
420 }
421
422 } // namespace
423
424
localContrast(Imagefloat * rgb)425 bool ImProcFunctions::localContrast(Imagefloat *rgb)
426 {
427 PlanarWhateverData<float> *editWhatever = nullptr;
428 EditUniqueID eid = pipetteBuffer ? pipetteBuffer->getEditID() : EUID_None;
429
430 if ((eid == EUID_LabMasks_H2 || eid == EUID_LabMasks_C2 || eid == EUID_LabMasks_L2) && pipetteBuffer->getDataProvider()->getCurrSubscriber()->getPipetteBufferType() == BT_SINGLEPLANE_FLOAT) {
431 editWhatever = pipetteBuffer->getSinglePlaneBuffer();
432 }
433
434 if (eid == EUID_LabMasks_DE2) {
435 if (getDeltaEColor(rgb, deltaE.x, deltaE.y, offset_x, offset_y, full_width, full_height, scale, deltaE.L, deltaE.C, deltaE.H)) {
436 deltaE.ok = true;
437 }
438 }
439
440 if (params->localContrast.enabled) {
441 rgb->setMode(Imagefloat::Mode::LAB, multiThread);
442
443 if (editWhatever) {
444 LabMasksEditID id = static_cast<LabMasksEditID>(int(eid) - EUID_LabMasks_H2);
445 fillPipetteLabMasks(rgb, editWhatever, id, multiThread);
446 }
447
448 int n = params->localContrast.regions.size();
449 int show_mask_idx = params->localContrast.showMask;
450 if (show_mask_idx >= n || (cur_pipeline != Pipeline::PREVIEW && cur_pipeline != Pipeline::OUTPUT)) {
451 show_mask_idx = -1;
452 }
453 std::vector<array2D<float>> mask(n);
454 if (!generateLabMasks(rgb, params->localContrast.labmasks, offset_x, offset_y, full_width, full_height, scale, multiThread, show_mask_idx, &mask, nullptr, cur_pipeline == Pipeline::NAVIGATOR ? plistener : nullptr)) {
455 return true; // show mask is active, nothing more to do
456 }
457
458 const int W = rgb->getWidth();
459 const int H = rgb->getHeight();
460
461 array2D<float> L(W, H, rgb->g.ptrs);
462
463 for (int i = 0; i < n; ++i) {
464 if (!params->localContrast.labmasks[i].enabled) {
465 continue;
466 }
467
468 auto &r = params->localContrast.regions[i];
469 local_contrast_wavelets(L, r, scale, multiThread);
470 const auto &blend = mask[i];
471 #ifdef _OPENMP
472 # pragma omp parallel for if (multiThread)
473 #endif
474 for (int y = 0; y < H; ++y) {
475 for (int x = 0; x < W; ++x) {
476 float l = rgb->g(y, x);
477 rgb->g(y, x) = intp(blend[y][x], L[y][x], l);
478 L[y][x] = rgb->g(y, x);
479 }
480 }
481 }
482 } else if (editWhatever) {
483 editWhatever->fill(0.f);
484 }
485
486 return false;
487 }
488
489 } // namespace rtengine
490