1 ////////////////////////////////////////////////////////////////
2 //
3 //
4 //
5 //
6 // code dated: December , 2014
7 //
8 // Ipwaveletcc is free software: you can redistribute it and/or modify
9 // it under the terms of the GNU General Public License as published by
10 // the Free Software Foundation, either version 3 of the License, or
11 // (at your option) any later version.
12 //
13 // This program is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with this program. If not, see <https://www.gnu.org/licenses/>.
20 // * 2014 Jacques Desmis <jdesmis@gmail.com>
21 // * 2014 Ingo Weyrich <heckflosse@i-weyrich.de>
22
23 //
24 ////////////////////////////////////////////////////////////////
25
26 #include <cassert>
27 #include <cmath>
28
29 #include "../rtgui/threadutils.h"
30
31 #include "array2D.h"
32 #include "color.h"
33 #include "curves.h"
34 #include "EdgePreservingDecomposition.h"
35 #include "iccstore.h"
36 #include "improcfun.h"
37 #include "labimage.h"
38 #include "LUT.h"
39 #include "median.h"
40 #include "opthelper.h"
41 #include "procparams.h"
42 #include "rt_math.h"
43 #include "rtengine.h"
44 #include "sleef.h"
45 #include "../rtgui/options.h"
46
47 #ifdef _OPENMP
48 #include <omp.h>
49 #endif
50
51 #include "cplx_wavelet_dec.h"
52
53 #define TS 64 // Tile size
54 #define offset 25 // shift between tiles
55 #define fTS ((TS/2+1)) // second dimension of Fourier tiles
56 #define blkrad 1 // radius of block averaging
57
58 #define epsilon 0.001f/(TS*TS) //tolerance
59
60 namespace rtengine
61 {
62
63 struct cont_params {
64 float mul[10];
65 int chrom;
66 int chro;
67 int contrast;
68 float th;
69 float thH;
70 float conres;
71 float conresH;
72 float chrores;
73 float hueres;
74 float sky;
75 float b_l, t_l, b_r, t_r;
76 float b_ly, t_ly, b_ry, t_ry;
77 float b_lsl, t_lsl, b_rsl, t_rsl;
78 float b_lhl, t_lhl, b_rhl, t_rhl;
79 float edg_low, edg_mean, edg_sd, edg_max;
80 float lev0s, lev0n, lev1s, lev1n, lev2s, lev2n, lev3s, lev3n;
81 float b_lpast, t_lpast, b_rpast, t_rpast;
82 float b_lsat, t_lsat, b_rsat, t_rsat;
83 int rad;
84 int val;
85 int til;
86 int numlevH, numlevS;
87 float mulC[9];
88 float mulopaRG[9];
89 float mulopaBY[9];
90 bool curv;
91 bool opaBY;
92 bool opaRG;
93 bool edgcurv;
94 bool diagcurv;
95 int CHmet;
96 int CHSLmet;
97 int EDmet;
98 bool HSmet;
99 bool avoi;
100 float strength;
101 int reinforce;
102 bool detectedge;
103 int backm;
104 float eddet;
105 float eddetthr;
106 float eddetthrHi;
107 bool link;
108 bool lip3;
109 bool tonemap;
110 bool diag;
111 int TMmeth;
112 float tmstrength;
113 float balan;
114 int ite;
115 int contmet;
116 bool opaW;
117 int BAmet;
118 bool bam;
119 float blhigh;
120 float grhigh;
121 float blmed;
122 float grmed;
123 float bllow;
124 float grlow;
125 bool cbena;
126 bool contena;
127 bool chromena;
128 bool edgeena;
129 bool resena;
130 bool finena;
131 bool toningena;
132 bool noiseena;
133 int maxilev;
134 float edgsens;
135 float edgampl;
136 int neigh;
137 bool lipp;
138 };
139
140 int wavNestedLevels = 1;
141
142
ip_wavelet(LabImage * lab,LabImage * dst,int kall,const procparams::WaveletParams & waparams,const WavCurve & wavCLVCcurve,const WavOpacityCurveRG & waOpacityCurveRG,const WavOpacityCurveBY & waOpacityCurveBY,const WavOpacityCurveW & waOpacityCurveW,const WavOpacityCurveWL & waOpacityCurveWL,LUTf & wavclCurve,int skip)143 void ImProcFunctions::ip_wavelet(LabImage * lab, LabImage * dst, int kall, const procparams::WaveletParams & waparams, const WavCurve & wavCLVCcurve, const WavOpacityCurveRG & waOpacityCurveRG, const WavOpacityCurveBY & waOpacityCurveBY, const WavOpacityCurveW & waOpacityCurveW, const WavOpacityCurveWL & waOpacityCurveWL, LUTf &wavclCurve, int skip)
144
145
146 {
147 #ifdef _DEBUG
148 // init variables to display Munsell corrections
149 MunsellDebugInfo* MunsDebugInfo = new MunsellDebugInfo();
150 #endif
151 TMatrix wiprof = ICCStore::getInstance()->workingSpaceInverseMatrix (params->icm.workingProfile);
152 const double wip[3][3] = {
153 {wiprof[0][0], wiprof[0][1], wiprof[0][2]},
154 {wiprof[1][0], wiprof[1][1], wiprof[1][2]},
155 {wiprof[2][0], wiprof[2][1], wiprof[2][2]}
156 };
157 const short int imheight = lab->H, imwidth = lab->W;
158 struct cont_params cp;
159 cp.avoi = params->wavelet.avoid;
160
161 if(params->wavelet.Medgreinf == "more") {
162 cp.reinforce = 1;
163 }
164
165 if(params->wavelet.Medgreinf == "none") {
166 cp.reinforce = 2;
167 }
168
169 if(params->wavelet.Medgreinf == "less") {
170 cp.reinforce = 3;
171 }
172
173 if(params->wavelet.NPmethod == "none") {
174 cp.lip3 = false;
175 }
176
177 if(params->wavelet.NPmethod == "low") {
178 cp.lip3 = true;
179 cp.neigh = 0;
180 }
181
182 if(params->wavelet.NPmethod == "high") {
183 cp.lip3 = true;
184 cp.neigh = 1;
185 }
186
187 cp.lipp = params->wavelet.lipst;
188 cp.diag = params->wavelet.tmr;
189 cp.balan = (float)params->wavelet.balance;
190 cp.ite = params->wavelet.iter;
191 cp.tonemap = false;
192 cp.bam = false;
193
194 if(params->wavelet.tmrs == 0) {
195 cp.tonemap = false;
196 } else {
197 cp.tonemap = true;
198 }
199
200 if(params->wavelet.TMmethod == "cont") {
201 cp.contmet = 1;
202 } else if(params->wavelet.TMmethod == "tm") {
203 cp.contmet = 2;
204 }
205
206 if(params->wavelet.BAmethod != "none") {
207 cp.bam = true;
208 }
209
210 if(params->wavelet.BAmethod == "sli") {
211 cp.BAmet = 1;
212 }
213
214 if(params->wavelet.BAmethod == "cur") {
215 cp.BAmet = 2;
216 }
217
218 cp.tmstrength = params->wavelet.tmrs;
219 //cp.tonemap = params->wavelet.tmr;
220 cp.contena = params->wavelet.expcontrast;
221 cp.chromena = params->wavelet.expchroma;
222 cp.edgeena = params->wavelet.expedge;
223 cp.resena = params->wavelet.expresid;
224 cp.finena = params->wavelet.expfinal;
225 cp.toningena = params->wavelet.exptoning;
226 cp.noiseena = params->wavelet.expnoise;
227
228 if(params->wavelet.Backmethod == "black") {
229 cp.backm = 0;
230 }
231
232 if(params->wavelet.Backmethod == "grey") {
233 cp.backm = 1;
234 }
235
236 if(params->wavelet.Backmethod == "resid") {
237 cp.backm = 2;
238 }
239
240 cp.link = params->wavelet.linkedg;
241 cp.eddet = (float) params->wavelet.edgedetect;
242 cp.eddetthr = (float) params->wavelet.edgedetectthr;
243 cp.eddetthrHi = (float) params->wavelet.edgedetectthr2;
244
245 cp.edgsens = 60.f;
246 cp.edgampl = 10.f;
247
248 if(cp.lipp) {
249 cp.edgsens = (float) params->wavelet.edgesensi;
250 cp.edgampl = (float) params->wavelet.edgeampli;
251 }
252
253 //int N = imheight * imwidth;
254 int maxmul = params->wavelet.thres;
255 cp.maxilev = maxmul;
256 static const float scales[10] = {1.f, 2.f, 4.f, 8.f, 16.f, 32.f, 64.f, 128.f, 256.f, 512.f};
257 float scaleskip[10];
258
259 for(int sc = 0; sc < 10; sc++) {
260 scaleskip[sc] = scales[sc] / skip;
261 }
262
263 float atten0 = 0.40f;
264 float atten123 = 0.90f;
265
266 //int DaubLen = settings->daubech ? 8 : 6;
267 int DaubLen;
268
269 if(params->wavelet.daubcoeffmethod == "2_") {
270 DaubLen = 4;
271 }
272
273 if(params->wavelet.daubcoeffmethod == "4_") {
274 DaubLen = 6;
275 }
276
277 if(params->wavelet.daubcoeffmethod == "6_") {
278 DaubLen = 8;
279 }
280
281 if(params->wavelet.daubcoeffmethod == "10_") {
282 DaubLen = 12;
283 }
284
285 if(params->wavelet.daubcoeffmethod == "14_") {
286 DaubLen = 16;
287 }
288
289 cp.CHSLmet = 1;
290 // if(params->wavelet.CHSLmethod=="SL") cp.CHSLmet=1;
291 // if(params->wavelet.CHSLmethod=="CU") cp.CHSLmet=2;
292 cp.EDmet = 1;
293
294 if(params->wavelet.EDmethod == "SL") {
295 cp.EDmet = 1;
296 }
297
298 if(params->wavelet.EDmethod == "CU") {
299 cp.EDmet = 2;
300 }
301
302 cp.cbena = params->wavelet.cbenab;
303 cp.blhigh = (float)params->wavelet.bluehigh;
304 cp.grhigh = (float)params->wavelet.greenhigh;
305 cp.blmed = (float)params->wavelet.bluemed;
306 cp.grmed = (float)params->wavelet.greenmed;
307 cp.bllow = (float)params->wavelet.bluelow;
308 cp.grlow = (float)params->wavelet.greenlow;
309 cp.curv = false;
310 cp.edgcurv = false;
311 cp.diagcurv = false;
312 cp.opaRG = false;
313 cp.opaBY = false;
314 cp.opaW = false;
315 cp.CHmet = 0;
316 cp.HSmet = false;
317
318 if(params->wavelet.CHmethod == "with") {
319 cp.CHmet = 1;
320 }
321
322 if(params->wavelet.CHmethod == "link") {
323 cp.CHmet = 2;
324 }
325
326 if(params->wavelet.HSmethod == "with") {
327 cp.HSmet = true;
328 }
329
330 cp.strength = min(1.f, max(0.f, ((float)params->wavelet.strength / 100.f)));
331
332 for(int m = 0; m < maxmul; m++) {
333 cp.mulC[m] = waparams.ch[m];
334 }
335
336 if(waOpacityCurveRG) {
337 cp.opaRG = true;
338 }
339
340 if(cp.opaRG) {
341 cp.mulopaRG[0] = 200.f * (waOpacityCurveRG[0] - 0.5f);
342 cp.mulopaRG[1] = 200.f * (waOpacityCurveRG[62] - 0.5f);
343 cp.mulopaRG[2] = 200.f * (waOpacityCurveRG[125] - 0.5f);
344 cp.mulopaRG[3] = 200.f * (waOpacityCurveRG[187] - 0.5f);
345 cp.mulopaRG[4] = 200.f * (waOpacityCurveRG[250] - 0.5f);
346 cp.mulopaRG[5] = 200.f * (waOpacityCurveRG[312] - 0.5f);
347 cp.mulopaRG[6] = 200.f * (waOpacityCurveRG[375] - 0.5f);
348 cp.mulopaRG[7] = 200.f * (waOpacityCurveRG[438] - 0.5f);
349 cp.mulopaRG[8] = 200.f * (waOpacityCurveRG[500] - 0.5f);
350 } else {
351 for(int level = 0; level < 9; level++) {
352 cp.mulopaRG[level] = 0.f;
353 }
354 }
355
356 if(waOpacityCurveBY) {
357 cp.opaBY = true;
358 }
359
360 if(cp.opaBY) {
361 cp.mulopaBY[0] = 200.f * (waOpacityCurveBY[0] - 0.5f);
362 cp.mulopaBY[1] = 200.f * (waOpacityCurveBY[62] - 0.5f);
363 cp.mulopaBY[2] = 200.f * (waOpacityCurveBY[125] - 0.5f);
364 cp.mulopaBY[3] = 200.f * (waOpacityCurveBY[187] - 0.5f);
365 cp.mulopaBY[4] = 200.f * (waOpacityCurveBY[250] - 0.5f);
366 cp.mulopaBY[5] = 200.f * (waOpacityCurveBY[312] - 0.5f);
367 cp.mulopaBY[6] = 200.f * (waOpacityCurveBY[375] - 0.5f);
368 cp.mulopaBY[7] = 200.f * (waOpacityCurveBY[438] - 0.5f);
369 cp.mulopaBY[8] = 200.f * (waOpacityCurveBY[500] - 0.5f);
370 } else {
371 for(int level = 0; level < 9; level++) {
372 cp.mulopaBY[level] = 0.f;
373 }
374 }
375
376 if(wavCLVCcurve) {
377 cp.edgcurv = true;
378 }
379
380 if(waOpacityCurveWL) {
381 cp.diagcurv = true;
382 }
383
384 for(int m = 0; m < maxmul; m++) {
385 cp.mul[m] = waparams.c[m];
386 }
387
388 cp.mul[9] = (float) waparams.sup;
389
390 for(int sc = 0; sc < 10; sc++) { //reduce strength if zoom < 100% for contrast
391 if(sc == 0) {
392 if(scaleskip[sc] < 1.f) {
393 cp.mul[sc] *= (atten0 * scaleskip[sc]);
394 }
395 } else {
396 if(scaleskip[sc] < 1.f) {
397 cp.mul[sc] *= (atten123 * scaleskip[sc]);
398 }
399 }
400 }
401
402 // if(settings->verbose) printf("Wav mul 0=%f 1=%f 2=%f 3=%f 4=%f 5=%f 6=%f 7=%f 8=%f 9=%f\n",cp.mul[0],cp.mul[1],cp.mul[2],cp.mul[3],cp.mul[4],cp.mul[5],cp.mul[6],cp.mul[7],cp.mul[8],cp.mul[9]);
403 for(int sc = 0; sc < 9; sc++) { //reduce strength if zoom < 100% for chroma and tuning
404 if(sc == 0) {
405 if(scaleskip[sc] < 1.f) {
406 cp.mulC[sc] *= (atten0 * scaleskip[sc]);
407 cp.mulopaRG[sc] *= (atten0 * scaleskip[sc]);
408 cp.mulopaBY[sc] *= (atten0 * scaleskip[sc]);
409 }
410 } else {
411 if(scaleskip[sc] < 1.f) {
412 cp.mulC[sc] *= (atten123 * scaleskip[sc]);
413 cp.mulopaRG[sc] *= (atten123 * scaleskip[sc]);
414 cp.mulopaBY[sc] *= (atten123 * scaleskip[sc]);
415 }
416 }
417 }
418
419 cp.chro = waparams.chro;
420 cp.chrom = waparams.chroma;
421 cp.contrast = waparams.contrast;
422 cp.rad = waparams.edgrad;
423 cp.val = waparams.edgval;
424 cp.til = waparams.edgthresh;
425
426 cp.conres = waparams.rescon;
427 cp.conresH = waparams.resconH;
428 cp.chrores = waparams.reschro;
429 //cp.hueres=waparams.reshue;
430 cp.hueres = 2.f;
431 cp.th = float(waparams.thr);
432 cp.thH = float(waparams.thrH);
433 cp.sky = waparams.sky;
434 //skin
435 cp.b_l = static_cast<float>(params->wavelet.hueskin.getBottomLeft()) / 100.0f;
436 cp.t_l = static_cast<float>(params->wavelet.hueskin.getTopLeft()) / 100.0f;
437 cp.b_r = static_cast<float>(params->wavelet.hueskin.getBottomRight()) / 100.0f;
438 cp.t_r = static_cast<float>(params->wavelet.hueskin.getTopRight()) / 100.0f;
439
440 cp.b_ly = static_cast<float>(params->wavelet.hueskin2.getBottomLeft()) / 100.0f;
441 cp.t_ly = static_cast<float>(params->wavelet.hueskin2.getTopLeft()) / 100.0f;
442 cp.b_ry = static_cast<float>(params->wavelet.hueskin2.getBottomRight()) / 100.0f;
443 cp.t_ry = static_cast<float>(params->wavelet.hueskin2.getTopRight()) / 100.0f;
444 cp.numlevH = params->wavelet.threshold;
445
446 //shadows
447 cp.b_lsl = static_cast<float>(params->wavelet.bllev.getBottomLeft());
448 cp.t_lsl = static_cast<float>(params->wavelet.bllev.getTopLeft());
449 cp.b_rsl = static_cast<float>(params->wavelet.bllev.getBottomRight());
450 cp.t_rsl = static_cast<float>(params->wavelet.bllev.getTopRight());
451 cp.numlevS = params->wavelet.threshold2;
452 int maxlevS = 9 - cp.numlevH;
453 cp.numlevS = MIN(cp.numlevS, maxlevS);
454 //printf("levHigh=%d levShad=%d\n",cp.numlevH,cp.numlevS);
455 //highlight
456 cp.b_lhl = static_cast<float>(params->wavelet.hllev.getBottomLeft());
457 cp.t_lhl = static_cast<float>(params->wavelet.hllev.getTopLeft());
458 cp.b_rhl = static_cast<float>(params->wavelet.hllev.getBottomRight());
459 cp.t_rhl = static_cast<float>(params->wavelet.hllev.getTopRight());
460 //printf("BL=%f TL=%f BR=%f TR=%f\n",cp.b_lhl,cp.t_lhl,cp.b_rhl,cp.t_rhl);
461 //pastel
462 cp.b_lpast = static_cast<float>(params->wavelet.pastlev.getBottomLeft());
463 cp.t_lpast = static_cast<float>(params->wavelet.pastlev.getTopLeft());
464 cp.b_rpast = static_cast<float>(params->wavelet.pastlev.getBottomRight());
465 cp.t_rpast = static_cast<float>(params->wavelet.pastlev.getTopRight());
466 //saturated
467 cp.b_lsat = static_cast<float>(params->wavelet.satlev.getBottomLeft());
468 cp.t_lsat = static_cast<float>(params->wavelet.satlev.getTopLeft());
469 cp.b_rsat = static_cast<float>(params->wavelet.satlev.getBottomRight());
470 cp.t_rsat = static_cast<float>(params->wavelet.satlev.getTopRight());
471 //edge local contrast
472 cp.edg_low = static_cast<float>(params->wavelet.edgcont.getBottomLeft());
473 cp.edg_mean = static_cast<float>(params->wavelet.edgcont.getTopLeft());
474 cp.edg_max = static_cast<float>(params->wavelet.edgcont.getBottomRight());
475 cp.edg_sd = static_cast<float>(params->wavelet.edgcont.getTopRight());
476 //level noise
477 cp.lev0s = static_cast<float>(params->wavelet.level0noise.getBottom());
478 cp.lev0n = static_cast<float>(params->wavelet.level0noise.getTop());
479 cp.lev1s = static_cast<float>(params->wavelet.level1noise.getBottom());
480 cp.lev1n = static_cast<float>(params->wavelet.level1noise.getTop());
481 cp.lev2s = static_cast<float>(params->wavelet.level2noise.getBottom());
482 cp.lev2n = static_cast<float>(params->wavelet.level2noise.getTop());
483 cp.lev3s = static_cast<float>(params->wavelet.level3noise.getBottom());
484 cp.lev3n = static_cast<float>(params->wavelet.level3noise.getTop());
485
486 cp.detectedge = params->wavelet.medianlev;
487 //printf("low=%f mean=%f sd=%f max=%f\n",cp.edg_low,cp.edg_mean,cp.edg_sd,cp.edg_max);
488 int minwin = min(imwidth, imheight);
489 int maxlevelcrop = 9;
490
491 if(cp.mul[9] != 0) {
492 maxlevelcrop = 10;
493 }
494
495 // adap maximum level wavelet to size of crop
496 if(minwin * skip < 1024) {
497 maxlevelcrop = 9; //sampling wavelet 512
498 }
499
500 if(minwin * skip < 512) {
501 maxlevelcrop = 8; //sampling wavelet 256
502 }
503
504 if(minwin * skip < 256) {
505 maxlevelcrop = 7; //sampling 128
506 }
507
508 if(minwin * skip < 128) {
509 maxlevelcrop = 6;
510 }
511
512 if(minwin < 64) {
513 maxlevelcrop = 5;
514 }
515
516 // printf("minwin=%d maxcrop=%d\n",minwin, maxlevelcrop);
517
518 int levwav = params->wavelet.thres;
519
520 if(levwav == 9 && cp.mul[9] != 0) {
521 levwav = 10;
522 }
523
524 levwav = min(maxlevelcrop, levwav);
525
526 // determine number of levels to process.
527 // for(levwav=min(maxlevelcrop,levwav);levwav>0;levwav--)
528 // if(cp.mul[levwav-1]!=0.f || cp.curv)
529 // if(cp.mul[levwav-1]!=0.f)
530 // break;
531 // I suppress this fonctionality ==> crash for level < 3
532 if(levwav < 1) {
533 return; // nothing to do
534 }
535
536 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
537 // begin tile processing of image
538
539 //output buffer
540 int realtile = 0;
541
542 if(params->wavelet.Tilesmethod == "big") {
543 realtile = 22;
544 }
545
546 if(params->wavelet.Tilesmethod == "lit") {
547 realtile = 12;
548 }
549
550 int tilesize = 128 * realtile;
551 int overlap = (int) tilesize * 0.125f;
552 int numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip;
553
554 if(params->wavelet.Tilesmethod == "full") {
555 kall = 0;
556 }
557
558 Tile_calc (tilesize, overlap, kall, imwidth, imheight, numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip);
559
560 const int numtiles = numtiles_W * numtiles_H;
561 LabImage * dsttmp;
562
563 if(numtiles == 1) {
564 dsttmp = dst;
565 } else {
566 dsttmp = new LabImage(imwidth, imheight);
567
568 for (int n = 0; n < 3 * imwidth * imheight; n++) {
569 dsttmp->data[n] = 0;
570 }
571 }
572
573 //now we have tile dimensions, overlaps
574 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
575 int minsizetile = min(tilewidth, tileheight);
576 int maxlev2 = 10;
577
578 if(minsizetile < 1024 && levwav == 10) {
579 maxlev2 = 9;
580 }
581
582 if(minsizetile < 512) {
583 maxlev2 = 8;
584 }
585
586 if(minsizetile < 256) {
587 maxlev2 = 7;
588 }
589
590 if(minsizetile < 128) {
591 maxlev2 = 6;
592 }
593
594 levwav = min(maxlev2, levwav);
595
596 //printf("levwav = %d\n",levwav);
597
598 #ifdef _OPENMP
599 int numthreads = 1;
600 int maxnumberofthreadsforwavelet = 0;
601
602 //reduce memory for big tile size
603 if(kall != 0) {
604 if(realtile <= 22) {
605 maxnumberofthreadsforwavelet = 2;
606 }
607
608 if(realtile <= 20) {
609 maxnumberofthreadsforwavelet = 3;
610 }
611
612 if(realtile <= 18) {
613 maxnumberofthreadsforwavelet = 4;
614 }
615
616 if(realtile <= 16) {
617 maxnumberofthreadsforwavelet = 6;
618 }
619
620 if(realtile <= 14) {
621 maxnumberofthreadsforwavelet = 8;
622 }
623
624 //printf("maxNRT=%d\n",maxnumberofthreadsforwavelet);
625 if((maxnumberofthreadsforwavelet == 6 || maxnumberofthreadsforwavelet == 8) && levwav == 10) {
626 maxnumberofthreadsforwavelet -= 2;
627 }
628
629 if(levwav <= 7 && maxnumberofthreadsforwavelet == 8) {
630 maxnumberofthreadsforwavelet = 0;
631 }
632 }
633
634 //printf("maxthre=%d\n",maxnumberofthreadsforwavelet);
635
636
637 // Calculate number of tiles. If less than omp_get_max_threads(), then limit num_threads to number of tiles
638 if( options.rgbDenoiseThreadLimit > 0) {
639 maxnumberofthreadsforwavelet = min(max(options.rgbDenoiseThreadLimit / 2, 1), maxnumberofthreadsforwavelet);
640 }
641
642 numthreads = MIN(numtiles, omp_get_max_threads());
643
644 if(maxnumberofthreadsforwavelet > 0) {
645 numthreads = MIN(numthreads, maxnumberofthreadsforwavelet);
646 }
647
648 #ifdef _OPENMP
649 wavNestedLevels = omp_get_max_threads() / numthreads;
650 bool oldNested = omp_get_nested();
651
652 if(wavNestedLevels < 2) {
653 wavNestedLevels = 1;
654 } else {
655 omp_set_nested(true);
656 }
657
658 if(maxnumberofthreadsforwavelet > 0)
659 while(wavNestedLevels * numthreads > maxnumberofthreadsforwavelet) {
660 wavNestedLevels--;
661 }
662
663 #endif
664
665 if(settings->verbose) {
666 printf("Ip Wavelet uses %d main thread(s) and up to %d nested thread(s) for each main thread\n", numthreads, wavNestedLevels);
667 }
668
669 #pragma omp parallel num_threads(numthreads)
670 #endif
671 {
672 float mean[10];
673 float meanN[10];
674 float sigma[10];
675 float sigmaN[10];
676 float MaxP[10];
677 float MaxN[10];
678
679 float** varhue = new float*[tileheight];
680
681 for (int i = 0; i < tileheight; i++) {
682 varhue[i] = new float[tilewidth];
683 }
684
685 float** varchro = new float*[tileheight];
686
687 for (int i = 0; i < tileheight; i++) {
688 varchro[i] = new float[tilewidth];
689 }
690
691 #ifdef _OPENMP
692 #pragma omp for schedule(dynamic) collapse(2)
693 #endif
694
695 for (int tiletop = 0; tiletop < imheight; tiletop += tileHskip) {
696 for (int tileleft = 0; tileleft < imwidth ; tileleft += tileWskip) {
697 int tileright = MIN(imwidth, tileleft + tilewidth);
698 int tilebottom = MIN(imheight, tiletop + tileheight);
699 int width = tileright - tileleft;
700 int height = tilebottom - tiletop;
701 LabImage * labco;
702 float **Lold = nullptr;
703 float *LoldBuffer = nullptr;
704
705 if(numtiles == 1) { // untiled processing => we can use output buffer for labco
706 labco = dst;
707
708 if(cp.avoi) { // we need a buffer to hold a copy of the L channel
709 Lold = new float*[tileheight];
710 LoldBuffer = new float[tilewidth * tileheight];
711 memcpy(LoldBuffer, lab->L[0], tilewidth * tileheight * sizeof(float));
712
713 for (int i = 0; i < tileheight; i++) {
714 Lold[i] = LoldBuffer + i * tilewidth;
715 }
716 }
717
718 } else {
719 labco = new LabImage(width, height);
720 Lold = lab->L;
721 }
722
723 #ifdef _OPENMP
724 #pragma omp parallel for num_threads(wavNestedLevels) if(wavNestedLevels>1)
725 #endif
726
727 for (int i = tiletop; i < tilebottom; i++) {
728 int i1 = i - tiletop;
729 int j;
730 #ifdef __SSE2__
731 __m128 c327d68v = _mm_set1_ps(327.68f);
732 __m128 av, bv, huev, chrov;
733
734 for (j = tileleft; j < tileright - 3; j += 4) {
735 int j1 = j - tileleft;
736 av = LVFU(lab->a[i][j]);
737 bv = LVFU(lab->b[i][j]);
738 huev = xatan2f(bv, av);
739 chrov = vsqrtf(SQRV(av) + SQRV(bv)) / c327d68v;
740 _mm_storeu_ps(&varhue[i1][j1], huev);
741 _mm_storeu_ps(&varchro[i1][j1], chrov);
742
743 if(labco != lab) {
744 _mm_storeu_ps(&(labco->L[i1][j1]), LVFU(lab->L[i][j]));
745 _mm_storeu_ps(&(labco->a[i1][j1]), av);
746 _mm_storeu_ps(&(labco->b[i1][j1]), bv);
747 }
748 }
749
750 #else
751 j = tileleft;
752 #endif
753
754 for (; j < tileright; j++) {
755 int j1 = j - tileleft;
756 float a = lab->a[i][j];
757 float b = lab->b[i][j];
758 varhue[i1][j1] = xatan2f(b, a);
759 varchro[i1][j1] = (sqrtf(a * a + b * b)) / 327.68f;
760
761 if(labco != lab) {
762 labco->L[i1][j1] = lab->L[i][j];
763 labco->a[i1][j1] = a;
764 labco->b[i1][j1] = b;
765 }
766 }
767 }
768
769 //to avoid artifacts in blue sky
770 if(params->wavelet.median) {
771 float** tmL;
772 int wid = labco->W;
773 int hei = labco->H;
774 int borderL = 1;
775 tmL = new float*[hei];
776
777 for (int i = 0; i < hei; i++) {
778 tmL[i] = new float[wid];
779 }
780
781 for(int i = borderL; i < hei - borderL; i++ ) {
782 for(int j = borderL; j < wid - borderL; j++) {
783 tmL[i][j] = labco->L[i][j];
784 }
785 }
786
787 #ifdef _OPENMP
788 #pragma omp parallel for num_threads(wavNestedLevels) if(wavNestedLevels>1)
789 #endif
790
791 for (int i = 1; i < hei - 1; i++) {
792 for (int j = 1; j < wid - 1; j++) {
793 if((varhue[i][j] < -1.3f && varhue[i][j] > - 2.5f) && (varchro[i][j] > 15.f && varchro[i][j] < 55.f) && labco->L[i][j] > 6000.f) { //blue sky + med3x3 ==> after for more effect use denoise
794 tmL[i][j] = median(labco->L[i][j] , labco->L[i - 1][j], labco->L[i + 1][j] , labco->L[i][j + 1], labco->L[i][j - 1], labco->L[i - 1][j - 1], labco->L[i - 1][j + 1], labco->L[i + 1][j - 1], labco->L[i + 1][j + 1]); //3x3
795 }
796 }
797 }
798
799 for(int i = borderL; i < hei - borderL; i++ ) {
800 for(int j = borderL; j < wid - borderL; j++) {
801 labco->L[i][j] = tmL[i][j];
802 }
803 }
804
805 for (int i = 0; i < hei; i++) {
806 delete [] tmL[i];
807 }
808
809 delete [] tmL;
810 // end blue sky
811 }
812
813 if(numtiles == 1) {
814 // reduce the varhue array to get faster access in following processing and reduce peak memory usage
815 float temphue[(tilewidth + 1) / 2] ALIGNED64;
816
817 for (int i = 0; i < (tileheight + 1) / 2; i++) {
818 for (int j = 0; j < (tilewidth + 1) / 2; j++) {
819 temphue[j] = varhue[i * 2][j * 2];
820 }
821
822 delete [] varhue[i];
823 varhue[i] = new float[(tilewidth + 1) / 2];
824 memcpy(varhue[i], temphue, ((tilewidth + 1) / 2) * sizeof(float));
825 }
826
827 for(int i = (tileheight + 1) / 2; i < tileheight; i++) {
828 delete [] varhue[i];
829 varhue[i] = nullptr;
830 }
831 } else { // reduce the varhue array to get faster access in following processing
832 for (int i = 0; i < (tileheight + 1) / 2; i++) {
833 for (int j = 0; j < (tilewidth + 1) / 2; j++) {
834 varhue[i][j] = varhue[i * 2][j * 2];
835 }
836 }
837 }
838
839 int datalen = labco->W * labco->H;
840
841 int levwavL = levwav;
842 bool ref0 = false;
843
844 if((cp.lev0s > 0.f || cp.lev1s > 0.f || cp.lev2s > 0.f || cp.lev3s > 0.f) && cp.noiseena) {
845 ref0 = true;
846 }
847
848 // printf("LevwavL before: %d\n",levwavL);
849 if(cp.contrast == 0.f && !cp.tonemap && cp.conres == 0.f && cp.conresH == 0.f && cp.val == 0 && !ref0 && params->wavelet.CLmethod == "all") { // no processing of residual L or edge=> we probably can reduce the number of levels
850 while(levwavL > 0 && cp.mul[levwavL - 1] == 0.f) { // cp.mul[level] == 0.f means no changes to level
851 levwavL--;
852 }
853 }
854
855 // printf("LevwavL after: %d\n",levwavL);
856 // if(cp.noiseena){
857 if(levwavL < 4 ) {
858 levwavL = 4; //to allow edge => I always allocate 3 (4) levels..because if user select wavelet it is to do something !!
859 }
860
861 // }
862 // else {
863 // if(levwavL < 3) levwavL=3;//to allow edge => I always allocate 3 (4) levels..because if user select wavelet it is to do something !!
864 // }
865 if(levwavL > 0) {
866 wavelet_decomposition* Ldecomp = new wavelet_decomposition (labco->data, labco->W, labco->H, levwavL, 1, skip, max(1, wavNestedLevels), DaubLen );
867
868 if(!Ldecomp->memoryAllocationFailed) {
869
870 float madL[8][3];
871 #ifdef _OPENMP
872 #pragma omp parallel for schedule(dynamic) collapse(2) num_threads(wavNestedLevels) if(wavNestedLevels>1)
873 #endif
874
875 for (int lvl = 0; lvl < 4; lvl++) {
876 for (int dir = 1; dir < 4; dir++) {
877 int Wlvl_L = Ldecomp->level_W(lvl);
878 int Hlvl_L = Ldecomp->level_H(lvl);
879
880 float ** WavCoeffs_L = Ldecomp->level_coeffs(lvl);
881
882 madL[lvl][dir - 1] = SQR(Mad(WavCoeffs_L[dir], Wlvl_L * Hlvl_L));
883 }
884 }
885
886 bool ref = false;
887
888 if((cp.lev0s > 0.f || cp.lev1s > 0.f || cp.lev2s > 0.f || cp.lev3s > 0.f) && cp.noiseena) {
889 ref = true;
890 }
891
892 bool contr = false;
893
894 for(int f = 0; f < levwavL; f++) {
895 if(cp.mul[f] != 0.f) {
896 contr = true;
897 }
898 }
899
900 if(cp.val > 0 || ref || contr) {//edge
901 Evaluate2(*Ldecomp, mean, meanN, sigma, sigmaN, MaxP, MaxN);
902 }
903
904 //init for edge and denoise
905 float vari[4];
906
907 vari[0] = 8.f * SQR((cp.lev0n / 125.0) * (1.0 + cp.lev0n / 25.0));
908 vari[1] = 8.f * SQR((cp.lev1n / 125.0) * (1.0 + cp.lev1n / 25.0));
909 vari[2] = 8.f * SQR((cp.lev2n / 125.0) * (1.0 + cp.lev2n / 25.0));
910 vari[3] = 8.f * SQR((cp.lev3n / 125.0) * (1.0 + cp.lev3n / 25.0));
911
912 if((cp.lev0n > 0.1f || cp.lev1n > 0.1f || cp.lev2n > 0.1f || cp.lev3n > 0.1f) && cp.noiseena) {
913 int edge = 1;
914 vari[0] = max(0.0001f, vari[0]);
915 vari[1] = max(0.0001f, vari[1]);
916 vari[2] = max(0.0001f, vari[2]);
917 vari[3] = max(0.0001f, vari[3]);
918 float* noisevarlum = nullptr; // we need a dummy to pass it to WaveletDenoiseAllL
919
920 WaveletDenoiseAllL(*Ldecomp, noisevarlum, madL, vari, edge);
921 }
922
923 //Flat curve for Contrast=f(H) in levels
924 FlatCurve* ChCurve = new FlatCurve(params->wavelet.Chcurve); //curve C=f(H)
925 bool Chutili = false;
926
927 if (!ChCurve || ChCurve->isIdentity()) {
928 if (ChCurve) {
929 delete ChCurve;
930 ChCurve = nullptr;
931 }
932 } else {
933 Chutili = true;
934 }
935
936
937 WaveletcontAllL(labco, varhue, varchro, *Ldecomp, cp, skip, mean, sigma, MaxP, MaxN, wavCLVCcurve, waOpacityCurveW, ChCurve, Chutili);
938
939 if(cp.val > 0 || ref || contr || cp.diagcurv) {//edge
940 Evaluate2(*Ldecomp, mean, meanN, sigma, sigmaN, MaxP, MaxN);
941 }
942
943 WaveletcontAllLfinal(*Ldecomp, cp, mean, sigma, MaxP, waOpacityCurveWL);
944 //Evaluate2(*Ldecomp, cp, ind, mean, meanN, sigma, sigmaN, MaxP, MaxN, madL);
945
946 Ldecomp->reconstruct(labco->data, cp.strength);
947 }
948
949 delete Ldecomp;
950 }
951
952 //Flat curve for H=f(H) in residual image
953 FlatCurve* hhCurve = new FlatCurve(params->wavelet.hhcurve); //curve H=f(H)
954 bool hhutili = false;
955
956 if (!hhCurve || hhCurve->isIdentity()) {
957 if (hhCurve) {
958 delete hhCurve;
959 hhCurve = nullptr;
960 }
961 } else {
962 hhutili = true;
963 }
964
965
966 if(!hhutili) {//always a or b
967 int levwava = levwav;
968
969 // printf("Levwava before: %d\n",levwava);
970 if(cp.chrores == 0.f && params->wavelet.CLmethod == "all" && !cp.cbena) { // no processing of residual ab => we probably can reduce the number of levels
971 while(levwava > 0 && !cp.diag && (((cp.CHmet == 2 && (cp.chro == 0.f || cp.mul[levwava - 1] == 0.f )) || (cp.CHmet != 2 && (levwava == 10 || (!cp.curv || cp.mulC[levwava - 1] == 0.f))))) && (!cp.opaRG || levwava == 10 || (cp.opaRG && cp.mulopaRG[levwava - 1] == 0.f)) && ((levwava == 10 || (cp.CHSLmet == 1 && cp.mulC[levwava - 1] == 0.f)))) {
972 levwava--;
973 }
974 }
975
976 //printf("Levwava after: %d\n",levwava);
977 if(levwava > 0) {
978 wavelet_decomposition* adecomp = new wavelet_decomposition (labco->data + datalen, labco->W, labco->H, levwava, 1, skip, max(1, wavNestedLevels), DaubLen );
979
980 if(!adecomp->memoryAllocationFailed) {
981 WaveletcontAllAB(labco, varhue, varchro, *adecomp, waOpacityCurveW, cp, true);
982 adecomp->reconstruct(labco->data + datalen, cp.strength);
983 }
984
985 delete adecomp;
986 }
987
988 int levwavb = levwav;
989
990 //printf("Levwavb before: %d\n",levwavb);
991 if(cp.chrores == 0.f && params->wavelet.CLmethod == "all" && !cp.cbena) { // no processing of residual ab => we probably can reduce the number of levels
992 while(levwavb > 0 && !cp.diag && (((cp.CHmet == 2 && (cp.chro == 0.f || cp.mul[levwavb - 1] == 0.f )) || (cp.CHmet != 2 && (levwavb == 10 || (!cp.curv || cp.mulC[levwavb - 1] == 0.f))))) && (!cp.opaBY || levwavb == 10 || (cp.opaBY && cp.mulopaBY[levwavb - 1] == 0.f)) && ((levwavb == 10 || (cp.CHSLmet == 1 && cp.mulC[levwavb - 1] == 0.f)))) {
993 levwavb--;
994 }
995 }
996
997 // printf("Levwavb after: %d\n",levwavb);
998 if(levwavb > 0) {
999 wavelet_decomposition* bdecomp = new wavelet_decomposition (labco->data + 2 * datalen, labco->W, labco->H, levwavb, 1, skip, max(1, wavNestedLevels), DaubLen );
1000
1001 if(!bdecomp->memoryAllocationFailed) {
1002 WaveletcontAllAB(labco, varhue, varchro, *bdecomp, waOpacityCurveW, cp, false);
1003 bdecomp->reconstruct(labco->data + 2 * datalen, cp.strength);
1004 }
1005
1006 delete bdecomp;
1007 }
1008 } else {// a and b
1009 int levwavab = levwav;
1010
1011 // printf("Levwavab before: %d\n",levwavab);
1012 if(cp.chrores == 0.f && !hhutili && params->wavelet.CLmethod == "all") { // no processing of residual ab => we probably can reduce the number of levels
1013 while(levwavab > 0 && (((cp.CHmet == 2 && (cp.chro == 0.f || cp.mul[levwavab - 1] == 0.f )) || (cp.CHmet != 2 && (levwavab == 10 || (!cp.curv || cp.mulC[levwavab - 1] == 0.f))))) && (!cp.opaRG || levwavab == 10 || (cp.opaRG && cp.mulopaRG[levwavab - 1] == 0.f)) && ((levwavab == 10 || (cp.CHSLmet == 1 && cp.mulC[levwavab - 1] == 0.f)))) {
1014 levwavab--;
1015 }
1016 }
1017
1018 // printf("Levwavab after: %d\n",levwavab);
1019 if(levwavab > 0) {
1020 wavelet_decomposition* adecomp = new wavelet_decomposition (labco->data + datalen, labco->W, labco->H, levwavab, 1, skip, max(1, wavNestedLevels), DaubLen );
1021 wavelet_decomposition* bdecomp = new wavelet_decomposition (labco->data + 2 * datalen, labco->W, labco->H, levwavab, 1, skip, max(1, wavNestedLevels), DaubLen );
1022
1023 if(!adecomp->memoryAllocationFailed && !bdecomp->memoryAllocationFailed) {
1024 WaveletcontAllAB(labco, varhue, varchro, *adecomp, waOpacityCurveW, cp, true);
1025 WaveletcontAllAB(labco, varhue, varchro, *bdecomp, waOpacityCurveW, cp, false);
1026 WaveletAandBAllAB(*adecomp, *bdecomp, cp, hhCurve, hhutili );
1027
1028 adecomp->reconstruct(labco->data + datalen, cp.strength);
1029 bdecomp->reconstruct(labco->data + 2 * datalen, cp.strength);
1030
1031 }
1032
1033 delete adecomp;
1034 delete bdecomp;
1035 }
1036 }
1037
1038 if (hhCurve) {
1039 delete hhCurve;
1040 }
1041
1042 if(numtiles > 1 || (numtiles == 1 /*&& cp.avoi*/)) {//in all case since I add contrast curve
1043 //calculate mask for feathering output tile overlaps
1044 float Vmask[height + overlap] ALIGNED16;
1045 float Hmask[width + overlap] ALIGNED16;
1046
1047 if(numtiles > 1) {
1048 for (int i = 0; i < height; i++) {
1049 Vmask[i] = 1;
1050 }
1051
1052 for (int j = 0; j < width; j++) {
1053 Hmask[j] = 1;
1054 }
1055
1056 for (int i = 0; i < overlap; i++) {
1057 float mask = SQR(sin((rtengine::RT_PI * i) / (2 * overlap)));
1058
1059 if (tiletop > 0) {
1060 Vmask[i] = mask;
1061 }
1062
1063 if (tilebottom < imheight) {
1064 Vmask[height - i] = mask;
1065 }
1066
1067 if (tileleft > 0) {
1068 Hmask[i] = mask;
1069 }
1070
1071 if (tileright < imwidth) {
1072 Hmask[width - i] = mask;
1073 }
1074 }
1075 }
1076
1077 bool highlight = params->toneCurve.hrenabled;
1078
1079 #ifdef _OPENMP
1080 #pragma omp parallel for schedule(dynamic,16) num_threads(wavNestedLevels) if(wavNestedLevels>1)
1081 #endif
1082
1083 for (int i = tiletop; i < tilebottom; i++) {
1084 int i1 = i - tiletop;
1085 float L, a, b;
1086 #ifdef __SSE2__
1087 int rowWidth = tileright - tileleft;
1088 float atan2Buffer[rowWidth] ALIGNED64;
1089 float chprovBuffer[rowWidth] ALIGNED64;
1090 float xBuffer[rowWidth] ALIGNED64;
1091 float yBuffer[rowWidth] ALIGNED64;
1092
1093 if(cp.avoi) {
1094 int col;
1095 __m128 av, bv;
1096 __m128 cv, yv, xv;
1097 __m128 zerov = _mm_setzero_ps();
1098 __m128 onev = _mm_set1_ps(1.f);
1099 __m128 c327d68v = _mm_set1_ps(327.68f);
1100 vmask xyMask;
1101
1102 for(col = 0; col < rowWidth - 3; col += 4) {
1103 av = LVFU(labco->a[i1][col]);
1104 bv = LVFU(labco->b[i1][col]);
1105 STVF(atan2Buffer[col], xatan2f(bv, av));
1106
1107 cv = vsqrtf(SQRV(av) + SQRV(bv));
1108 yv = av / cv;
1109 xv = bv / cv;
1110 xyMask = vmaskf_eq(zerov, cv);
1111 yv = vself(xyMask, onev, yv);
1112 xv = vself(xyMask, zerov, xv);
1113 STVF(yBuffer[col], yv);
1114 STVF(xBuffer[col], xv);
1115 STVF(chprovBuffer[col], cv / c327d68v);
1116
1117 }
1118
1119 for(; col < rowWidth; col++) {
1120 float a = labco->a[i1][col];
1121 float b = labco->b[i1][col];
1122 atan2Buffer[col] = xatan2f(b, a);
1123 float Chprov1 = sqrtf(SQR(a) + SQR(b));
1124 yBuffer[col] = (Chprov1 == 0.f) ? 1.f : a / Chprov1;
1125 xBuffer[col] = (Chprov1 == 0.f) ? 0.f : b / Chprov1;
1126 chprovBuffer[col] = Chprov1 / 327.68;
1127 }
1128 }
1129
1130 #endif
1131
1132 for (int j = tileleft; j < tileright; j++) {
1133 int j1 = j - tileleft;
1134
1135 if(cp.avoi) { //Gamut and Munsell
1136 #ifdef __SSE2__
1137 float HH = atan2Buffer[j1];
1138 float Chprov1 = chprovBuffer[j1];
1139 float2 sincosv;
1140 sincosv.y = yBuffer[j1];
1141 sincosv.x = xBuffer[j1];
1142 #else
1143 a = labco->a[i1][j1];
1144 b = labco->b[i1][j1];
1145 float HH = xatan2f(b, a);
1146 float Chprov1 = sqrtf(SQR(a) + SQR(b));
1147 float2 sincosv;
1148 sincosv.y = (Chprov1 == 0.0f) ? 1.f : a / (Chprov1);
1149 sincosv.x = (Chprov1 == 0.0f) ? 0.f : b / (Chprov1);
1150 Chprov1 /= 327.68f;
1151 #endif
1152 L = labco->L[i1][j1];
1153 const float Lin = labco->L[i1][j1];
1154
1155 if(wavclCurve && cp.finena) {
1156 labco->L[i1][j1] = (0.5f * Lin + 1.5f * wavclCurve[Lin]) / 2.f; //apply contrast curve
1157 }
1158
1159 L = labco->L[i1][j1];
1160
1161 float Lprov1 = L / 327.68f;
1162 float Lprov2 = Lold[i][j] / 327.68f;
1163 float memChprov = varchro[i1][j1];
1164 float R, G, B;
1165 #ifdef _DEBUG
1166 bool neg = false;
1167 bool more_rgb = false;
1168 Color::gamutLchonly(HH, sincosv, Lprov1, Chprov1, R, G, B, wip, highlight, 0.15f, 0.96f, neg, more_rgb);
1169 #else
1170 Color::gamutLchonly(HH, sincosv, Lprov1, Chprov1, R, G, B, wip, highlight, 0.15f, 0.96f);
1171 #endif
1172 L = Lprov1 * 327.68f;
1173
1174 a = 327.68f * Chprov1 * sincosv.y; //gamut
1175 b = 327.68f * Chprov1 * sincosv.x; //gamut
1176 float correctionHue = 0.0f; // Munsell's correction
1177 float correctlum = 0.0f;
1178 Lprov1 = L / 327.68f;
1179 float Chprov = sqrtf(SQR(a) + SQR(b)) / 327.68f;
1180 #ifdef _DEBUG
1181 Color::AllMunsellLch(true, Lprov1, Lprov2, HH, Chprov, memChprov, correctionHue, correctlum, MunsDebugInfo);
1182 #else
1183 Color::AllMunsellLch(true, Lprov1, Lprov2, HH, Chprov, memChprov, correctionHue, correctlum);
1184 #endif
1185
1186 if(correctionHue != 0.f || correctlum != 0.f) { // only calculate sin and cos if HH changed
1187 if(fabs(correctionHue) < 0.015f) {
1188 HH += correctlum; // correct only if correct Munsell chroma very little.
1189 }
1190
1191 sincosv = xsincosf(HH + correctionHue);
1192 }
1193
1194 a = 327.68f * Chprov * sincosv.y; // apply Munsell
1195 b = 327.68f * Chprov * sincosv.x; //aply Munsell
1196 } else {//general case
1197 L = labco->L[i1][j1];
1198 const float Lin = std::max(0.f, L);
1199
1200 if(wavclCurve && cp.finena) {
1201 labco->L[i1][j1] = (0.5f * Lin + 1.5f * wavclCurve[Lin]) / 2.f; //apply contrast curve
1202 }
1203
1204 L = labco->L[i1][j1];
1205 a = labco->a[i1][j1];
1206 b = labco->b[i1][j1];
1207 }
1208
1209 if(numtiles > 1) {
1210 float factor = Vmask[i1] * Hmask[j1];
1211 dsttmp->L[i][j] += factor * L;
1212 dsttmp->a[i][j] += factor * a;
1213 dsttmp->b[i][j] += factor * b;
1214 } else {
1215 dsttmp->L[i][j] = L;
1216 dsttmp->a[i][j] = a;
1217 dsttmp->b[i][j] = b;
1218
1219 }
1220 }
1221 }
1222 }
1223
1224 if(LoldBuffer != nullptr) {
1225 delete [] LoldBuffer;
1226 delete [] Lold;
1227 }
1228
1229 if(numtiles > 1) {
1230 delete labco;
1231 }
1232 }
1233 }
1234
1235 for (int i = 0; i < tileheight; i++)
1236 if(varhue[i] != nullptr) {
1237 delete [] varhue[i];
1238 }
1239
1240 delete [] varhue;
1241
1242 for (int i = 0; i < tileheight; i++) {
1243 delete [] varchro[i];
1244 }
1245
1246 delete [] varchro;
1247
1248 }
1249 #ifdef _OPENMP
1250 omp_set_nested(oldNested);
1251 #endif
1252
1253 if(numtiles != 1) {
1254 dst->CopyFrom(dsttmp);
1255 delete dsttmp;
1256 }
1257
1258 #ifdef _DEBUG
1259 delete MunsDebugInfo;
1260 #endif
1261
1262 }
1263
1264 #undef TS
1265 #undef fTS
1266 #undef offset
1267 #undef epsilon
1268
Aver(float * RESTRICT DataList,int datalen,float & averagePlus,float & averageNeg,float & max,float & min)1269 void ImProcFunctions::Aver( float * RESTRICT DataList, int datalen, float &averagePlus, float &averageNeg, float &max, float &min)
1270 {
1271
1272 //find absolute mean
1273 int countP = 0, countN = 0;
1274 double averaP = 0.0, averaN = 0.0; // use double precision for large summations
1275
1276 float thres = 5.f;//different fom zero to take into account only data large enough
1277 max = 0.f;
1278 min = 0.f;
1279 #ifdef _OPENMP
1280 #pragma omp parallel num_threads(wavNestedLevels) if(wavNestedLevels>1)
1281 #endif
1282 {
1283 float lmax = 0.f, lmin = 0.f;
1284 #ifdef _OPENMP
1285 #pragma omp for reduction(+:averaP,averaN,countP,countN) nowait
1286 #endif
1287
1288 for(int i = 0; i < datalen; i++) {
1289 if(DataList[i] >= thres) {
1290 averaP += DataList[i];
1291
1292 if(DataList[i] > lmax) {
1293 lmax = DataList[i];
1294 }
1295
1296 countP++;
1297 } else if(DataList[i] < -thres) {
1298 averaN += DataList[i];
1299
1300 if(DataList[i] < lmin) {
1301 lmin = DataList[i];
1302 }
1303
1304 countN++;
1305 }
1306 }
1307
1308 #ifdef _OPENMP
1309 #pragma omp critical
1310 #endif
1311 {
1312 max = max > lmax ? max : lmax;
1313 min = min < lmin ? min : lmin;
1314 }
1315 }
1316
1317 if(countP > 0) {
1318 averagePlus = averaP / countP;
1319 } else {
1320 averagePlus = 0;
1321 }
1322
1323 if(countN > 0) {
1324 averageNeg = averaN / countN;
1325 } else {
1326 averageNeg = 0;
1327 }
1328
1329 }
1330
1331
Sigma(float * RESTRICT DataList,int datalen,float averagePlus,float averageNeg,float & sigmaPlus,float & sigmaNeg)1332 void ImProcFunctions::Sigma( float * RESTRICT DataList, int datalen, float averagePlus, float averageNeg, float &sigmaPlus, float &sigmaNeg)
1333 {
1334 int countP = 0, countN = 0;
1335 double variP = 0.0, variN = 0.0; // use double precision for large summations
1336 float thres = 5.f;//different fom zero to take into account only data large enough
1337
1338 #ifdef _OPENMP
1339 #pragma omp parallel for reduction(+:variP,variN,countP,countN) num_threads(wavNestedLevels) if(wavNestedLevels>1)
1340 #endif
1341
1342 for(int i = 0; i < datalen; i++) {
1343 if(DataList[i] >= thres) {
1344 variP += SQR(DataList[i] - averagePlus);
1345 countP++;
1346 } else if(DataList[i] <= -thres) {
1347 variN += SQR(DataList[i] - averageNeg);
1348 countN++;
1349 }
1350 }
1351
1352 if(countP > 0) {
1353 sigmaPlus = sqrt(variP / countP);
1354 } else {
1355 sigmaPlus = 0;
1356 }
1357
1358 if(countN > 0) {
1359 sigmaNeg = sqrt(variN / countN);
1360 } else {
1361 sigmaNeg = 0;
1362 }
1363
1364 }
1365
Evaluate2(wavelet_decomposition & WaveletCoeffs_L,float * mean,float * meanN,float * sigma,float * sigmaN,float * MaxP,float * MaxN)1366 void ImProcFunctions::Evaluate2(wavelet_decomposition &WaveletCoeffs_L,
1367 float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN)
1368 {
1369 //StopWatch Stop1("Evaluate2");
1370 int maxlvl = WaveletCoeffs_L.maxlevel();
1371
1372 for (int lvl = 0; lvl < maxlvl; lvl++) {
1373
1374 int Wlvl_L = WaveletCoeffs_L.level_W(lvl);
1375 int Hlvl_L = WaveletCoeffs_L.level_H(lvl);
1376
1377 float ** WavCoeffs_L = WaveletCoeffs_L.level_coeffs(lvl);
1378
1379 Eval2 (WavCoeffs_L, lvl, Wlvl_L, Hlvl_L, mean, meanN, sigma, sigmaN, MaxP, MaxN);
1380 }
1381
1382 }
Eval2(float ** WavCoeffs_L,int level,int W_L,int H_L,float * mean,float * meanN,float * sigma,float * sigmaN,float * MaxP,float * MaxN)1383 void ImProcFunctions::Eval2 (float ** WavCoeffs_L, int level,
1384 int W_L, int H_L, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN)
1385 {
1386
1387 float avLP[4], avLN[4];
1388 float maxL[4], minL[4];
1389 float sigP[4], sigN[4];
1390 float AvL, AvN, SL, SN, maxLP, maxLN;
1391
1392 for (int dir = 1; dir < 4; dir++) {
1393 Aver(WavCoeffs_L[dir], W_L * H_L, avLP[dir], avLN[dir], maxL[dir], minL[dir]);
1394 Sigma(WavCoeffs_L[dir], W_L * H_L, avLP[dir], avLN[dir], sigP[dir], sigN[dir]);
1395 }
1396
1397 AvL = 0.f;
1398 AvN = 0.f;
1399 SL = 0.f;
1400 SN = 0.f;
1401 maxLP = 0.f;
1402 maxLN = 0.f;
1403
1404 for (int dir = 1; dir < 4; dir++) {
1405 AvL += avLP[dir];
1406 AvN += avLN[dir];
1407 SL += sigP[dir];
1408 SN += sigN[dir];
1409 maxLP += maxL[dir];
1410 maxLN += minL[dir];
1411 }
1412
1413 AvL /= 3;
1414 AvN /= 3;
1415 SL /= 3;
1416 SN /= 3;
1417 maxLP /= 3;
1418 maxLN /= 3;
1419
1420 mean[level] = AvL;
1421 meanN[level] = AvN;
1422 sigma[level] = SL;
1423 sigmaN[level] = SN;
1424 MaxP[level] = maxLP;
1425 MaxN[level] = maxLN;
1426 }
1427
CompressDR(float * Source,int W_L,int H_L,float Compression,float DetailBoost)1428 void ImProcFunctions::CompressDR(float *Source, int W_L, int H_L, float Compression, float DetailBoost)
1429 {
1430 const int n = W_L * H_L;
1431
1432 float exponent;
1433
1434 if(DetailBoost > 0.f && DetailBoost < 0.05f ) {
1435 float betemp = expf(-(2.f - DetailBoost + 0.694f)) - 1.f; //0.694 = log(2)
1436 exponent = 1.2f * xlogf( -betemp);
1437 exponent /= 20.f;
1438 } else if(DetailBoost >= 0.05f && DetailBoost < 0.25f ) {
1439 float betemp = expf(-(2.f - DetailBoost + 0.694f)) - 1.f; //0.694 = log(2)
1440 exponent = 1.2f * xlogf( -betemp);
1441 exponent /= (-75.f * DetailBoost + 23.75f);
1442 } else if(DetailBoost >= 0.25f) {
1443 float betemp = expf(-(2.f - DetailBoost + 0.694f)) - 1.f; //0.694 = log(2)
1444 exponent = 1.2f * xlogf( -betemp);
1445 exponent /= (-2.f * DetailBoost + 5.5f);
1446 } else {
1447 exponent = (Compression - 1.0f) / 20.f;
1448 }
1449
1450 exponent += 1.f;
1451
1452 // now calculate Source = pow(Source, exponent)
1453 #ifdef __SSE2__
1454 #ifdef _OPENMP
1455 #pragma omp parallel
1456 #endif
1457 {
1458 vfloat exponentv = F2V(exponent);
1459 #ifdef _OPENMP
1460 #pragma omp for
1461 #endif
1462
1463 for(int i = 0; i < n - 3; i += 4) {
1464 STVFU(Source[i], xexpf(xlogf(LVFU(Source[i])) * exponentv));
1465 }
1466 }
1467
1468 for(int i = n - (n % 4); i < n; i++) {
1469 Source[i] = xexpf(xlogf(Source[i]) * exponent);
1470 }
1471
1472 #else
1473 #ifdef _OPENMP
1474 #pragma omp parallel for
1475 #endif
1476
1477 for(int i = 0; i < n; i++) {
1478 Source[i] = xexpf(xlogf(Source[i]) * exponent);
1479 }
1480
1481 #endif
1482
1483 }
1484
ContrastResid(float * WavCoeffs_L0,struct cont_params & cp,int W_L,int H_L,float max0,float min0)1485 void ImProcFunctions::ContrastResid(float * WavCoeffs_L0, struct cont_params &cp, int W_L, int H_L, float max0, float min0)
1486 {
1487 float stren = cp.tmstrength;
1488 float gamm = params->wavelet.gamma;
1489 cp.TMmeth = 2; //default after testing
1490
1491 if(cp.TMmeth == 1) {
1492 min0 = 0.0f;
1493 max0 = 32768.f;
1494 } else if (cp.TMmeth == 2) {
1495 min0 = 0.0f;
1496 }
1497
1498 #ifdef _OPENMP
1499 #pragma omp parallel for
1500 #endif
1501
1502 for(int i = 0; i < W_L * H_L; i++) {
1503 WavCoeffs_L0[i] = (WavCoeffs_L0[i] - min0) / max0;
1504 WavCoeffs_L0[i] *= gamm;
1505 }
1506
1507 float Compression = expf(-stren); //This modification turns numbers symmetric around 0 into exponents.
1508 float DetailBoost = stren;
1509
1510 if(stren < 0.0f) {
1511 DetailBoost = 0.0f; //Go with effect of exponent only if uncompressing.
1512 }
1513
1514
1515 CompressDR(WavCoeffs_L0, W_L, H_L, Compression, DetailBoost);
1516
1517
1518 #ifdef _OPENMP
1519 #pragma omp parallel for // removed schedule(dynamic,10)
1520 #endif
1521
1522 for(int ii = 0; ii < W_L * H_L; ii++) {
1523 WavCoeffs_L0[ii] = WavCoeffs_L0[ii] * max0 * (1.f / gamm) + min0;
1524 }
1525 }
1526
1527
1528
1529
EPDToneMapResid(float * WavCoeffs_L0,unsigned int Iterates,int skip,struct cont_params & cp,int W_L,int H_L,float max0,float min0)1530 void ImProcFunctions::EPDToneMapResid(float * WavCoeffs_L0, unsigned int Iterates, int skip, struct cont_params& cp, int W_L, int H_L, float max0, float min0)
1531 {
1532
1533
1534 float stren = cp.tmstrength;
1535 float edgest = params->epd.edgeStopping;
1536 float sca = params->epd.scale;
1537 float gamm = params->wavelet.gamma;
1538 float rew = params->epd.reweightingIterates;
1539 EdgePreservingDecomposition epd2(W_L, H_L);
1540 cp.TMmeth = 2; //default after testing
1541
1542 if(cp.TMmeth == 1) {
1543 min0 = 0.0f;
1544 max0 = 32768.f;
1545 } else if (cp.TMmeth == 2) {
1546 min0 = 0.0f;
1547 }
1548
1549 // max0=32768.f;
1550 #ifdef _OPENMP
1551 #pragma omp parallel for
1552 #endif
1553
1554 for(int i = 0; i < W_L * H_L; i++) {
1555 WavCoeffs_L0[i] = (WavCoeffs_L0[i] - min0) / max0;
1556 WavCoeffs_L0[i] *= gamm;
1557 }
1558
1559 float Compression = expf(-stren); //This modification turns numbers symmetric around 0 into exponents.
1560 float DetailBoost = stren;
1561
1562 if(stren < 0.0f) {
1563 DetailBoost = 0.0f; //Go with effect of exponent only if uncompressing.
1564 }
1565
1566 //Auto select number of iterates. Note that p->EdgeStopping = 0 makes a Gaussian blur.
1567 if(Iterates == 0) {
1568 Iterates = (unsigned int)(edgest * 15.0f);
1569 }
1570
1571
1572 epd2.CompressDynamicRange(WavCoeffs_L0, (float)sca / skip, edgest, Compression, DetailBoost, Iterates, rew);
1573
1574 //Restore past range, also desaturate a bit per Mantiuk's Color correction for tone mapping.
1575 #ifdef _OPENMP
1576 #pragma omp parallel for // removed schedule(dynamic,10)
1577 #endif
1578
1579 for(int ii = 0; ii < W_L * H_L; ii++) {
1580 WavCoeffs_L0[ii] = WavCoeffs_L0[ii] * max0 * (1.f / gamm) + min0;
1581 }
1582 }
1583
WaveletcontAllLfinal(wavelet_decomposition & WaveletCoeffs_L,struct cont_params & cp,float * mean,float * sigma,float * MaxP,const WavOpacityCurveWL & waOpacityCurveWL)1584 void ImProcFunctions::WaveletcontAllLfinal(wavelet_decomposition &WaveletCoeffs_L, struct cont_params &cp, float *mean, float *sigma, float *MaxP, const WavOpacityCurveWL & waOpacityCurveWL)
1585 {
1586 int maxlvl = WaveletCoeffs_L.maxlevel();
1587 float * WavCoeffs_L0 = WaveletCoeffs_L.coeff0;
1588
1589 for (int dir = 1; dir < 4; dir++) {
1590 for (int lvl = 0; lvl < maxlvl; lvl++) {
1591 int Wlvl_L = WaveletCoeffs_L.level_W(lvl);
1592 int Hlvl_L = WaveletCoeffs_L.level_H(lvl);
1593 float ** WavCoeffs_L = WaveletCoeffs_L.level_coeffs(lvl);
1594 finalContAllL (WavCoeffs_L, WavCoeffs_L0, lvl, dir, cp, Wlvl_L, Hlvl_L, mean, sigma, MaxP, waOpacityCurveWL);
1595 }
1596 }
1597 }
1598
1599
WaveletcontAllL(LabImage * labco,float ** varhue,float ** varchrom,wavelet_decomposition & WaveletCoeffs_L,struct cont_params & cp,int skip,float * mean,float * sigma,float * MaxP,float * MaxN,const WavCurve & wavCLVCcurve,const WavOpacityCurveW & waOpacityCurveW,FlatCurve * ChCurve,bool Chutili)1600 void ImProcFunctions::WaveletcontAllL(LabImage * labco, float ** varhue, float **varchrom, wavelet_decomposition &WaveletCoeffs_L,
1601 struct cont_params &cp, int skip, float *mean, float *sigma, float *MaxP, float *MaxN, const WavCurve & wavCLVCcurve, const WavOpacityCurveW & waOpacityCurveW, FlatCurve* ChCurve, bool Chutili)
1602 {
1603 int maxlvl = WaveletCoeffs_L.maxlevel();
1604 int W_L = WaveletCoeffs_L.level_W(0);
1605 int H_L = WaveletCoeffs_L.level_H(0);
1606 float * WavCoeffs_L0 = WaveletCoeffs_L.coeff0;
1607
1608 float maxh = 2.5f; //amplification contrast above mean
1609 float maxl = 2.5f; //reduction contrast under mean
1610 float contrast = cp.contrast;
1611 float multL = (float)contrast * (maxl - 1.f) / 100.f + 1.f;
1612 float multH = (float) contrast * (maxh - 1.f) / 100.f + 1.f;
1613 double avedbl = 0.0; // use double precision for large summations
1614 float max0 = 0.f;
1615 float min0 = FLT_MAX;
1616
1617 if(contrast != 0.f || (cp.tonemap && cp.resena)) { // contrast = 0.f means that all will be multiplied by 1.f, so we can skip this step
1618 #ifdef _OPENMP
1619 #pragma omp parallel for reduction(+:avedbl) num_threads(wavNestedLevels) if(wavNestedLevels>1)
1620 #endif
1621
1622 for (int i = 0; i < W_L * H_L; i++) {
1623 avedbl += WavCoeffs_L0[i];
1624 }
1625
1626 #ifdef _OPENMP
1627 #pragma omp parallel num_threads(wavNestedLevels) if(wavNestedLevels>1)
1628 #endif
1629 {
1630 float lminL = FLT_MAX;
1631 float lmaxL = 0.f;
1632
1633 #ifdef _OPENMP
1634 #pragma omp for
1635 #endif
1636
1637 for(int i = 0; i < W_L * H_L; i++) {
1638 if(WavCoeffs_L0[i] < lminL) {
1639 lminL = WavCoeffs_L0[i];
1640 }
1641
1642 if(WavCoeffs_L0[i] > lmaxL) {
1643 lmaxL = WavCoeffs_L0[i];
1644 }
1645
1646 }
1647
1648 #ifdef _OPENMP
1649 #pragma omp critical
1650 #endif
1651 {
1652 if(lminL < min0) {
1653 min0 = lminL;
1654 }
1655
1656 if(lmaxL > max0) {
1657 max0 = lmaxL;
1658 }
1659 }
1660
1661 }
1662
1663 }
1664
1665 // printf("MAXmax0=%f MINmin0=%f\n",max0,min0);
1666
1667 //tone mapping
1668 if(cp.tonemap && cp.contmet == 2 && cp.resena) {
1669 //iterate = 5
1670 EPDToneMapResid(WavCoeffs_L0, 0, skip, cp, W_L, H_L, max0, min0);
1671
1672 }
1673
1674 //end tonemapping
1675
1676
1677 max0 /= 327.68f;
1678 min0 /= 327.68f;
1679 float ave = avedbl / (double)(W_L * H_L);
1680 float av = ave / 327.68f;
1681 float ah = (multH - 1.f) / (av - max0); //
1682 float bh = 1.f - max0 * ah;
1683 float al = (multL - 1.f) / (av - min0);
1684 float bl = 1.f - min0 * al;
1685 float factorx = 1.f;
1686 // float *koeLi[9];
1687 // float maxkoeLi[9];
1688 float *koeLi[12];
1689 float maxkoeLi[12];
1690
1691 float *koeLibuffer = nullptr;
1692
1693 for(int y = 0; y < 12; y++) {
1694 maxkoeLi[y] = 0.f; //9
1695 }
1696
1697 koeLibuffer = new float[12 * H_L * W_L]; //12
1698
1699 for (int i = 0; i < 12; i++) { //9
1700 koeLi[i] = &koeLibuffer[i * W_L * H_L];
1701 }
1702
1703 for(int j = 0; j < 12; j++) //9
1704 for (int i = 0; i < W_L * H_L; i++) {
1705 koeLi[j][i] = 0.f;
1706 }
1707
1708 #ifdef _OPENMP
1709 #pragma omp parallel num_threads(wavNestedLevels) if(wavNestedLevels>1)
1710 #endif
1711 {
1712 if(contrast != 0.f && cp.resena && max0 > 0.0) { // contrast = 0.f means that all will be multiplied by 1.f, so we can skip this step
1713 {
1714 #ifdef _OPENMP
1715 #pragma omp for
1716 #endif
1717
1718 for (int i = 0; i < W_L * H_L; i++) { //contrast
1719 if(WavCoeffs_L0[i] < 32768.f) {
1720 float prov;
1721
1722 if( WavCoeffs_L0[i] > ave) {
1723 float kh = ah * (WavCoeffs_L0[i] / 327.68f) + bh;
1724 prov = WavCoeffs_L0[i];
1725 WavCoeffs_L0[i] = ave + kh * (WavCoeffs_L0[i] - ave);
1726 } else {
1727 float kl = al * (WavCoeffs_L0[i] / 327.68f) + bl;
1728 prov = WavCoeffs_L0[i];
1729 WavCoeffs_L0[i] = ave - kl * (ave - WavCoeffs_L0[i]);
1730 }
1731
1732 float diflc = WavCoeffs_L0[i] - prov;
1733 diflc *= factorx;
1734 WavCoeffs_L0[i] = prov + diflc;
1735 }
1736 }
1737 }
1738 }
1739
1740 if(cp.tonemap && cp.contmet == 1 && cp.resena) {
1741 float maxp = max0 * 256.f;
1742 float minp = min0 * 256.f;
1743 #ifdef _OPENMP
1744 #pragma omp single
1745 #endif
1746 ContrastResid(WavCoeffs_L0, cp, W_L, H_L, maxp, minp);
1747 }
1748
1749 #ifdef _OPENMP
1750 #pragma omp barrier
1751 #endif
1752
1753 if((cp.conres != 0.f || cp.conresH != 0.f) && cp.resena) { // cp.conres = 0.f and cp.comresH = 0.f means that all will be multiplied by 1.f, so we can skip this step
1754 #ifdef _OPENMP
1755 #pragma omp for nowait
1756 #endif
1757
1758 for (int i = 0; i < W_L * H_L; i++) {
1759 float LL = WavCoeffs_L0[i];
1760 float LL100 = LL / 327.68f;
1761 float tran = 5.f;//transition
1762 //shadow
1763 float alp = 3.f; //increase contrast sahdow in lowlights between 1 and ??
1764
1765 if(cp.th > (100.f - tran)) {
1766 tran = 100.f - cp.th;
1767 }
1768
1769 if(LL100 < cp.th) {
1770 float aalp = (1.f - alp) / cp.th; //no changes for LL100 = cp.th
1771 float kk = aalp * LL100 + alp;
1772 WavCoeffs_L0[i] *= (1.f + kk * cp.conres / 200.f);
1773 } else if(LL100 < cp.th + tran) {
1774 float ath = -cp.conres / tran;
1775 float bth = cp.conres - ath * cp.th;
1776 WavCoeffs_L0[i] *= (1.f + (LL100 * ath + bth) / 200.f);
1777 }
1778
1779 //highlight
1780 tran = 5.f;
1781
1782 if(cp.thH < (tran)) {
1783 tran = cp.thH;
1784 }
1785
1786 if(LL100 > cp.thH) {
1787 WavCoeffs_L0[i] *= (1.f + cp.conresH / 200.f);
1788 } else if(LL100 > (cp.thH - tran)) {
1789 float athH = cp.conresH / tran;
1790 float bthH = cp.conresH - athH * cp.thH;
1791 WavCoeffs_L0[i] *= (1.f + (LL100 * athH + bthH) / 200.f);
1792 }
1793 }
1794 }
1795
1796 //enabled Lipschitz..replace simple by complex edge detection
1797 // I found this concept on the web (doctoral thesis on medical Imaging)
1798 // I was inspired by the principle of Canny and Lipschitz (continuity and derivability)
1799 // I adapted the principle but have profoundly changed the algorithm
1800 // One can 1) change all parameters and found good parameters;
1801 //one can also change in calckoe
1802 float edd = 3.f;
1803 float eddlow = 15.f;
1804 float eddlipinfl = 0.005f * cp.edgsens + 0.4f;
1805 float eddlipampl = 1.f + cp.edgampl / 50.f;
1806
1807
1808 if(cp.detectedge) { //enabled Lipschitz control...more memory..more time...
1809 float *tmCBuffer = new float[H_L * W_L];
1810 float *tmC[H_L];
1811
1812 for (int i = 0; i < H_L; i++) {
1813 tmC[i] = &tmCBuffer[i * W_L];
1814 }
1815
1816 #ifdef _OPENMP
1817 #pragma omp for schedule(dynamic) collapse(2)
1818 #endif
1819
1820 for (int lvl = 0; lvl < 4; lvl++) {
1821 for (int dir = 1; dir < 4; dir++) {
1822 int W_L = WaveletCoeffs_L.level_W(lvl);
1823 int H_L = WaveletCoeffs_L.level_H(lvl);
1824
1825 float ** WavCoeffs_LL = WaveletCoeffs_L.level_coeffs(lvl);
1826 calckoe(WavCoeffs_LL, cp, koeLi, lvl , dir, W_L, H_L, edd, maxkoeLi, tmC);
1827 // return convolution KoeLi and maxkoeLi of level 0 1 2 3 and Dir Horiz, Vert, Diag
1828 }
1829 }
1830
1831 delete [] tmCBuffer;
1832
1833 float aamp = 1.f + cp.eddetthrHi / 100.f;
1834
1835 for (int lvl = 0; lvl < 4; lvl++) {
1836 #ifdef _OPENMP
1837 #pragma omp for schedule(dynamic,16)
1838 #endif
1839
1840 for (int i = 1; i < H_L - 1; i++) {
1841 for (int j = 1; j < W_L - 1; j++) {
1842 //treatment of koeLi and maxkoeLi
1843 float interm = 0.f;
1844
1845 if(cp.lip3 && cp.lipp) {
1846 // comparison between pixel and neighbours
1847 const auto neigh = cp.neigh == 1;
1848 const auto kneigh = neigh ? 28.f : 38.f;
1849 const auto somm = neigh ? 40.f : 50.f;
1850
1851 for (int dir = 1; dir < 4; dir++) { //neighbours proxi
1852 koeLi[lvl * 3 + dir - 1][i * W_L + j] = (kneigh * koeLi[lvl * 3 + dir - 1][i * W_L + j] + 2.f * koeLi[lvl * 3 + dir - 1][(i - 1) * W_L + j] + 2.f * koeLi[lvl * 3 + dir - 1][(i + 1) * W_L + j]
1853 + 2.f * koeLi[lvl * 3 + dir - 1][i * W_L + j + 1] + 2.f * koeLi[lvl * 3 + dir - 1][i * W_L + j - 1] + koeLi[lvl * 3 + dir - 1][(i - 1) * W_L + j - 1]
1854 + koeLi[lvl * 3 + dir - 1][(i - 1) * W_L + j + 1] + koeLi[lvl * 3 + dir - 1][(i + 1) * W_L + j - 1] + koeLi[lvl * 3 + dir - 1][(i + 1) * W_L + j + 1]) / somm;
1855 }
1856 }
1857
1858 for (int dir = 1; dir < 4; dir++) {
1859 //here I evaluate combinaison of vert / diag / horiz...we are with multiplicators of the signal
1860 interm += SQR(koeLi[lvl * 3 + dir - 1][i * W_L + j]);
1861 }
1862
1863 interm = sqrt(interm);
1864
1865 // interm /= 1.732f;//interm = pseudo variance koeLi
1866 interm *= 0.57736721f;
1867 float kampli = 1.f;
1868 float eps = 0.0001f;
1869 // I think this double ratio (alph, beta) is better than arctg
1870
1871 float alph = koeLi[lvl * 3][i * W_L + j] / (koeLi[lvl * 3 + 1][i * W_L + j] + eps); //ratio between horizontal and vertical
1872 float beta = koeLi[lvl * 3 + 2][i * W_L + j] / (koeLi[lvl * 3 + 1][i * W_L + j] + eps); //ratio between diagonal and horizontal
1873
1874 float alipinfl = (eddlipampl - 1.f) / (1.f - eddlipinfl);
1875 float blipinfl = eddlipampl - alipinfl;
1876
1877 //alph evaluate the direction of the gradient regularity Lipschitz
1878 // if = 1 we are on an edge
1879 // if 0 we are not
1880 // we can change and use log..or Arctg but why ?? we can change if need ...
1881 //Liamp=1 for eddlipinfl
1882 //liamp > 1 for alp >eddlipinfl and alph < 1
1883 //Liamp < 1 for alp < eddlipinfl and alph > 0
1884 if(alph > 1.f) {
1885 alph = 1.f / alph;
1886 }
1887
1888 if(beta > 1.f) {
1889 beta = 1.f / beta;
1890 }
1891
1892 //take into account diagonal
1893 //if in same value OK
1894 //if not no edge or reduction
1895 float bet = 1.f;
1896
1897 //if(cp.lip3) {//enhance algorithm
1898 if(alph > eddlipinfl && beta < 0.85f * eddlipinfl) { //0.85 arbitrary value ==> eliminate from edge if H V D too different
1899 bet = beta;
1900 }
1901
1902 //}
1903 float AmpLip = 1.f;
1904
1905 if(alph > eddlipinfl) {
1906 AmpLip = alipinfl * alph + blipinfl; //If beta low reduce kampli
1907 kampli = SQR(bet) * AmpLip * aamp;
1908 } else {
1909 AmpLip = (1.f / eddlipinfl) * SQR(SQR(alph * bet)); //Strong Reduce if beta low
1910 kampli = AmpLip / aamp;
1911 }
1912
1913 // comparison betwwen pixel and neighbours to do ==> I think 3 dir above is better
1914 /* if(cp.lip3){
1915 koeLi[lvl*3][i*W_L + j] = (koeLi[lvl*3][i*W_L + j] + koeLi[lvl*3][(i-1)*W_L + j] + koeLi[lvl*3][(i+1)*W_L + j]
1916 + koeLi[lvl*3][i*W_L + j+1] + koeLi[lvl*3][i*W_L + j-1] + koeLi[lvl*3][(i-1)*W_L + j-1]
1917 + koeLi[lvl*3][(i-1)*W_L + j+1] +koeLi[lvl*3][(i+1)*W_L + j-1] +koeLi[lvl*3][(i+1)*W_L + j+1])/9.f;
1918 }
1919 */
1920 // apply to each direction Wavelet level : horizontal / vertiacle / diagonal
1921 //interm += SQR(koeLi[lvl*3 + dir-1][i*W_L + j]);
1922
1923 interm *= kampli;
1924
1925 if(interm < cp.eddetthr / eddlow) {
1926 interm = 0.01f; //eliminate too low values
1927 }
1928
1929 //we can change this part of algo==> not equal but ponderate
1930 koeLi[lvl * 3][i * W_L + j] = koeLi[lvl * 3 + 1][i * W_L + j] = koeLi[lvl * 3 + 2][i * W_L + j] = interm; //new value
1931 //here KoeLi contains values where gradient is high and coef high, and eliminate low values...
1932 }
1933 }
1934 }
1935
1936 // end
1937 }
1938
1939 #ifdef _OPENMP
1940 #pragma omp for schedule(dynamic) collapse(2)
1941 #endif
1942
1943 for (int dir = 1; dir < 4; dir++) {
1944 for (int lvl = 0; lvl < maxlvl; lvl++) {
1945
1946 int Wlvl_L = WaveletCoeffs_L.level_W(lvl);
1947 int Hlvl_L = WaveletCoeffs_L.level_H(lvl);
1948
1949 float ** WavCoeffs_L = WaveletCoeffs_L.level_coeffs(lvl);
1950
1951 ContAllL (koeLi, maxkoeLi, true, maxlvl, labco, varhue, varchrom, WavCoeffs_L, WavCoeffs_L0, lvl, dir, cp, Wlvl_L, Hlvl_L, skip, mean, sigma, MaxP, MaxN, wavCLVCcurve, waOpacityCurveW, ChCurve, Chutili);
1952
1953
1954 }
1955 }
1956 }
1957
1958 //delete edge detection
1959 if(koeLibuffer) {
1960 delete [] koeLibuffer;
1961 }
1962 }
1963
WaveletAandBAllAB(wavelet_decomposition & WaveletCoeffs_a,wavelet_decomposition & WaveletCoeffs_b,struct cont_params & cp,FlatCurve * hhCurve,bool hhutili)1964 void ImProcFunctions::WaveletAandBAllAB(wavelet_decomposition &WaveletCoeffs_a, wavelet_decomposition &WaveletCoeffs_b,
1965 struct cont_params &cp, FlatCurve* hhCurve, bool hhutili)
1966 {
1967 // StopWatch Stop1("WaveletAandBAllAB");
1968 if (hhutili && cp.resena) { // H=f(H)
1969 int W_L = WaveletCoeffs_a.level_W(0);
1970 int H_L = WaveletCoeffs_a.level_H(0);
1971
1972 float * WavCoeffs_a0 = WaveletCoeffs_a.coeff0;
1973 float * WavCoeffs_b0 = WaveletCoeffs_b.coeff0;
1974 #ifdef _OPENMP
1975 #pragma omp parallel num_threads(wavNestedLevels) if(wavNestedLevels>1)
1976 #endif
1977 {
1978 #ifdef __SSE2__
1979 float huebuffer[W_L] ALIGNED64;
1980 float chrbuffer[W_L] ALIGNED64;
1981 #endif // __SSE2__
1982 #ifdef _OPENMP
1983 #pragma omp for schedule(dynamic,16)
1984 #endif
1985
1986 for (int i = 0; i < H_L; i++) {
1987 #ifdef __SSE2__
1988 // precalculate hue and chr
1989 int k;
1990
1991 for (k = 0; k < W_L - 3; k += 4) {
1992 __m128 av = LVFU(WavCoeffs_a0[i * W_L + k]);
1993 __m128 bv = LVFU(WavCoeffs_b0[i * W_L + k]);
1994 __m128 huev = xatan2f(bv, av);
1995 __m128 chrv = vsqrtf(SQRV(av) + SQRV(bv));
1996 STVF(huebuffer[k], huev);
1997 STVF(chrbuffer[k], chrv);
1998 }
1999
2000 for(; k < W_L; k++) {
2001 huebuffer[k] = xatan2f(WavCoeffs_b0[i * W_L + k], WavCoeffs_a0[i * W_L + k]);
2002 chrbuffer[k] = sqrtf(SQR(WavCoeffs_b0[i * W_L + k]) + SQR(WavCoeffs_a0[i * W_L + k])) / 327.68f;
2003 }
2004
2005 #endif // __SSE2__
2006
2007 for (int j = 0; j < W_L; j++) {
2008
2009 #ifdef __SSE2__
2010 float hueR = huebuffer[j];
2011 float chR = chrbuffer[j];
2012 #else
2013 float hueR = xatan2f(WavCoeffs_b0[i * W_L + j], WavCoeffs_a0[i * W_L + j]);
2014 float chR = sqrtf(SQR(WavCoeffs_b0[i * W_L + j]) + SQR(WavCoeffs_a0[i * W_L + j]));
2015 #endif
2016 /* if (editID == EUID_WW_HHCurve) {//H pipette
2017 float valpar =Color::huelab_to_huehsv2(hueR);
2018 editWhatever->v(i,j) = valpar;
2019 }
2020 */
2021 float valparam = float((hhCurve->getVal(Color::huelab_to_huehsv2(hueR)) - 0.5f) * 1.7f) + hueR; //get H=f(H) 1.7 optimisation !
2022 float2 sincosval = xsincosf(valparam);
2023 WavCoeffs_a0[i * W_L + j] = chR * sincosval.y;
2024 WavCoeffs_b0[i * W_L + j] = chR * sincosval.x;
2025 }
2026 }
2027 }
2028 }
2029
2030 }
2031
WaveletcontAllAB(LabImage * labco,float ** varhue,float ** varchrom,wavelet_decomposition & WaveletCoeffs_ab,const WavOpacityCurveW & waOpacityCurveW,struct cont_params & cp,const bool useChannelA)2032 void ImProcFunctions::WaveletcontAllAB(LabImage * labco, float ** varhue, float **varchrom, wavelet_decomposition &WaveletCoeffs_ab, const WavOpacityCurveW & waOpacityCurveW,
2033 struct cont_params &cp, const bool useChannelA)
2034 {
2035
2036 int maxlvl = WaveletCoeffs_ab.maxlevel();
2037 int W_L = WaveletCoeffs_ab.level_W(0);
2038 int H_L = WaveletCoeffs_ab.level_H(0);
2039
2040 float * WavCoeffs_ab0 = WaveletCoeffs_ab.coeff0;
2041
2042 #ifdef _OPENMP
2043 #pragma omp parallel num_threads(wavNestedLevels) if(wavNestedLevels>1)
2044 #endif
2045 {
2046 if(cp.chrores != 0.f && cp.resena) { // cp.chrores == 0.f means all will be multiplied by 1.f, so we can skip the processing of residual
2047
2048 #ifdef _OPENMP
2049 #pragma omp for nowait
2050 #endif
2051
2052 for (int i = 0; i < W_L * H_L; i++) {
2053 const float skyprot = cp.sky;
2054 //chroma
2055 int ii = i / W_L;
2056 int jj = i - ii * W_L;
2057 float modhue = varhue[ii][jj];
2058 float scale = 1.f;
2059
2060 if(skyprot > 0.f) {
2061 if((modhue < cp.t_ry && modhue > cp.t_ly)) {
2062 scale = (100.f - cp.sky) / 100.1f;
2063 } else if((modhue >= cp.t_ry && modhue < cp.b_ry)) {
2064 scale = (100.f - cp.sky) / 100.1f;
2065 float ar = (scale - 1.f) / (cp.t_ry - cp.b_ry);
2066 float br = scale - cp.t_ry * ar;
2067 scale = ar * modhue + br;
2068 } else if((modhue > cp.b_ly && modhue < cp.t_ly)) {
2069 scale = (100.f - cp.sky) / 100.1f;
2070 float al = (scale - 1.f) / (-cp.b_ly + cp.t_ly);
2071 float bl = scale - cp.t_ly * al;
2072 scale = al * modhue + bl;
2073 }
2074 } else if(skyprot < 0.f) {
2075 if((modhue > cp.t_ry || modhue < cp.t_ly)) {
2076 scale = (100.f + cp.sky) / 100.1f;
2077 }
2078
2079 /* else if((modhue >= cp.t_ry && modhue < cp.b_ry)) {
2080 scale=(100.f+cp.sky)/100.1f;
2081 float ar=(scale-1.f)/(cp.t_ry- cp.b_ry);
2082 float br=scale-cp.t_ry*ar;
2083 scale=ar*modhue+br;
2084 }
2085 else if((modhue > cp.b_ly && modhue < cp.t_ly)) {
2086 scale=(100.f+cp.sky)/100.1f;
2087 float al=(scale-1.f)/(-cp.b_ly + cp.t_ly);
2088 float bl=scale-cp.t_ly*al;
2089 scale=al*modhue+bl;
2090 }
2091 */
2092 }
2093
2094 WavCoeffs_ab0[i] *= (1.f + cp.chrores * (scale) / 100.f);
2095
2096 }
2097 }
2098
2099 if(cp.cbena && cp.resena) {//if user select Toning and color balance
2100
2101 #ifdef _OPENMP
2102 #pragma omp for nowait
2103 #endif
2104
2105 for (int i = 0; i < W_L * H_L; i++) {
2106 int ii = i / W_L;
2107 int jj = i - ii * W_L;
2108 float LL = (labco->L[ii * 2][jj * 2]) / 327.68f; //I use labco but I can use also WavCoeffs_L0 (more exact but more memory)
2109
2110 float sca = 1.f; //amplifer - reducter...about 1, but perhaps 0.6 or 1.3
2111
2112 if(useChannelA) {//green red (little magenta)
2113 //transition to avoid artifacts with 6 between 30 to 36 and 63 to 69
2114 float aa = (cp.grmed - cp.grlow) / 6.f;
2115 float bb = cp.grlow - 30.f * aa;
2116 float aaa = (cp.grhigh - cp.grmed) / 6.f;
2117 float bbb = cp.grmed - 63.f * aaa;
2118
2119 if(LL < 30.f) { //shadows
2120 WavCoeffs_ab0[i] += cp.grlow * (sca) * 300.f;
2121 } else if(LL >= 30.f && LL < 36.f) { //transition
2122 float tr = aa * LL + bb;
2123 WavCoeffs_ab0[i] += tr * (sca) * 300.f;
2124 } else if(LL >= 36.f && LL < 63.f) { //midtones
2125 WavCoeffs_ab0[i] += cp.grmed * (sca) * 300.f;
2126 } else if(LL >= 63.f && LL < 69.f) { //transition
2127 float trh = aaa * LL + bbb;
2128 WavCoeffs_ab0[i] += trh * (sca) * 300.f;
2129 } else if(LL >= 69.f) { //highlights
2130 WavCoeffs_ab0[i] += cp.grhigh * (sca) * 300.f;
2131 }
2132 } else { //blue yellow
2133 //transition with 6 between 30 to 36 and 63 to 69
2134 float aa1 = (cp.blmed - cp.bllow) / 6.f;
2135 float bb1 = cp.bllow - 30.f * aa1;
2136 float aaa1 = (cp.blhigh - cp.blmed) / 6.f;
2137 float bbb1 = cp.blmed - 63.f * aaa1;
2138
2139 if(LL < 30.f) {
2140 WavCoeffs_ab0[i] += cp.bllow * (sca) * 300.f;
2141 } else if(LL >= 30.f && LL < 36.f) {
2142 float tr1 = aa1 * LL + bb1;
2143 WavCoeffs_ab0[i] += tr1 * (sca) * 300.f;
2144 } else if(LL >= 36.f && LL < 63.f) {
2145 WavCoeffs_ab0[i] += cp.blmed * (sca) * 300.f;
2146 } else if(LL >= 63.f && LL < 69.f) {
2147 float trh1 = aaa1 * LL + bbb1;
2148 WavCoeffs_ab0[i] += trh1 * (sca) * 300.f;
2149 } else if(LL >= 69.f) {
2150 WavCoeffs_ab0[i] += cp.blhigh * (sca) * 300.f;
2151 }
2152 }
2153 }
2154 }
2155
2156 #ifdef _OPENMP
2157 #pragma omp for schedule(dynamic) collapse(2)
2158 #endif
2159
2160 for (int dir = 1; dir < 4; dir++) {
2161 for (int lvl = 0; lvl < maxlvl; lvl++) {
2162
2163 int Wlvl_ab = WaveletCoeffs_ab.level_W(lvl);
2164 int Hlvl_ab = WaveletCoeffs_ab.level_H(lvl);
2165
2166 float ** WavCoeffs_ab = WaveletCoeffs_ab.level_coeffs(lvl);
2167 ContAllAB (labco, maxlvl, varhue, varchrom, WavCoeffs_ab, WavCoeffs_ab0, lvl, dir, waOpacityCurveW, cp, Wlvl_ab, Hlvl_ab, useChannelA);
2168 }
2169 }
2170
2171
2172 }
2173 }
2174
2175 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2176 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2177
calckoe(float ** WavCoeffs_LL,const struct cont_params & cp,float * koeLi[12],int level,int dir,int W_L,int H_L,float edd,float * maxkoeLi,float ** tmC)2178 void ImProcFunctions::calckoe(float ** WavCoeffs_LL, const struct cont_params& cp, float *koeLi[12], int level, int dir, int W_L, int H_L, float edd, float *maxkoeLi, float **tmC)
2179 {
2180 int borderL = 2;
2181
2182 // printf("cpedth=%f\n",cp.eddetthr);
2183 if(cp.eddetthr < 30.f) {
2184 borderL = 1;
2185
2186 // I calculate coefficients with r size matrix 3x3 r=1 ; 5x5 r=2; 7x7 r=3
2187 /*
2188 float k[2*r][2*r];
2189 for(int i=1;i<=(2*r+1);i++) {
2190 for(int j=1;j<=(2*r+1);j++) {
2191 k[i][j]=(1.f/6.283*sigma*sigma)*exp(-SQR(i-r-1)+SQR(j-r-1)/2.f*SQR(sigma));
2192 }
2193 }
2194 //I could also use Gauss.h for 3x3
2195 // If necessary I can put a 7x7 matrix
2196 */
2197 for (int i = 1; i < H_L - 1; i++) { //sigma=0.55
2198 for (int j = 1; j < W_L - 1; j++) {
2199 tmC[i][j] = (8.94f * WavCoeffs_LL[dir][i * W_L + j] + 1.71f * (WavCoeffs_LL[dir][(i - 1) * W_L + j] + 1.71f * WavCoeffs_LL[dir][(i + 1) * W_L + j]
2200 + 1.71f * WavCoeffs_LL[dir][i * W_L + j + 1] + 1.71f * WavCoeffs_LL[dir][i * W_L + j - 1]) + 0.33f * WavCoeffs_LL[dir][(i - 1) * W_L + j - 1]
2201 + 0.33f * WavCoeffs_LL[dir][(i - 1) * W_L + j + 1] + 0.33f * WavCoeffs_LL[dir][(i + 1) * W_L + j - 1] + 0.33f * WavCoeffs_LL[dir][(i + 1) * W_L + j + 1]) * 0.0584795f;
2202 // apply to each direction Wavelet level : horizontal / vertiacle / diagonal
2203
2204
2205 }
2206 }
2207 } else if(cp.eddetthr >= 30.f && cp.eddetthr < 50.f) {
2208 borderL = 1;
2209
2210 for (int i = 1; i < H_L - 1; i++) { //sigma=0.85
2211 for (int j = 1; j < W_L - 1; j++) {
2212 tmC[i][j] = (4.0091f * WavCoeffs_LL[dir][i * W_L + j] + 2.0068f * (WavCoeffs_LL[dir][(i - 1) * W_L + j] + 2.0068f * WavCoeffs_LL[dir][(i + 1) * W_L + j]
2213 + 2.0068f * WavCoeffs_LL[dir][i * W_L + j + 1] + 2.0068f * WavCoeffs_LL[dir][i * W_L + j - 1]) + 1.0045f * WavCoeffs_LL[dir][(i - 1) * W_L + j - 1]
2214 + 1.0045f * WavCoeffs_LL[dir][(i - 1) * W_L + j + 1] + 1.0045f * WavCoeffs_LL[dir][(i + 1) * W_L + j - 1] + 1.0045f * WavCoeffs_LL[dir][(i + 1) * W_L + j + 1]) * 0.062288f;
2215 // apply to each direction Wavelet level : horizontal / vertiacle / diagonal
2216
2217
2218 }
2219 }
2220 }
2221
2222
2223 else if(cp.eddetthr >= 50.f && cp.eddetthr < 75.f) {
2224 borderL = 1;
2225
2226 for (int i = 1; i < H_L - 1; i++) {
2227 for (int j = 1; j < W_L - 1; j++) { //sigma=1.1
2228 tmC[i][j] = (3.025f * WavCoeffs_LL[dir][i * W_L + j] + 2.001f * (WavCoeffs_LL[dir][(i - 1) * W_L + j] + 2.001f * WavCoeffs_LL[dir][(i + 1) * W_L + j]
2229 + 2.001f * WavCoeffs_LL[dir][i * W_L + j + 1] + 2.001f * WavCoeffs_LL[dir][i * W_L + j - 1]) + 1.323f * WavCoeffs_LL[dir][(i - 1) * W_L + j - 1]
2230 + 1.323f * WavCoeffs_LL[dir][(i - 1) * W_L + j + 1] + 1.323f * WavCoeffs_LL[dir][(i + 1) * W_L + j - 1] + 1.323f * WavCoeffs_LL[dir][(i + 1) * W_L + j + 1]) * 0.06127f;
2231 }
2232 }
2233 }
2234
2235 else if(cp.eddetthr >= 75.f) {
2236 borderL = 2;
2237
2238 //if(cp.lip3 && level > 1) {
2239 if(level > 1) {// do not activate 5x5 if level 0 or 1
2240
2241 for (int i = 2; i < H_L - 2; i++) {
2242 for (int j = 2; j < W_L - 2; j++) {
2243 // Gaussian 1.1
2244 // 0.5 2 3 2 0.5
2245 // 2 7 10 7 2
2246 // 3 10 15 10 3
2247 // 2 7 10 7 2
2248 // 0.5 2 3 2 0.5
2249 // divi 113
2250 //Gaussian 1.4
2251 // 2 4 5 4 2
2252 // 4 9 12 9 4
2253 // 5 12 15 12 5
2254 // 4 9 12 9 4
2255 // 2 4 5 4 2
2256 // divi 159
2257 if(cp.eddetthr < 85.f) { //sigma=1.1
2258 tmC[i][j] = (15.f * WavCoeffs_LL[dir][i * W_L + j] + 10.f * WavCoeffs_LL[dir][(i - 1) * W_L + j] + 10.f * WavCoeffs_LL[dir][(i + 1) * W_L + j]
2259 + 10.f * WavCoeffs_LL[dir][i * W_L + j + 1] + 10.f * WavCoeffs_LL[dir][i * W_L + j - 1] + 7.f * WavCoeffs_LL[dir][(i - 1) * W_L + j - 1]
2260 + 7.f * WavCoeffs_LL[dir][(i - 1) * W_L + j + 1] + 7.f * WavCoeffs_LL[dir][(i + 1) * W_L + j - 1] + 7.f * WavCoeffs_LL[dir][(i + 1) * W_L + j + 1]
2261 + 3.f * WavCoeffs_LL[dir][(i - 2) * W_L + j] + 3.f * WavCoeffs_LL[dir][(i + 2) * W_L + j] + 3.f * WavCoeffs_LL[dir][i * W_L + j - 2] + 3.f * WavCoeffs_LL[dir][i * W_L + j + 2]
2262 + 2.f * WavCoeffs_LL[dir][(i - 2) * W_L + j - 1] + 2.f * WavCoeffs_LL[dir][(i - 2) * W_L + j + 1] + 2.f * WavCoeffs_LL[dir][(i + 2) * W_L + j + 1] + 2.f * WavCoeffs_LL[dir][(i + 2) * W_L + j - 1]
2263 + 2.f * WavCoeffs_LL[dir][(i - 1) * W_L + j - 2] + 2.f * WavCoeffs_LL[dir][(i - 1) * W_L + j + 2] + 2.f * WavCoeffs_LL[dir][(i + 1) * W_L + j + 2] + 2.f * WavCoeffs_LL[dir][(i + 1) * W_L + j - 2]
2264 + 0.5f * WavCoeffs_LL[dir][(i - 2) * W_L + j - 2] + 0.5f * WavCoeffs_LL[dir][(i - 2) * W_L + j + 2] + 0.5f * WavCoeffs_LL[dir][(i + 2) * W_L + j - 2] + 0.5f * WavCoeffs_LL[dir][(i + 2) * W_L + j + 2]
2265 ) * 0.0088495f;
2266
2267 }
2268
2269 else {//sigma=1.4
2270 tmC[i][j] = (15.f * WavCoeffs_LL[dir][i * W_L + j] + 12.f * WavCoeffs_LL[dir][(i - 1) * W_L + j] + 12.f * WavCoeffs_LL[dir][(i + 1) * W_L + j]
2271 + 12.f * WavCoeffs_LL[dir][i * W_L + j + 1] + 12.f * WavCoeffs_LL[dir][i * W_L + j - 1] + 9.f * WavCoeffs_LL[dir][(i - 1) * W_L + j - 1]
2272 + 9.f * WavCoeffs_LL[dir][(i - 1) * W_L + j + 1] + 9.f * WavCoeffs_LL[dir][(i + 1) * W_L + j - 1] + 9.f * WavCoeffs_LL[dir][(i + 1) * W_L + j + 1]
2273 + 5.f * WavCoeffs_LL[dir][(i - 2) * W_L + j] + 5.f * WavCoeffs_LL[dir][(i + 2) * W_L + j] + 5.f * WavCoeffs_LL[dir][i * W_L + j - 2] + 5.f * WavCoeffs_LL[dir][i * W_L + j + 2]
2274 + 4.f * WavCoeffs_LL[dir][(i - 2) * W_L + j - 1] + 4.f * WavCoeffs_LL[dir][(i - 2) * W_L + j + 1] + 4.f * WavCoeffs_LL[dir][(i + 2) * W_L + j + 1] + 4.f * WavCoeffs_LL[dir][(i + 2) * W_L + j - 1]
2275 + 4.f * WavCoeffs_LL[dir][(i - 1) * W_L + j - 2] + 4.f * WavCoeffs_LL[dir][(i - 1) * W_L + j + 2] + 4.f * WavCoeffs_LL[dir][(i + 1) * W_L + j + 2] + 4.f * WavCoeffs_LL[dir][(i + 1) * W_L + j - 2]
2276 + 2.f * WavCoeffs_LL[dir][(i - 2) * W_L + j - 2] + 2.f * WavCoeffs_LL[dir][(i - 2) * W_L + j + 2] + 2.f * WavCoeffs_LL[dir][(i + 2) * W_L + j - 2] + 2.f * WavCoeffs_LL[dir][(i + 2) * W_L + j + 2]
2277 ) * 0.0062893f;
2278 }
2279
2280
2281 // apply to each direction Wavelet level : horizontal / vertiacle / diagonal
2282 }
2283 }
2284 }
2285
2286 }
2287
2288
2289 /*
2290 // I suppress these 2 convolutions ==> lees good results==> probably because structure data different and also I compare to original value which have + and -
2291 for(int i = borderL; i < H_L-borderL; i++ ) {//[-1 0 1] x==>j
2292 for(int j = borderL; j < W_L-borderL; j++) {
2293 tmC[i][j]=- WavCoeffs_LL[dir][(i)*W_L + j-1] + WavCoeffs_LL[dir][(i)*W_L + j+1];
2294 }
2295 }
2296 for(int i = borderL; i < H_L-borderL; i++ ) {//[1 0 -1] y==>i
2297 for(int j = borderL; j < W_L-borderL; j++) {
2298 tmC[i][j]= - WavCoeffs_LL[dir][(i-1)*W_L + j] + WavCoeffs_LL[dir][(i+1)*W_L + j];
2299 }
2300 }
2301 */
2302
2303 float thr = 40.f; //avoid artifact eg. noise...to test
2304 float thr2 = 1.5f * edd; //edd can be modified in option ed_detect
2305 thr2 += cp.eddet / 30.f; //to test
2306 float diffFactor = (cp.eddet / 100.f);
2307
2308 for(int i = 0; i < H_L; i++ ) {
2309 for(int j = 0; j < W_L; j++) {
2310 koeLi[level * 3 + dir - 1][i * W_L + j] = 1.f;
2311 }
2312 }
2313
2314 for(int i = borderL; i < H_L - borderL; i++ ) {
2315 for(int j = borderL; j < W_L - borderL; j++) {
2316 // my own algo : probably a little false, but simpler as Lipschitz !
2317 // Thr2 = maximum of the function ==> Lipsitch says = probably edge
2318 // float temp = WavCoeffs_LL[dir][i*W_L + j];
2319 // if(temp>=0.f && temp < thr) temp = thr;
2320 // if(temp < 0.f && temp > -thr) temp = -thr;
2321 float temp = max(fabsf(WavCoeffs_LL[dir][i * W_L + j]), thr );
2322 koeLi[level * 3 + dir - 1][i * W_L + j] = min(thr2, fabs(tmC[i][j] / temp)); // limit maxi
2323
2324 //it will be more complicated to calculate both Wh and Wv, but we have also Wd==> pseudo Lipschitz
2325 if(koeLi[level * 3 + dir - 1][i * W_L + j] > maxkoeLi[level * 3 + dir - 1]) {
2326 maxkoeLi[level * 3 + dir - 1] = koeLi[level * 3 + dir - 1][i * W_L + j];
2327 }
2328
2329 float diff = maxkoeLi[level * 3 + dir - 1] - koeLi[level * 3 + dir - 1][i * W_L + j];
2330 diff *= diffFactor;
2331 koeLi[level * 3 + dir - 1][i * W_L + j] = maxkoeLi[level * 3 + dir - 1] - diff;
2332 }
2333 }
2334
2335 }
2336
finalContAllL(float ** WavCoeffs_L,float * WavCoeffs_L0,int level,int dir,struct cont_params & cp,int W_L,int H_L,float * mean,float * sigma,float * MaxP,const WavOpacityCurveWL & waOpacityCurveWL)2337 void ImProcFunctions::finalContAllL (float ** WavCoeffs_L, float * WavCoeffs_L0, int level, int dir, struct cont_params &cp,
2338 int W_L, int H_L, float *mean, float *sigma, float *MaxP, const WavOpacityCurveWL & waOpacityCurveWL)
2339 {
2340 if(cp.diagcurv && cp.finena && MaxP[level] > 0.f && mean[level] != 0.f && sigma[level] != 0.f ) { //curve
2341 float insigma = 0.666f; //SD
2342 float logmax = log(MaxP[level]); //log Max
2343 float rapX = (mean[level] + sigma[level]) / MaxP[level]; //rapport between sD / max
2344 float inx = log(insigma);
2345 float iny = log(rapX);
2346 float rap = inx / iny; //koef
2347 float asig = 0.166f / sigma[level];
2348 float bsig = 0.5f - asig * mean[level];
2349 float amean = 0.5f / mean[level];
2350
2351 #ifdef _OPENMP
2352 #pragma omp parallel for schedule(dynamic, W_L * 16) num_threads(wavNestedLevels) if(wavNestedLevels>1)
2353 #endif
2354
2355 for (int i = 0; i < W_L * H_L; i++) {
2356 float absciss;
2357
2358 if(fabsf(WavCoeffs_L[dir][i]) >= (mean[level] + sigma[level])) { //for max
2359 float valcour = xlogf(fabsf(WavCoeffs_L[dir][i]));
2360 float valc = valcour - logmax;
2361 float vald = valc * rap;
2362 absciss = xexpf(vald);
2363 } else if(fabsf(WavCoeffs_L[dir][i]) >= mean[level]) {
2364 absciss = asig * fabsf(WavCoeffs_L[dir][i]) + bsig;
2365 } else {
2366 absciss = amean * fabsf(WavCoeffs_L[dir][i]);
2367 }
2368
2369 float kc = waOpacityCurveWL[absciss * 500.f] - 0.5f;
2370 float reduceeffect = kc <= 0.f ? 1.f : 1.5f;
2371
2372 float kinterm = 1.f + reduceeffect * kc;
2373 kinterm = kinterm <= 0.f ? 0.01f : kinterm;
2374
2375 WavCoeffs_L[dir][i] *= kinterm;
2376 }
2377 }
2378
2379 int choicelevel = params->wavelet.Lmethod - 1;
2380 choicelevel = choicelevel == -1 ? 4 : choicelevel;
2381
2382 int choiceClevel = 0;
2383
2384 if(params->wavelet.CLmethod == "one") {
2385 choiceClevel = 0;
2386 } else if(params->wavelet.CLmethod == "inf") {
2387 choiceClevel = 1;
2388 } else if(params->wavelet.CLmethod == "sup") {
2389 choiceClevel = 2;
2390 } else if(params->wavelet.CLmethod == "all") {
2391 choiceClevel = 3;
2392 }
2393
2394 int choiceDir = 0;
2395
2396 if(params->wavelet.Dirmethod == "one") {
2397 choiceDir = 1;
2398 } else if(params->wavelet.Dirmethod == "two") {
2399 choiceDir = 2;
2400 } else if(params->wavelet.Dirmethod == "thr") {
2401 choiceDir = 3;
2402 } else if(params->wavelet.Dirmethod == "all") {
2403 choiceDir = 0;
2404 }
2405
2406 int dir1 = (choiceDir == 2) ? 1 : 2;
2407 int dir2 = (choiceDir == 3) ? 1 : 3;
2408
2409 if(choiceClevel < 3) { // not all levels visible, paint residual
2410 if(level == 0) {
2411 if(cp.backm != 2) { // nothing to change when residual is used as background
2412 float backGroundColor = (cp.backm == 1) ? 12000.f : 0.f;
2413
2414 for (int i = 0; i < W_L * H_L; i++) {
2415 WavCoeffs_L0[i] = backGroundColor;
2416 }
2417 }
2418 }
2419 }
2420
2421 if(choiceClevel == 0) { // Only one level
2422
2423 if(choiceDir == 0) { // All directions
2424 if(level != choicelevel) { // zero all for the levels != choicelevel
2425 for (int dir = 1; dir < 4; dir++) {
2426 for (int i = 0; i < W_L * H_L; i++) {
2427 WavCoeffs_L[dir][i] = 0.f;
2428 }
2429 }
2430 }
2431 } else { // zero the unwanted directions for level == choicelevel
2432
2433 if(choicelevel >= cp.maxilev) {
2434 for (int dir = 1; dir < 4; dir++) {
2435 for (int i = 0; i < W_L * H_L; i++) {
2436 WavCoeffs_L[dir][i] = 0.f;
2437 }
2438 }
2439 } else if(level != choicelevel) { // zero all for the levels != choicelevel
2440 for (int i = 0; i < W_L * H_L; i++) {
2441 WavCoeffs_L[dir1][i] = WavCoeffs_L[dir2][i] = 0.f;
2442 }
2443 }
2444 }
2445 } else if(choiceClevel == 1) { // Only below level
2446 if(choiceDir == 0) { // All directions
2447 if(level > choicelevel) {
2448 for (int dir = 1; dir < 4; dir++) {
2449 for (int i = 0; i < W_L * H_L; i++) {
2450 WavCoeffs_L[dir][i] = 0.f;
2451 }
2452 }
2453 }
2454 } else { // zero the unwanted directions for level >= choicelevel
2455 if(level > choicelevel) {
2456 for (int i = 0; i < W_L * H_L; i++) {
2457 WavCoeffs_L[dir1][i] = WavCoeffs_L[dir2][i] = 0.f;
2458 }
2459 }
2460 }
2461 } else if(choiceClevel == 2) { // Only above level
2462 if(choiceDir == 0) { // All directions
2463 if(level <= choicelevel) {
2464 for (int dir = 1; dir < 4; dir++) {
2465 for (int i = 0; i < W_L * H_L; i++) {
2466 WavCoeffs_L[dir][i] = 0.f;
2467 }
2468 }
2469 }
2470 } else { // zero the unwanted directions for level >= choicelevel
2471 if(choicelevel >= cp.maxilev) {
2472 for (int dir = 1; dir < 4; dir++) {
2473 for (int i = 0; i < W_L * H_L; i++) {
2474 WavCoeffs_L[dir][i] = 0.f;
2475 }
2476 }
2477 }
2478
2479
2480 else if(level <= choicelevel) {
2481 for (int i = 0; i < W_L * H_L; i++) {
2482 WavCoeffs_L[dir1][i] = WavCoeffs_L[dir2][i] = 0.f;
2483 }
2484 }
2485 }
2486 }
2487
2488
2489 }
2490
ContAllL(float * koeLi[12],float * maxkoeLi,bool lipschitz,int maxlvl,LabImage * labco,float ** varhue,float ** varchrom,float ** WavCoeffs_L,float * WavCoeffs_L0,int level,int dir,struct cont_params & cp,int W_L,int H_L,int skip,float * mean,float * sigma,float * MaxP,float * MaxN,const WavCurve & wavCLVCcurve,const WavOpacityCurveW & waOpacityCurveW,FlatCurve * ChCurve,bool Chutili)2491 void ImProcFunctions::ContAllL (float *koeLi[12], float *maxkoeLi, bool lipschitz, int maxlvl, LabImage * labco, float ** varhue, float **varchrom, float ** WavCoeffs_L, float * WavCoeffs_L0, int level, int dir, struct cont_params &cp,
2492 int W_L, int H_L, int skip, float *mean, float *sigma, float *MaxP, float *MaxN, const WavCurve & wavCLVCcurve, const WavOpacityCurveW & waOpacityCurveW, FlatCurve* ChCurve, bool Chutili)
2493 {
2494 assert (level >= 0);
2495 assert (maxlvl > level);
2496
2497 static const float scales[10] = {1.f, 2.f, 4.f, 8.f, 16.f, 32.f, 64.f, 128.f, 256.f, 512.f};
2498 float scaleskip[10];
2499
2500 for(int sc = 0; sc < 10; sc++) {
2501 scaleskip[sc] = scales[sc] / skip;
2502 }
2503
2504 float t_r = 40.f;
2505 float t_l = 10.f;
2506 float b_r = 75.f;
2507 float edd = 3.f;
2508 float eddstrength = 1.3f;
2509 float aedstr = (eddstrength - 1.f) / 90.f;
2510 float bedstr = 1.f - 10.f * aedstr;
2511
2512 if(cp.val > 0 && cp.edgeena) {
2513 float * koe = nullptr;
2514 float maxkoe = 0.f;
2515
2516 if(!lipschitz) {
2517 koe = new float [H_L * W_L];
2518
2519 for (int i = 0; i < W_L * H_L; i++) {
2520 koe[i] = 0.f;
2521 }
2522
2523 maxkoe = 0.f;
2524
2525 if(cp.detectedge) {
2526 float** tmC;
2527 int borderL = 1;
2528 tmC = new float*[H_L];
2529
2530 for (int i = 0; i < H_L; i++) {
2531 tmC[i] = new float[W_L];
2532 }
2533
2534 {
2535 for (int i = 1; i < H_L - 1; i++) {
2536 for (int j = 1; j < W_L - 1; j++) {
2537 //edge detection wavelet TMC Canny
2538 // also possible to detect noise with 5x5 instead of 3x3
2539 tmC[i][j] = (4.f * WavCoeffs_L[dir][i * W_L + j] + 2.f * WavCoeffs_L[dir][(i - 1) * W_L + j] + 2.f * WavCoeffs_L[dir][(i + 1) * W_L + j]
2540 + 2.f * WavCoeffs_L[dir][i * W_L + j + 1] + 2.f * WavCoeffs_L[dir][i * W_L + j - 1] + WavCoeffs_L[dir][(i - 1) * W_L + j - 1]
2541 + WavCoeffs_L[dir][(i - 1) * W_L + j + 1] + WavCoeffs_L[dir][(i + 1) * W_L + j - 1] + WavCoeffs_L[dir][(i + 1) * W_L + j + 1]) / 16.f;
2542
2543 // apply to each direction Wavelet level : horizontal / vertiacle / diagonal
2544 }
2545 }
2546 }
2547
2548
2549
2550
2551 for(int i = borderL; i < H_L - borderL; i++ ) {
2552 for(int j = borderL; j < W_L - borderL; j++) {
2553 // my own algo : probably a little false, but simpler as Lipschitz !
2554 float thr = 40.f; //avoid artifact eg. noise...to test
2555 float thr2 = edd; //edd can be modified in option ed_detect
2556 thr2 += cp.eddet / 30.f; //to test
2557 float temp = WavCoeffs_L[dir][i * W_L + j];
2558
2559 if(temp >= 0.f && temp < thr) {
2560 temp = thr;
2561 }
2562
2563 if(temp < 0.f && temp > -thr) {
2564 temp = -thr;
2565 }
2566
2567 koe[i * W_L + j] = min(thr2, fabs(tmC[i][j] / temp));
2568
2569 if(koe[i * W_L + j] > maxkoe) {
2570 maxkoe = koe[i * W_L + j];
2571 }
2572
2573 float diff = maxkoe - koe[i * W_L + j];
2574 diff *= (cp.eddet / 100.f);
2575 float interm = maxkoe - diff;
2576
2577 if(interm < cp.eddetthr / 30.f) {
2578 interm = 0.01f;
2579 }
2580
2581 koe[i * W_L + j] = interm;
2582
2583 }
2584 }
2585
2586 for (int i = 0; i < H_L; i++) {
2587 delete [] tmC[i];
2588 }
2589
2590 delete [] tmC;
2591
2592 }
2593 }
2594
2595 //end detect edge
2596 float rad = ((float)cp.rad) / 60.f; //radius ==> not too high value to avoid artifacts
2597 float value = ((float)cp.val) / 8.f; //strength
2598
2599 if (scaleskip[1] < 1.f) {
2600 float atten01234 = 0.80f;
2601 value *= (atten01234 * scaleskip[1]); //for zoom < 100% reduce strength...I choose level 1...but!!
2602 }
2603
2604 float edge = 1.f;
2605 float lim0 = 20.f; //arbitrary limit for low radius and level between 2 or 3 to 30 maxi
2606 float lev = float (level);
2607 float repart = (float)cp.til;
2608 float brepart;
2609
2610 if(cp.reinforce == 1) {
2611 brepart = 3.f;
2612 }
2613
2614 if(cp.reinforce == 3) {
2615 brepart = 0.5f; //arbitrary value to increase / decrease repart, between 1 and 0
2616 }
2617
2618 float arepart = -(brepart - 1.f) / (lim0 / 60.f);
2619
2620 if (cp.reinforce != 2) {
2621 if(rad < lim0 / 60.f) {
2622 repart *= (arepart * rad + brepart); //linear repartition of repart
2623 }
2624 }
2625
2626 float al0 = 1.f + (repart) / 50.f;
2627 float al10 = 1.0f; //arbitrary value ==> less = take into account high levels
2628 // float ak =-(al0-al10)/10.f;//10 = maximum levels
2629 float ak = -(al0 - al10) / 10.f; //10 = maximum levels
2630 float bk = al0;
2631 float koef = ak * level + bk; //modulate for levels : more levels high, more koef low ==> concentrated action on low levels, without or near for high levels
2632 float expkoef = -pow(fabs(rad - lev), koef); //reduce effect for high levels
2633
2634 if (cp.reinforce == 3) {
2635 if(rad < lim0 / 60.f && level == 0) {
2636 expkoef *= abs(repart); //reduce effect for low values of rad and level=0==> quasi only level 1 is effective
2637 }
2638 }
2639
2640 if (cp.reinforce == 1) {
2641 if(rad < lim0 / 60.f && level == 1) {
2642 expkoef /= repart; //increase effect for low values of rad and level=1==> quasi only level 0 is effective
2643 }
2644 }
2645
2646 //take into account local contrast
2647 float refin = value * exp (expkoef);
2648
2649 if(cp.link && cp.noiseena) { //combi
2650 {
2651 if(level == 0) {
2652 refin *= (1.f + cp.lev0s / 50.f); // we can change this sensibility!
2653 }
2654
2655 if(level == 1) {
2656 refin *= (1.f + cp.lev1s / 50.f);
2657 }
2658
2659 if(level == 2) {
2660 refin *= (1.f + cp.lev2s / 50.f);
2661 }
2662
2663 if(level == 3) {
2664 refin *= (1.f + cp.lev3s / 50.f);
2665 }
2666 }
2667 }
2668
2669 float edgePrecalc = 1.f + refin; //estimate edge "pseudo variance"
2670
2671 if(cp.EDmet == 2 && MaxP[level] > 0.f) { //curve
2672 // if(exa) {//curve
2673 float insigma = 0.666f; //SD
2674 float logmax = log(MaxP[level]); //log Max
2675 float rapX = (mean[level] + sigma[level]) / MaxP[level]; //rapport between sD / max
2676 float inx = log(insigma);
2677 float iny = log(rapX);
2678 float rap = inx / iny; //koef
2679 float asig = 0.166f / sigma[level];
2680 float bsig = 0.5f - asig * mean[level];
2681 float amean = 0.5f / mean[level];
2682 float absciss = 0.f;
2683 float kinterm;
2684 float kmul;
2685 int borderL = 1;
2686
2687 for(int i = borderL; i < H_L - borderL; i++ ) {
2688 for(int j = borderL; j < W_L - borderL; j++) {
2689 int k = i * W_L + j;
2690
2691 if(cp.detectedge) {
2692 if(!lipschitz) {
2693 if(cp.eddet > 10.f) {
2694 edge = (aedstr * cp.eddet + bedstr) * (edgePrecalc * (1.f + koe[k])) / (1.f + 0.9f * maxkoe);
2695 } else {
2696 edge = (edgePrecalc * (1.f + koe[k])) / (1.f + 0.9f * maxkoe);
2697 }
2698 }
2699
2700 if(lipschitz) {
2701 if(level < 4) {
2702 edge = 1.f + (edgePrecalc - 1.f) * (koeLi[level * 3][k]) / (1.f + 0.9f * maxkoeLi[level * 3 + dir - 1]);
2703 } else {
2704 edge = edgePrecalc;
2705 }
2706 }
2707 } else {
2708 edge = edgePrecalc;
2709 }
2710
2711 if(cp.edgcurv) {
2712 if(fabs(WavCoeffs_L[dir][k]) >= (mean[level] + sigma[level])) { //for max
2713 float valcour = log(fabs(WavCoeffs_L[dir][k]));
2714 float valc = valcour - logmax;
2715 float vald = valc * rap;
2716 absciss = exp(vald);
2717
2718 } else if(fabs(WavCoeffs_L[dir][k]) >= mean[level] && fabs(WavCoeffs_L[dir][k]) < (mean[level] + sigma[level])) {
2719 absciss = asig * fabs(WavCoeffs_L[dir][k]) + bsig;
2720 } else if(fabs(WavCoeffs_L[dir][k]) < mean[level]) {
2721 absciss = amean * fabs(WavCoeffs_L[dir][k]);
2722 }
2723
2724 // Threshold adjuster settings==> approximative for curve
2725 //kmul about average cbrt(3--40 / 10)==>1.5 to 2.5
2726 //kmul about SD 10--60 / 35 ==> 2
2727 // kmul about low cbrt((5.f+cp.edg_low)/5.f);==> 1.5
2728 // kmul about max ==> 9
2729 // we can change these values
2730 // result is different not best or bad than threshold slider...but similar
2731 float abssd = 4.f; //amplification reference
2732 float bbssd = 2.f; //mini ampli
2733 float maxamp = 2.5f; //maxi ampli at end
2734 float maxampd = 10.f; //maxi ampli at end
2735 float a_abssd = (maxamp - abssd) / 0.333f;
2736 float b_abssd = maxamp - a_abssd;
2737 float da_abssd = (maxampd - abssd) / 0.333f;
2738 float db_abssd = maxampd - da_abssd;
2739 float am = (abssd - bbssd) / 0.666f;
2740 float kmuld = 0.f;
2741
2742 if(absciss > 0.666f && absciss < 1.f) {
2743 kmul = a_abssd * absciss + b_abssd; //about max ==> kinterm
2744 kmuld = da_abssd * absciss + db_abssd;
2745 } else {
2746 kmul = kmuld = absciss * am + bbssd;
2747 }
2748
2749 kinterm = 1.f;
2750 float kc = kmul * (wavCLVCcurve[absciss * 500.f] - 0.5f);
2751 float kcd = kmuld * (wavCLVCcurve[absciss * 500.f] - 0.5f);
2752
2753 if(kc >= 0.f) {
2754 float reduceeffect = 0.6f;
2755 kinterm = 1.f + reduceeffect * kmul * (wavCLVCcurve[absciss * 500.f] - 0.5f); //about 1 to 3 general and big amplification for max (under 0)
2756 } else {
2757 kinterm = 1.f - (SQR(kcd)) / 10.f;
2758 }
2759
2760 if(kinterm < 0.f) {
2761 kinterm = 0.01f;
2762 }
2763
2764 edge *= kinterm;
2765
2766 if(edge < 1.f) {
2767 edge = 1.f;
2768 }
2769 }
2770
2771 WavCoeffs_L[dir][k] *= edge;
2772 }
2773 }
2774 } else if(cp.EDmet == 1) { //threshold adjuster
2775 float MaxPCompare = MaxP[level] * SQR(cp.edg_max / 100.f); //100 instead of b_r...case if b_r < 100
2776 float MaxNCompare = MaxN[level] * SQR(cp.edg_max / 100.f); //always reduce a little edge for near max values
2777 float edgeSdCompare = (mean[level] + 1.5f * sigma[level]) * SQR(cp.edg_sd / t_r); // 1.5 standard deviation #80% range between mean 50% and 80%
2778 float edgeMeanCompare = mean[level] * SQR(cp.edg_mean / t_l);
2779 float edgeLowCompare = (5.f + SQR(cp.edg_low));
2780 float edgeMeanFactor = cbrt(cp.edg_mean / t_l);
2781 float interm;
2782
2783 if(cp.edg_low < 10.f) {
2784 interm = cbrt((5.f + cp.edg_low) / 5.f);
2785 } else {
2786 interm = 1.437f; //cbrt(3);
2787 }
2788
2789 float edgeLowFactor = interm;
2790 float edgeSdFactor = cp.edg_sd / t_r;
2791 float edgeMaxFactor = SQR(cp.edg_max / b_r);
2792 float edgMaxFsup = (cp.edg_max / b_r); //reduce increase of effect for high values contrast..if slider > b_r
2793
2794 //for (int i=0; i<W_L*H_L; i++) {
2795 int borderL = 1;
2796
2797 for(int i = borderL; i < H_L - borderL; i++ ) {
2798 for(int j = borderL; j < W_L - borderL; j++) {
2799 int k = i * W_L + j;
2800
2801 if(cp.detectedge) {
2802 if(!lipschitz) {
2803 if(cp.eddet > 10.f) {
2804 edge = (aedstr * cp.eddet + bedstr) * (edgePrecalc * (1.f + koe[k])) / (1.f + 0.9f * maxkoe);
2805 } else {
2806 edge = (edgePrecalc * (1.f + koe[k])) / (1.f + 0.9f * maxkoe);
2807 }
2808 }
2809
2810 if(lipschitz) {
2811 if(level < 4) {
2812 edge = 1.f + (edgePrecalc - 1.f) * (koeLi[level * 3][k]) / (1.f + 0.9f * maxkoeLi[level * 3 + dir - 1]);
2813 } else {
2814 edge = edgePrecalc;
2815 }
2816 }
2817 } else {
2818 edge = edgePrecalc;
2819 }
2820
2821 //algorithm that takes into account local contrast
2822 // I use a thresholdadjuster with
2823 // Bottom left ==> minimal low value for local contrast (not 0, but 5...we can change)
2824 // 0 10*10 35*35 100*100 substantially correspond to the true distribution of low value, mean, standard-deviation and max (ed 5, 50, 400, 4000
2825 // Top left ==> mean reference value (for each level), we can change cbrt(cp.edg_mean/10.f)
2826 // Top Right==> standard deviation (for each level) we can change (cp.edg_sd/35.f)
2827 // bottom right ==> Max for positif and negatif contrast we can change cp.edg_max/100.f
2828 // If we move sliders to the left, local contrast is reduced
2829 // if we move sliders to the right local contrast is increased
2830 // MaxP, MaxN, mean, sigma are calculated if necessary (val > 0) by evaluate2(), eval2(), aver() , sigma()
2831 if(b_r < 100.f && cp.edg_max / b_r > 1.f) { //in case of b_r < 100 and slider move to right
2832 if (WavCoeffs_L[dir][k] > MaxPCompare * cp.edg_max / b_r) {
2833 edge *= edgMaxFsup;
2834
2835 if(edge < 1.f) {
2836 edge = 1.f;
2837 }
2838 } else if (WavCoeffs_L[dir][k] < MaxNCompare * cp.edg_max / b_r) {
2839 edge *= edgMaxFsup;
2840
2841 if(edge < 1.f) {
2842 edge = 1.f;
2843 }
2844 }
2845 }
2846
2847 if (WavCoeffs_L[dir][k] > MaxPCompare) {
2848 edge *= edgeMaxFactor;
2849
2850 if(edge < 1.f) {
2851 edge = 1.f;
2852 }
2853 }//reduce edge if > new max
2854 else if (WavCoeffs_L[dir][k] < MaxNCompare) {
2855 edge *= edgeMaxFactor;
2856
2857 if(edge < 1.f) {
2858 edge = 1.f;
2859 }
2860 }
2861
2862 if (fabs(WavCoeffs_L[dir][k]) >= edgeMeanCompare && fabs(WavCoeffs_L[dir][k]) < edgeSdCompare) {
2863 //if (fabs(WavCoeffs_L[dir][i]) > edgeSdCompare) {
2864 edge *= edgeSdFactor;
2865
2866 if(edge < 1.f) {
2867 edge = 1.f;
2868 }
2869 }//modify effect if sd change
2870
2871 if (fabs(WavCoeffs_L[dir][k]) < edgeMeanCompare) {
2872 edge *= edgeMeanFactor;
2873
2874 if(edge < 1.f) {
2875 edge = 1.f;
2876 }
2877 } // modify effect if mean change
2878
2879 if (fabs(WavCoeffs_L[dir][k]) < edgeLowCompare) {
2880 edge *= edgeLowFactor;
2881
2882 if(edge < 1.f) {
2883 edge = 1.f;
2884 }
2885 }
2886
2887 WavCoeffs_L[dir][k] *= edge;
2888 }
2889 }
2890 }
2891
2892 if(!lipschitz) {
2893 delete [] koe;
2894 }
2895 }
2896
2897
2898 if(!cp.link && cp.noiseena) { //used both with denoise 1 2 3
2899 float refine = 0.f;
2900
2901 for (int i = 0; i < W_L * H_L; i++) {
2902 if(level == 0) {
2903 refine = cp.lev0s / 40.f;
2904 }
2905
2906 if(level == 1) {
2907 refine = cp.lev1s / 40.f;
2908 }
2909
2910 if(level == 2) {
2911 refine = cp.lev2s / 40.f;
2912 }
2913
2914 if(level == 3) {
2915 refine = cp.lev3s / 40.f;
2916 }
2917
2918 WavCoeffs_L[dir][i] *= (1.f + refine);
2919 }
2920 }
2921
2922
2923 float cpMul = cp.mul[level];
2924
2925 if(cpMul != 0.f && cp.contena) { // cpMul == 0.f means all will be multiplied by 1.f, so we can skip this
2926
2927 const float skinprot = params->wavelet.skinprotect;
2928 const float skinprotneg = -skinprot;
2929 const float factorHard = (1.f - skinprotneg / 100.f);
2930
2931 //to adjust increase contrast with local contrast
2932
2933 //for each pixel and each level
2934 float beta;
2935 float mea[9];
2936 mea[0] = mean[level] / 6.f;
2937 mea[1] = mean[level] / 2.f;
2938 mea[2] = mean[level]; // 50% data
2939 mea[3] = mean[level] + sigma[level] / 2.f;
2940 mea[4] = mean[level] + sigma[level]; //66%
2941 mea[5] = mean[level] + 1.2f * sigma[level];
2942 mea[6] = mean[level] + 1.5f * sigma[level]; //
2943 mea[7] = mean[level] + 2.f * sigma[level]; //95%
2944 mea[8] = mean[level] + 2.5f * sigma[level]; //99%
2945
2946 bool useChromAndHue = (skinprot != 0.f || cp.HSmet);
2947 float modchro;
2948
2949 for (int i = 0; i < W_L * H_L; i++) {
2950 float kLlev = 1.f;
2951
2952 if(cpMul < 0.f) {
2953 beta = 1.f; // disabled for negatives values "less contrast"
2954 } else {
2955 float WavCL = fabsf(WavCoeffs_L[dir][i]);
2956
2957 //reduction amplification: max action between mean / 2 and mean + sigma
2958 // arbitrary coefficient, we can add a slider !!
2959 if(WavCL < mea[0]) {
2960 beta = 0.6f; //preserve very low contrast (sky...)
2961 } else if(WavCL < mea[1]) {
2962 beta = 0.8f;
2963 } else if(WavCL < mea[2]) {
2964 beta = 1.f; //standard
2965 } else if(WavCL < mea[3]) {
2966 beta = 1.f;
2967 } else if(WavCL < mea[4]) {
2968 beta = 0.8f; //+sigma
2969 } else if(WavCL < mea[5]) {
2970 beta = 0.6f;
2971 } else if(WavCL < mea[6]) {
2972 beta = 0.4f;
2973 } else if(WavCL < mea[7]) {
2974 beta = 0.2f; // + 2 sigma
2975 } else if(WavCL < mea[8]) {
2976 beta = 0.1f;
2977 } else {
2978 beta = 0.0f;
2979 }
2980 }
2981
2982 float scale = 1.f;
2983 float scale2 = 1.f;
2984
2985 float LL100, LL100res, LL100init, kH[maxlvl];
2986
2987 int ii = i / W_L;
2988 int jj = i - ii * W_L;
2989 float LL = labco->L[ii * 2][jj * 2];
2990 LL100 = LL100init = LL / 327.68f;
2991 LL100res = WavCoeffs_L0[i] / 327.68f;
2992 float delta = fabs(LL100init - LL100res) / (maxlvl / 2);
2993
2994 for(int ml = 0; ml < maxlvl; ml++) {
2995 if(ml < maxlvl / 2) {
2996 kH[ml] = (LL100res + ml * delta) / LL100res; // fixed a priori max to level middle
2997 } else {
2998 kH[ml] = (LL100init - ml * delta) / LL100res;
2999 }
3000 }
3001
3002
3003 if(useChromAndHue) {
3004 float modhue = varhue[ii][jj];
3005 modchro = varchrom[ii * 2][jj * 2];
3006 // hue chroma skin with initial lab data
3007 scale = 1.f;
3008
3009 if(skinprot > 0.f) {
3010 Color::SkinSatCbdl2 (LL100, modhue, modchro, skinprot, scale, true, cp.b_l, cp.t_l, cp.t_r, cp.b_r, 0); //0 for skin and extand
3011 } else if(skinprot < 0.f) {
3012 Color::SkinSatCbdl2 (LL100, modhue, modchro, skinprotneg, scale, false, cp.b_l, cp.t_l, cp.t_r, cp.b_r, 0);
3013
3014 if (scale == 1.f) {
3015 scale = factorHard;
3016 } else {
3017 scale = 1.f;
3018 }
3019 }
3020
3021 }
3022
3023 if(Chutili) {
3024 int i_i = i / W_L;
3025 int j_j = i - i_i * W_L;
3026 double lr;
3027 float modhue2 = varhue[i_i][j_j];
3028 float valparam = float((ChCurve->getVal(lr = Color::huelab_to_huehsv2(modhue2)) - 0.5f)); //get valparam=f(H)
3029
3030 if(valparam > 0.f) {
3031 scale2 = 1.f + 3.f * valparam; //arbitrary value
3032 } else {
3033 scale2 = 1.f + 1.9f * valparam; //near 0 but not zero if curve # 0
3034 }
3035 }
3036
3037 //linear transition HL
3038 float diagacc = 1.f;
3039 float alpha = (1024.f + 15.f * (float) cpMul * scale * scale2 * beta * diagacc) / 1024.f ;
3040
3041 if(cp.HSmet && cp.contena) {
3042 float aaal = (1.f - alpha) / ((cp.b_lhl - cp.t_lhl) * kH[level]);
3043 float bbal = 1.f - aaal * cp.b_lhl * kH[level];
3044 float aaar = (alpha - 1.f) / (cp.t_rhl - cp.b_rhl) * kH[level];
3045 float bbbr = 1.f - cp.b_rhl * aaar * kH[level];
3046 //linear transition Shadows
3047 float aaalS = (1.f - alpha) / (cp.b_lsl - cp.t_lsl);
3048 float bbalS = 1.f - aaalS * cp.b_lsl;
3049 float aaarS = (alpha - 1.f) / (cp.t_rsl - cp.b_rsl);
3050 float bbbrS = 1.f - cp.b_rsl * aaarS;
3051
3052 if(level <= cp.numlevH) { //in function of levels
3053 if((LL100 > cp.t_lhl * kH[level] && LL100 < cp.t_rhl * kH[level])) {
3054 kLlev = alpha;
3055 } else if((LL100 > cp.b_lhl * kH[level] && LL100 <= cp.t_lhl * kH[level])) {
3056 kLlev = aaal * LL100 + bbal;
3057 } else if((LL100 > cp.t_rhl * kH[level] && LL100 <= cp.b_rhl * kH[level])) {
3058 kLlev = aaar * LL100 + bbbr;
3059 } else {
3060 kLlev = 1.f;
3061 }
3062 }
3063
3064 if(level >= (9 - cp.numlevS)) {
3065 if((LL100 > cp.t_lsl && LL100 < cp.t_rsl)) {
3066 kLlev = alpha;
3067 } else if((LL100 > cp.b_lsl && LL100 <= cp.t_lsl)) {
3068 kLlev = aaalS * LL100 + bbalS;
3069 } else if((LL100 > cp.t_rsl && LL100 <= cp.b_rsl)) {
3070 kLlev = aaarS * LL100 + bbbrS;
3071 } else {
3072 kLlev = 1.f;
3073 }
3074 }
3075
3076 } else {
3077 kLlev = alpha;
3078 }
3079
3080 WavCoeffs_L[dir][i] *= (kLlev);
3081 }
3082 }
3083
3084 if(waOpacityCurveW) {
3085 cp.opaW = true;
3086 }
3087
3088 if(cp.bam && cp.finena) {
3089 if(cp.opaW && cp.BAmet == 2) {
3090 int iteration = cp.ite;
3091 int itplus = 7 + iteration;
3092 int itmoins = 7 - iteration;
3093 int med = maxlvl / 2;
3094 int it;
3095
3096 if(level < med) {
3097 it = itmoins;
3098 } else if(level == med) {
3099 it = 7;
3100 } else /*if(level > med)*/ {
3101 it = itplus;
3102 }
3103
3104 for(int j = 0; j < it; j++) {
3105 //float bal = cp.balan;//-100 +100
3106 float kba = 1.f;
3107
3108 // if(dir <3) kba= 1.f + bal/600.f;
3109 // if(dir==3) kba = 1.f - bal/300.f;
3110 for (int i = 0; i < W_L * H_L; i++) {
3111 int ii = i / W_L;
3112 int jj = i - ii * W_L;
3113 float LL100 = labco->L[ii * 2][jj * 2] / 327.68f;
3114 float k1 = 0.3f * (waOpacityCurveW[6.f * LL100] - 0.5f); //k1 between 0 and 0.5 0.5==> 1/6=0.16
3115 float k2 = k1 * 2.f;
3116
3117 if(dir < 3) {
3118 kba = 1.f + k1;
3119 }
3120
3121 if(dir == 3) {
3122 kba = 1.f - k2;
3123 }
3124
3125 WavCoeffs_L[dir][i] *= (kba);
3126 }
3127 }
3128 }
3129
3130 if(cp.BAmet == 1) {
3131 int iteration = cp.ite;
3132 int itplus = 7 + iteration;
3133 int itmoins = 7 - iteration;
3134 int med = maxlvl / 2;
3135 int it;
3136
3137 if(level < med) {
3138 it = itmoins;
3139 } else if(level == med) {
3140 it = 7;
3141 } else /*if(level > med)*/ {
3142 it = itplus;
3143 }
3144
3145 for(int j = 0; j < it; j++) {
3146 float bal = cp.balan;//-100 +100
3147 float kba = 1.f;
3148
3149 // if(dir <3) kba= 1.f + bal/600.f;
3150 // if(dir==3) kba = 1.f - bal/300.f;
3151 for (int i = 0; i < W_L * H_L; i++) {
3152 int ii = i / W_L;
3153 int jj = i - ii * W_L;
3154 float k1 = 600.f;
3155 float k2 = 300.f;
3156 float LL100 = labco->L[ii * 2][jj * 2] / 327.68f;
3157 float aa = 4970.f;
3158 float bb = -397000.f;
3159 float b0 = 100000.f;
3160 float a0 = -4970.f;
3161
3162 if(LL100 > 80.f) {
3163 k1 = aa * LL100 + bb;
3164 k2 = 0.5f * k1;
3165 }
3166
3167 if(LL100 < 20.f) {
3168 k1 = a0 * LL100 + b0;
3169 k2 = 0.5f * k1;
3170 }
3171
3172 //k1=600.f;
3173 //k2=300.f;
3174 //k1=0.3f*(waOpacityCurveW[6.f*LL100]-0.5f);//k1 between 0 and 0.5 0.5==> 1/6=0.16
3175 //k2=k1*2.f;
3176 if(dir < 3) {
3177 kba = 1.f + bal / k1;
3178 }
3179
3180 if(dir == 3) {
3181 kba = 1.f - bal / k2;
3182 }
3183
3184 WavCoeffs_L[dir][i] *= (kba);
3185 }
3186 }
3187 }
3188
3189 }
3190
3191 // to see each level of wavelet ...level from 0 to 8
3192 int choicelevel = params->wavelet.Lmethod - 1;
3193 choicelevel = choicelevel == -1 ? 4 : choicelevel;
3194 }
3195
ContAllAB(LabImage * labco,int maxlvl,float ** varhue,float ** varchrom,float ** WavCoeffs_ab,float * WavCoeffs_ab0,int level,int dir,const WavOpacityCurveW & waOpacityCurveW,struct cont_params & cp,int W_ab,int H_ab,const bool useChannelA)3196 void ImProcFunctions::ContAllAB (LabImage * labco, int maxlvl, float ** varhue, float **varchrom, float ** WavCoeffs_ab, float * WavCoeffs_ab0, int level, int dir, const WavOpacityCurveW & waOpacityCurveW, struct cont_params &cp,
3197 int W_ab, int H_ab, const bool useChannelA)
3198 {
3199 float cpMul = cp.mul[level];
3200
3201 if(cpMul != 0.f && cp.CHmet == 2 && cp.chro != 0.f && cp.chromena) { // cpMul == 0.f or cp.chro = 0.f means all will be multiplied by 1.f, so we can skip this
3202 const float skinprot = params->wavelet.skinprotect;
3203 const float skinprotneg = -skinprot;
3204 const float factorHard = (1.f - skinprotneg / 100.f);
3205 const float cpChrom = cp.chro;
3206
3207 //to adjust increase contrast with local contrast
3208 bool useSkinControl = (skinprot != 0.f);
3209 float alphaC = (1024.f + 15.f * cpMul * cpChrom / 50.f) / 1024.f ;
3210
3211 for (int i = 0; i < W_ab * H_ab; i++) {
3212 if(useSkinControl) {
3213 int ii = i / W_ab;
3214 int jj = i - ii * W_ab;
3215 float LL100 = labco->L[ii * 2][jj * 2] / 327.68f;
3216 float modhue = varhue[ii][jj];
3217 float modchro = varchrom[ii * 2][jj * 2];
3218 // hue chroma skin with initial lab data
3219 float scale = 1.f;
3220
3221 if(skinprot > 0.f) {
3222 Color::SkinSatCbdl2 (LL100, modhue, modchro, skinprot, scale, true, cp.b_l, cp.t_l, cp.t_r, cp.b_r, 0); //0 for skin and extand
3223 } else if(skinprot < 0.f) {
3224 Color::SkinSatCbdl2 (LL100, modhue, modchro, skinprotneg, scale, false, cp.b_l, cp.t_l, cp.t_r, cp.b_r, 0);
3225 scale = (scale == 1.f) ? factorHard : 1.f;
3226 }
3227
3228 alphaC = (1024.f + 15.f * cpMul * cpChrom * scale / 50.f) / 1024.f ;
3229 }
3230
3231 WavCoeffs_ab[dir][i] *= alphaC;
3232 }
3233 }
3234
3235 //Curve chro
3236
3237 float cpMulC = cp.mulC[level];
3238
3239 // if( (cp.curv || cp.CHSLmet==1) && cp.CHmet!=2 && level < 9 && cpMulC != 0.f) { // cpMulC == 0.f means all will be multiplied by 1.f, so we can skip
3240 if( cp.CHmet != 2 && level < 9 && cpMulC != 0.f && cp.chromena) { // cpMulC == 0.f means all will be multiplied by 1.f, so we can skip
3241 const float skinprot = params->wavelet.skinprotect;
3242 const float skinprotneg = -skinprot;
3243 const float factorHard = (1.f - skinprotneg / 100.f);
3244 bool useSkinControl = (skinprot != 0.f);
3245
3246 for (int i = 0; i < W_ab * H_ab; i++) {
3247 int ii = i / W_ab;
3248 int jj = i - ii * W_ab;
3249 //WL and W_ab are identical
3250 float scale = 1.f;
3251 float modchro = varchrom[ii * 2][jj * 2];
3252
3253 if(useSkinControl) {
3254 // hue chroma skin with initial lab data
3255 float LL100 = labco->L[ii * 2][jj * 2] / 327.68f;
3256 float modhue = varhue[ii][jj];
3257
3258 if(skinprot > 0.f) {
3259 Color::SkinSatCbdl2 (LL100, modhue, modchro, skinprot, scale, true, cp.b_l, cp.t_l, cp.t_r, cp.b_r, 1); //1 for curve
3260 } else if(skinprot < 0.f) {
3261 Color::SkinSatCbdl2 (LL100, modhue, modchro, skinprotneg, scale, false, cp.b_l, cp.t_l, cp.t_r, cp.b_r, 1);
3262 scale = (scale == 1.f) ? factorHard : 1.f;
3263 }
3264 }
3265
3266 float beta = (1024.f + 20.f * cpMulC * scale) / 1024.f ;
3267
3268 if(beta < 0.02f) {
3269 beta = 0.02f;
3270 }
3271
3272 float kClev = beta;
3273
3274 if(cp.CHmet == 1) {
3275 if(level < cp.chrom) {
3276 //linear for saturated
3277 if((modchro > cp.t_lsat && modchro < cp.t_rsat)) {
3278 kClev = beta;
3279 } else if((modchro > cp.b_lsat && modchro <= cp.t_lsat)) {
3280 float aaal = (1.f - beta) / (cp.b_lsat - cp.t_lsat);
3281 float bbal = 1.f - aaal * cp.b_lsat;
3282 kClev = aaal * modchro + bbal;
3283 } else if((modchro > cp.t_rsat && modchro <= cp.b_rsat)) {
3284 float aaar = (beta - 1.f) / (cp.t_rsat - cp.b_rsat);
3285 float bbbr = 1.f - cp.b_rsat * aaar;
3286 kClev = aaar * modchro + bbbr;
3287 } else {
3288 kClev = 1.f;
3289 }
3290 } else {
3291 //linear for pastel
3292 if((modchro > cp.t_lpast && modchro < cp.t_rpast)) {
3293 kClev = beta;
3294 } else if((modchro > cp.b_lpast && modchro <= cp.t_lpast)) {
3295 float aaalS = (1.f - beta) / (cp.b_lpast - cp.t_lpast);
3296 float bbalS = 1.f - aaalS * cp.b_lpast;
3297 kClev = aaalS * modchro + bbalS;
3298 } else if((modchro > cp.t_rpast && modchro <= cp.b_rpast)) {
3299 float aaarS = (beta - 1.f) / (cp.t_rpast - cp.b_rpast);
3300 float bbbrS = 1.f - cp.b_rpast * aaarS;
3301 kClev = aaarS * modchro + bbbrS;
3302 } else {
3303 kClev = 1.f;
3304 }
3305 }
3306 } else if(cp.CHmet == 0) {
3307 kClev = beta;
3308 }
3309
3310 WavCoeffs_ab[dir][i] *= kClev;
3311 }
3312 }
3313
3314 bool useOpacity;
3315 float mulOpacity = 0.f;
3316
3317 if(useChannelA) {
3318 useOpacity = cp.opaRG;
3319 if(level < 9) {
3320 mulOpacity = cp.mulopaRG[level];
3321 }
3322 } else {
3323 useOpacity = cp.opaBY;
3324 if(level < 9) {
3325 mulOpacity = cp.mulopaBY[level];
3326 }
3327 }
3328
3329 if((useOpacity && level < 9 && mulOpacity != 0.f) && cp.toningena) { //toning
3330
3331 float beta = (1024.f + 20.f * mulOpacity) / 1024.f ;
3332
3333 //float beta = (1000.f * mulOpacity);
3334 for (int i = 0; i < W_ab * H_ab; i++) {
3335 WavCoeffs_ab[dir][i] *= beta;
3336 }
3337
3338 // WavCoeffs_ab[dir][i] += beta;
3339 }
3340
3341 if(waOpacityCurveW) {
3342 cp.opaW = true;
3343 }
3344
3345 if(cp.bam && cp.diag) {
3346 //printf("OK Chroma\n");
3347 if(cp.opaW && cp.BAmet == 2) {
3348 int iteration = cp.ite;
3349 int itplus = 7 + iteration;
3350 int itmoins = 7 - iteration;
3351 int med = maxlvl / 2;
3352 int it;
3353
3354 if(level < med) {
3355 it = itmoins;
3356 } else if(level == med) {
3357 it = 7;
3358 } else /*if(level > med)*/ {
3359 it = itplus;
3360 }
3361
3362 for(int j = 0; j < it; j++) {
3363 //float bal = cp.balan;//-100 +100
3364 float kba = 1.f;
3365
3366 // if(dir <3) kba= 1.f + bal/600.f;
3367 // if(dir==3) kba = 1.f - bal/300.f;
3368 for (int i = 0; i < W_ab * H_ab; i++) {
3369 int ii = i / W_ab;
3370 int jj = i - ii * W_ab;
3371 float LL100 = labco->L[ii * 2][jj * 2] / 327.68f;
3372 float k1 = 0.3f * (waOpacityCurveW[6.f * LL100] - 0.5f); //k1 between 0 and 0.5 0.5==> 1/6=0.16
3373 float k2 = k1 * 2.f;
3374
3375 if(dir < 3) {
3376 kba = 1.f + k1;
3377 }
3378
3379 if(dir == 3) {
3380 kba = 1.f - k2;
3381 }
3382
3383 WavCoeffs_ab[dir][i] *= (kba);
3384 }
3385 }
3386 }
3387
3388 if(cp.BAmet == 1) {
3389 int iteration = cp.ite;
3390 int itplus = 7 + iteration;
3391 int itmoins = 7 - iteration;
3392 int med = maxlvl / 2;
3393 int it;
3394
3395 if(level < med) {
3396 it = itmoins;
3397 } else if(level == med) {
3398 it = 7;
3399 } else /*if(level > med)*/ {
3400 it = itplus;
3401 }
3402
3403 for(int j = 0; j < it; j++) {
3404 float bal = cp.balan;//-100 +100
3405 float kba = 1.f;
3406
3407 // if(dir <3) kba= 1.f + bal/600.f;
3408 // if(dir==3) kba = 1.f - bal/300.f;
3409 for (int i = 0; i < W_ab * H_ab; i++) {
3410 int ii = i / W_ab;
3411 int jj = i - ii * W_ab;
3412 float k1 = 600.f;
3413 float k2 = 300.f;
3414 float LL100 = labco->L[ii * 2][jj * 2] / 327.68f;
3415 float aa = 4970.f;
3416 float bb = -397000.f;
3417 float b0 = 100000.f;
3418 float a0 = -4970.f;
3419
3420 if(LL100 > 80.f) {
3421 k1 = aa * LL100 + bb;
3422 k2 = 0.5f * k1;
3423 }
3424
3425 if(LL100 < 20.f) {
3426 k1 = a0 * LL100 + b0;
3427 k2 = 0.5f * k1;
3428 }
3429
3430 //k1=600.f;
3431 //k2=300.f;
3432 //k1=0.3f*(waOpacityCurveW[6.f*LL100]-0.5f);//k1 between 0 and 0.5 0.5==> 1/6=0.16
3433 //k2=k1*2.f;
3434 if(dir < 3) {
3435 kba = 1.f + bal / k1;
3436 }
3437
3438 if(dir == 3) {
3439 kba = 1.f - bal / k2;
3440 }
3441
3442 WavCoeffs_ab[dir][i] *= (kba);
3443 }
3444 }
3445 }
3446
3447 }
3448
3449 // to see each level of wavelet ...level from 0 to 8
3450 int choicelevel = params->wavelet.Lmethod - 1;
3451 choicelevel = choicelevel == -1 ? 4 : choicelevel;
3452 int choiceClevel = 0;
3453
3454 if(params->wavelet.CLmethod == "one") {
3455 choiceClevel = 0;
3456 } else if(params->wavelet.CLmethod == "inf") {
3457 choiceClevel = 1;
3458 } else if(params->wavelet.CLmethod == "sup") {
3459 choiceClevel = 2;
3460 } else if(params->wavelet.CLmethod == "all") {
3461 choiceClevel = 3;
3462 }
3463
3464 int choiceDir = 0;
3465
3466 if(params->wavelet.Dirmethod == "one") {
3467 choiceDir = 1;
3468 } else if(params->wavelet.Dirmethod == "two") {
3469 choiceDir = 2;
3470 } else if(params->wavelet.Dirmethod == "thr") {
3471 choiceDir = 3;
3472 } else if(params->wavelet.Dirmethod == "all") {
3473 choiceDir = 0;
3474 }
3475
3476 int dir1 = (choiceDir == 2) ? 1 : 2;
3477 int dir2 = (choiceDir == 3) ? 1 : 3;
3478
3479 if(choiceClevel < 3) { // not all levels visible, paint residual
3480 if(level == 0) {
3481 if(cp.backm != 2) { // nothing to change when residual is used as background
3482 for (int i = 0; i < W_ab * H_ab; i++) {
3483 WavCoeffs_ab0[i] = 0.f;
3484 }
3485 }
3486 }
3487 }
3488
3489 if(choiceClevel == 0) { // Only one level
3490 if(choiceDir == 0) { // All directions
3491 if(level != choicelevel) { // zero all for the levels != choicelevel
3492 for (int dir = 1; dir < 4; dir++) {
3493 for (int i = 0; i < W_ab * H_ab; i++) {
3494 WavCoeffs_ab[dir][i] = 0.f;
3495 }
3496 }
3497 }
3498 } else { // zero the unwanted directions for level == choicelevel
3499 if(choicelevel >= cp.maxilev) {
3500 for (int dir = 1; dir < 4; dir++) {
3501 for (int i = 0; i < W_ab * H_ab; i++) {
3502 WavCoeffs_ab[dir][i] = 0.f;
3503 }
3504 }
3505 } else if(level != choicelevel) { // zero all for the levels != choicelevel
3506 for (int i = 0; i < W_ab * H_ab; i++) {
3507 WavCoeffs_ab[dir1][i] = WavCoeffs_ab[dir2][i] = 0.f;
3508 }
3509 }
3510 }
3511 } else if(choiceClevel == 1) { // Only below level
3512 if(choiceDir == 0) { // All directions
3513 if(level > choicelevel) {
3514 for (int dir = 1; dir < 4; dir++) {
3515 for (int i = 0; i < W_ab * H_ab; i++) {
3516 WavCoeffs_ab[dir][i] = 0.f;
3517 }
3518 }
3519 }
3520 } else { // zero the unwanted directions for level >= choicelevel
3521 if(level > choicelevel) {
3522 for (int i = 0; i < W_ab * H_ab; i++) {
3523 WavCoeffs_ab[dir1][i] = WavCoeffs_ab[dir2][i] = 0.f;
3524 }
3525 }
3526 }
3527 } else if(choiceClevel == 2) { // Only above level
3528 if(choiceDir == 0) { // All directions
3529 if(level <= choicelevel) {
3530 for (int dir = 1; dir < 4; dir++) {
3531 for (int i = 0; i < W_ab * H_ab; i++) {
3532 WavCoeffs_ab[dir][i] = 0.f;
3533 }
3534 }
3535 }
3536 } else { // zero the unwanted directions for level >= choicelevel
3537 if(choicelevel >= cp.maxilev) {
3538 for (int dir = 1; dir < 4; dir++) {
3539 for (int i = 0; i < W_ab * H_ab; i++) {
3540 WavCoeffs_ab[dir][i] = 0.f;
3541 }
3542 }
3543 } else if(level <= choicelevel) {
3544 for (int i = 0; i < W_ab * H_ab; i++) {
3545 WavCoeffs_ab[dir1][i] = WavCoeffs_ab[dir2][i] = 0.f;
3546 }
3547 }
3548 }
3549 }
3550 }
3551 }
3552