1 /*******************************************************************************
2 **
3 ** Photivo
4 **
5 ** Copyright (C) 2008,2009 Jos De Laender <jos.de_laender@telenet.be>
6 ** Copyright (C) 2009-2011 Michael Munzert <mail@mm-log.com>
7 **
8 ** This file is part of Photivo.
9 **
10 ** Photivo is free software: you can redistribute it and/or modify
11 ** it under the terms of the GNU General Public License version 3
12 ** as published by the Free Software Foundation.
13 **
14 ** Photivo is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 ** GNU General Public License for more details.
18 **
19 ** You should have received a copy of the GNU General Public License
20 ** along with Photivo. If not, see <http://www.gnu.org/licenses/>.
21 **
22 *******************************************************************************/
23
24 #include "ptDcRaw.h"
25
26 #include "ptCalloc.h"
27 #include "ptConstants.h"
28 #include "ptError.h"
29 #include "ptInfo.h"
30 #include "ptImage.h"
31 #include "ptMessageBox.h"
32 #include "ptResizeFilters.h"
33 #include "ptCurve.h"
34 #include "ptKernel.h"
35 #include "ptConstants.h"
36 #include "ptRefocusMatrix.h"
37 #include "ptCimg.h"
38 #include "fastbilateral/fast_lbf.h"
39
40 #include <QString>
41 #include <QTime>
42
43 #ifdef __GLIBCXX__
44 #include <parallel/algorithm>
45 #else
46 #include <algorithm>
47 #endif
48 #include <stack>
49
50 // std stuff needs to be declared apparently for jpeglib
51 // which seems a bug in the jpeglib header ?
52 #include <cstdlib>
53 #include <cstdio>
54
55 #include <functional>
56
57 #ifdef _OPENMP
58 #include <omp.h>
59 #endif
60
61 #ifdef WIN32
62 #define NO_JPEG
63 #endif
64
65 #ifndef NO_JPEG
66 #ifdef __cplusplus
67 // This hack copes with jpeglib.h that does or doesnt provide the
68 // extern internally.
69 #pragma push_macro("__cplusplus")
70 #undef __cplusplus
71 extern "C" {
72 #include <jpeglib.h>
73 }
74 #pragma pop_macro("__cplusplus")
75 #else
76 #include <jpeglib.h>
77 #endif
78 #endif
79
80 extern cmsCIExyY D65;
81 extern cmsCIExyY D50;
82
83 extern cmsHPROFILE PreviewColorProfile;
84 extern cmsHTRANSFORM ToPreviewTransform;
85
86 // Lut
87 extern float ToFloatTable[0x10000];
88 extern float ToFloatABNeutral[0x10000];
89 extern uint16_t ToInvertTable[0x10000];
90 extern uint16_t ToSRGBTable[0x10000];
91
92 short ptImage::CurrentRGBMode = 0;
93
94 //==============================================================================
95
96 //Helping functions
97
98 //==============================================================================
99
100 // Convert an RGB pixel to fake luminance
RGB_2_L(const uint16_t APixel[3])101 inline uint16_t RGB_2_L(const uint16_t APixel[3]) {
102 return CLIP((int32_t) (0.30f*APixel[0] +
103 0.59f*APixel[1] +
104 0.11f*APixel[2]));
105 }
106
107 //==============================================================================
108
109 // calculate the hue for A and B differences
ToHue(const float ADiffA,const float ADiffB)110 inline float ToHue(const float ADiffA, const float ADiffB) {
111 if (ADiffA == 0.0f && ADiffB == 0.0f) {
112 return 0.0f; // value for grey pixel
113 } else {
114 float hHue = atan2f(ADiffB,ADiffA);
115 while (hHue < 0.0f) hHue += pt2PI;
116 return hHue;
117 }
118 }
119
120 //==============================================================================
121
Sigmoidal_4_Value(const uint16_t AValue,const float APosContrast)122 inline uint16_t Sigmoidal_4_Value(const uint16_t AValue, const float APosContrast) {
123 float hContrastHalfExp = exp(0.5f*APosContrast);
124 float hOffset = -1.0f/(1.0f + hContrastHalfExp);
125 float hScaling = 1.0f/(1.0f + 1.0f/hContrastHalfExp) + hOffset;
126 return CLIP((int32_t)((((1.0f/(1.0f + exp(APosContrast*(0.5f - ToFloatTable[AValue])))) + hOffset)/hScaling)*ptWPf));
127 }
128
129 //==============================================================================
130
SigmoidalTable(uint16_t (& ATable)[0x10000],const float AContrast,const float AThreshold)131 void SigmoidalTable(uint16_t (&ATable)[0x10000], const float AContrast, const float AThreshold) {
132 float hInvContrast = 1.0f/AContrast;
133 float hContrastHalfExp = exp(0.5f*AContrast);
134 float hOffset = -1.0f/(1.0f + hContrastHalfExp);
135 float hScaling = 1.0f/(1.0f + 1.0f/hContrastHalfExp) + hOffset;
136 float hInvScaling = 1.0f/hScaling;
137 float logtf = -logf(AThreshold)/logf(2.0f);
138 float logft = -logf(2.0f)/logf(AThreshold);
139
140 ATable[0] = 0;
141 if (AContrast > 0) {
142 #pragma omp parallel for
143 for (uint32_t i = 1; i < 0x10000; i++) {
144 ATable[i] = CLIP((int32_t)(powf((((1.0f/(1.0f +
145 expf(AContrast*(0.5f - powf(ToFloatTable[i],logft))))) + hOffset)*hInvScaling),logtf)*ptWPf));
146 }
147 } else {
148 #pragma omp parallel for
149 for (uint32_t i = 1; i < 0x10000; i++) {
150 ATable[i] = CLIP((int32_t)(powf(0.5f - hInvContrast*
151 logf(1.0f/(hScaling*powf(ToFloatTable[i],logft) - hOffset) - 1.0f),logtf)*ptWPf));
152 }
153 }
154 }
155
156 ////////////////////////////////////////////////////////////////////////////////
157 //
158 // RGBToRGB.
159 // Conversion from RGB space to RGB space.
160 //
161 ////////////////////////////////////////////////////////////////////////////////
162
RGBToRGB(const short To,const short EvenIfEqual)163 ptImage* ptImage::RGBToRGB(const short To,
164 const short EvenIfEqual) {
165
166 assert ((m_ColorSpace>0) && (m_ColorSpace<5));
167 assert ((To>0) && (To<5));
168 assert (3 == m_Colors);
169
170 if ((m_ColorSpace == To) && !EvenIfEqual) return this;
171
172 // Matrix for conversion is a multiplication via XYZ
173 double Matrix[3][3];
174
175 for(short i=0;i<3;i++) {
176 for (short j=0;j<3;j++) {
177 Matrix[i][j] = 0.0;
178 for (short k=0;k<3;k++) {
179 Matrix[i][j] +=
180 MatrixXYZToRGB[To][i][k] * MatrixRGBToXYZ[m_ColorSpace][k][j];
181 }
182 }
183 }
184
185 // Convert the image.
186 #pragma omp parallel for
187 for (uint32_t i=0; i<(uint32_t)m_Height*m_Width; i++) {
188 int32_t Value[3];
189 for (short c=0; c<3; c++) {
190 Value[c] = 0;
191 for (short k=0; k<3; k++) {
192 Value[c] += (uint32_t) (Matrix[c][k]*m_Image[i][k]);
193 }
194 }
195 for (short c=0; c<3; c++) {
196 m_Image[i][c] = (uint16_t) CLIP(Value[c]);
197 }
198 }
199
200 m_ColorSpace = To;
201
202 return this;
203 }
204
205 ////////////////////////////////////////////////////////////////////////////////
206 //
207 // lcmsRGBToRGB.
208 // Conversion from RGB space to RGB space using lcms.
209 //
210 ////////////////////////////////////////////////////////////////////////////////
211
lcmsRGBToRGB(const short To,const short EvenIfEqual,const int Intent)212 ptImage* ptImage::lcmsRGBToRGB(const short To,
213 const short EvenIfEqual,
214 const int Intent) {
215
216 // assert ((m_ColorSpace>0) && (m_ColorSpace<5));
217 if (!((m_ColorSpace>0) && (m_ColorSpace<5))) {
218 ptMessageBox::critical(0,"Error","Too fast! Keep cool ;-)");
219 return this;
220 }
221 assert ((To>0) && (To<5));
222 assert (3 == m_Colors);
223
224 if ((m_ColorSpace == To) && !EvenIfEqual) return this;
225
226 cmsToneCurve* Gamma = cmsBuildGamma(NULL, 1.0);
227 cmsToneCurve* Gamma3[3];
228 Gamma3[0] = Gamma3[1] = Gamma3[2] = Gamma;
229
230 cmsHPROFILE InProfile = 0;
231
232 cmsCIExyY DFromReference;
233
234 switch (m_ColorSpace) {
235 case ptSpace_sRGB_D65 :
236 case ptSpace_AdobeRGB_D65 :
237 DFromReference = D65;
238 break;
239 case ptSpace_WideGamutRGB_D50 :
240 case ptSpace_ProPhotoRGB_D50 :
241 DFromReference = D50;
242 break;
243 default:
244 assert(0);
245 }
246
247 InProfile = cmsCreateRGBProfile(&DFromReference,
248 (cmsCIExyYTRIPLE*)&RGBPrimaries[m_ColorSpace],
249 Gamma3);
250
251 if (!InProfile) {
252 ptLogError(ptError_Profile,"Could not open InProfile profile.");
253 return NULL;
254 }
255
256 cmsHPROFILE OutProfile = 0;
257
258 cmsCIExyY DToReference;
259
260 switch (To) {
261 case ptSpace_sRGB_D65 :
262 case ptSpace_AdobeRGB_D65 :
263 DToReference = D65;
264 break;
265 case ptSpace_WideGamutRGB_D50 :
266 case ptSpace_ProPhotoRGB_D50 :
267 DToReference = D50;
268 break;
269 default:
270 assert(0);
271 }
272
273 OutProfile = cmsCreateRGBProfile(&DToReference,
274 (cmsCIExyYTRIPLE*)&RGBPrimaries[To],
275 Gamma3);
276
277 if (!OutProfile) {
278 ptLogError(ptError_Profile,"Could not open OutProfile profile.");
279 cmsCloseProfile (InProfile);
280 return NULL;
281 }
282
283 cmsFreeToneCurve(Gamma);
284
285 cmsHTRANSFORM Transform;
286 Transform = cmsCreateTransform(InProfile,
287 TYPE_RGB_16,
288 OutProfile,
289 TYPE_RGB_16,
290 Intent,
291 cmsFLAGS_NOOPTIMIZE | cmsFLAGS_BLACKPOINTCOMPENSATION);
292
293 int32_t Size = m_Width*m_Height;
294 int32_t Step = 100000;
295 #pragma omp parallel for schedule(static)
296 for (int32_t i = 0; i < Size; i+=Step) {
297 int32_t Length = (i+Step)<Size ? Step : Size - i;
298 uint16_t* Tile = &(m_Image[i][0]);
299 cmsDoTransform(Transform,Tile,Tile,Length);
300 }
301 cmsDeleteTransform(Transform);
302 cmsCloseProfile(OutProfile);
303 cmsCloseProfile(InProfile);
304
305 m_ColorSpace = To;
306
307 return this;
308 }
309
310 ////////////////////////////////////////////////////////////////////////////////
311 //
312 // lcmsRGBToRGB.
313 // Conversion from RGB space to RGB space using lcms.
314 // This variant uses an external profile for its output transformation.
315 // It would be used for instance to map onto an physical display.
316 //
317 ////////////////////////////////////////////////////////////////////////////////
318
lcmsRGBToRGB(cmsHPROFILE OutProfile,const int Intent,const short Quality)319 ptImage* ptImage::lcmsRGBToRGB(cmsHPROFILE OutProfile,
320 const int Intent,
321 const short Quality){
322
323 // assert ((m_ColorSpace>0) && (m_ColorSpace<5));
324 if (!((m_ColorSpace>0) && (m_ColorSpace<5))) {
325 ptMessageBox::critical(0,"Error","Too fast! Keep cool ;-)");
326 return this;
327 }
328 assert (3 == m_Colors);
329 assert (OutProfile);
330
331 cmsToneCurve* Gamma = cmsBuildGamma(NULL, 1.0);
332 cmsToneCurve* Gamma3[3];
333 Gamma3[0] = Gamma3[1] = Gamma3[2] = Gamma;
334
335 cmsHPROFILE InProfile = 0;
336
337 cmsCIExyY DFromReference;
338
339 switch (m_ColorSpace) {
340 case ptSpace_sRGB_D65 :
341 case ptSpace_AdobeRGB_D65 :
342 DFromReference = D65;
343 break;
344 case ptSpace_WideGamutRGB_D50 :
345 case ptSpace_ProPhotoRGB_D50 :
346 DFromReference = D50;
347 break;
348 default:
349 assert(0);
350 }
351
352 InProfile = cmsCreateRGBProfile(&DFromReference,
353 (cmsCIExyYTRIPLE*)&RGBPrimaries[m_ColorSpace],
354 Gamma3);
355
356 if (!InProfile) {
357 ptLogError(ptError_Profile,"Could not open InProfile profile.");
358 return NULL;
359 }
360
361 cmsFreeToneCurve(Gamma);
362
363 cmsHTRANSFORM Transform;
364 if (Quality == ptCMQuality_HighResPreCalc) {
365 Transform =
366 cmsCreateTransform(InProfile,
367 TYPE_RGB_16,
368 OutProfile,
369 TYPE_RGB_16,
370 Intent,
371 cmsFLAGS_HIGHRESPRECALC | cmsFLAGS_BLACKPOINTCOMPENSATION);
372 } else { // fast sRGB preview also uses the not optimized profile for output
373 Transform =
374 cmsCreateTransform(InProfile,
375 TYPE_RGB_16,
376 OutProfile,
377 TYPE_RGB_16,
378 Intent,
379 cmsFLAGS_NOOPTIMIZE | cmsFLAGS_BLACKPOINTCOMPENSATION);
380 }
381
382 int32_t Size = m_Width*m_Height;
383 int32_t Step = 10000;
384 #pragma omp parallel for schedule(static)
385 for (int32_t i = 0; i < Size; i+=Step) {
386 int32_t Length = (i+Step)<Size ? Step : Size - i;
387 uint16_t* Image = &m_Image[i][0];
388 cmsDoTransform(Transform,Image,Image,Length);
389 }
390
391 cmsDeleteTransform(Transform);
392 cmsCloseProfile(InProfile);
393
394 m_ColorSpace = ptSpace_Profiled;
395
396 return this;
397 }
398
lcmsRGBToPreviewRGB(const bool Fast)399 ptImage* ptImage::lcmsRGBToPreviewRGB(const bool Fast /* = false */){
400 if (Fast) {
401 // Fast transform to sRGB by lookup
402 RGBToRGB(ptSpace_sRGB_D65);
403 #pragma omp parallel for schedule(static)
404 for (uint32_t i=0; i<(uint32_t)m_Height*m_Width; i++) {
405 m_Image[i][0] = ToSRGBTable[m_Image[i][0]];
406 m_Image[i][1] = ToSRGBTable[m_Image[i][1]];
407 m_Image[i][2] = ToSRGBTable[m_Image[i][2]];
408 }
409 } else {
410 // assert ((m_ColorSpace>0) && (m_ColorSpace<5));
411 if (!((m_ColorSpace>0) && (m_ColorSpace<5))) {
412 ptMessageBox::critical(0,"Error","Too fast! Keep cool ;-)");
413 return this;
414 }
415 assert (3 == m_Colors);
416
417 int32_t Size = m_Width*m_Height;
418 int32_t Step = 10000;
419 #pragma omp parallel for schedule(static)
420 for (int32_t i = 0; i < Size; i+=Step) {
421 int32_t Length = (i+Step)<Size ? Step : Size - i;
422 uint16_t* Image = &m_Image[i][0];
423 cmsDoTransform(ToPreviewTransform,Image,Image,Length);
424 }
425 }
426 m_ColorSpace = ptSpace_Profiled;
427
428 return this;
429 }
430
431 ////////////////////////////////////////////////////////////////////////////////
432 //
433 // RGBToXYZ.
434 // Conversion from RGB space to XYZ space.
435 //
436 ////////////////////////////////////////////////////////////////////////////////
437
RGBToXYZ()438 ptImage* ptImage::RGBToXYZ() {
439
440 assert ((m_ColorSpace>0) && (m_ColorSpace<5));
441
442 // Convert the image.
443 #pragma omp parallel for
444 for (uint32_t i=0; i<(uint32_t)m_Height*m_Width; i++) {
445 int32_t Value[3];
446 for (short c=0; c<3; c++) {
447 Value[c] = 0;
448 for (short k=0; k<3; k++) {
449 Value[c] += (int32_t)(MatrixRGBToXYZ[m_ColorSpace][c][k]*m_Image[i][k]);
450 }
451 }
452 for (short c=0; c<3; c++) {
453 m_Image[i][c] = (uint16_t) CLIP(Value[c]>>1); // XYZ 1.0 maps on 8000H
454 }
455 }
456
457 m_ColorSpace = ptSpace_XYZ;
458
459 return this;
460 }
461
462 ////////////////////////////////////////////////////////////////////////////////
463 //
464 // lcmsRGBToXYZ.
465 // Conversion from RGB space to XYZ space using lcms.
466 //
467 ////////////////////////////////////////////////////////////////////////////////
468
lcmsRGBToXYZ(const int Intent)469 ptImage* ptImage::lcmsRGBToXYZ(const int Intent) {
470
471 assert ((m_ColorSpace>0) && (m_ColorSpace<5));
472
473 cmsToneCurve* Gamma = cmsBuildGamma(NULL, 1.0);
474 cmsToneCurve* Gamma3[3];
475 Gamma3[0] = Gamma3[1] = Gamma3[2] = Gamma;
476
477 cmsHPROFILE InProfile = 0;
478
479 cmsCIExyY DFromReference;
480
481 switch (m_ColorSpace) {
482 case ptSpace_sRGB_D65 :
483 case ptSpace_AdobeRGB_D65 :
484 DFromReference = D65;
485 break;
486 case ptSpace_WideGamutRGB_D50 :
487 case ptSpace_ProPhotoRGB_D50 :
488 DFromReference = D50;
489 break;
490 default:
491 assert(0);
492 }
493
494 InProfile = cmsCreateRGBProfile(&DFromReference,
495 (cmsCIExyYTRIPLE*)&RGBPrimaries[m_ColorSpace],
496 Gamma3);
497
498 if (!InProfile) {
499 ptLogError(ptError_Profile,"Could not open InProfile profile.");
500 return NULL;
501 }
502
503 cmsHPROFILE OutProfile = 0;
504
505 OutProfile = cmsCreateXYZProfile();
506
507 if (!OutProfile) {
508 ptLogError(ptError_Profile,"Could not open OutProfile profile.");
509 cmsCloseProfile (InProfile);
510 return NULL;
511 }
512
513 cmsFreeToneCurve(Gamma);
514
515 cmsHTRANSFORM Transform;
516 Transform = cmsCreateTransform(InProfile,
517 TYPE_RGB_16,
518 OutProfile,
519 TYPE_XYZ_16,
520 Intent,
521 cmsFLAGS_NOOPTIMIZE | cmsFLAGS_BLACKPOINTCOMPENSATION);
522
523 int32_t Size = m_Width*m_Height;
524 int32_t Step = 100000;
525 #pragma omp parallel for schedule(static)
526 for (int32_t i = 0; i < Size; i+=Step) {
527 int32_t Length = (i+Step)<Size ? Step : Size - i;
528 uint16_t* Tile = &(m_Image[i][0]);
529 cmsDoTransform(Transform,Tile,Tile,Length);
530 }
531 cmsDeleteTransform(Transform);
532 cmsCloseProfile(OutProfile);
533 cmsCloseProfile(InProfile);
534
535 m_ColorSpace = ptSpace_XYZ;
536
537 return this;
538 }
539
540 ////////////////////////////////////////////////////////////////////////////////
541 //
542 // XYZToRGB.
543 // Conversion from XYZ space to RGB space.
544 //
545 ////////////////////////////////////////////////////////////////////////////////
546
XYZToRGB(const short To)547 ptImage* ptImage::XYZToRGB(const short To) {
548
549 assert ((To>0) && (To<5));
550 assert (m_ColorSpace = ptSpace_XYZ);
551
552 // Convert the image.
553 #pragma omp parallel for
554 for (uint32_t i=0; i<(uint32_t)m_Height*m_Width; i++) {
555 int32_t Value[3];
556 for (short c=0; c<3; c++) {
557 Value[c] = 0;
558 for (short k=0; k<3; k++) {
559 Value[c] += (int32_t) (MatrixXYZToRGB[To][c][k]*m_Image[i][k]);
560 }
561 }
562 for (short c=0; c<3; c++) {
563 m_Image[i][c] = (uint16_t) CLIP(Value[c]<<1); // XYZ 1.0 maps on 8000H
564 }
565 }
566
567 m_ColorSpace = To;
568
569 return this;
570 }
571
572
573 ////////////////////////////////////////////////////////////////////////////////
574 //
575 // lcmsXYZToRGB.
576 // Conversion from XYZ space to RGB space using lcms.
577 //
578 ////////////////////////////////////////////////////////////////////////////////
579
lcmsXYZToRGB(const short To,const int Intent)580 ptImage* ptImage::lcmsXYZToRGB(const short To,
581 const int Intent) {
582
583 assert ((To>0) && (To<5));
584 assert (m_ColorSpace = ptSpace_XYZ);
585
586 cmsToneCurve* Gamma = cmsBuildGamma(NULL, 1.0);
587 cmsToneCurve* Gamma3[3];
588 Gamma3[0] = Gamma3[1] = Gamma3[2] = Gamma;
589
590 cmsHPROFILE InProfile = 0;
591
592 InProfile = cmsCreateXYZProfile();
593
594 if (!InProfile) {
595 ptLogError(ptError_Profile,"Could not open InProfile profile.");
596 return NULL;
597 }
598
599 cmsHPROFILE OutProfile = 0;
600
601 cmsCIExyY DToReference;
602
603 switch (To) {
604 case ptSpace_sRGB_D65 :
605 case ptSpace_AdobeRGB_D65 :
606 DToReference = D65;
607 break;
608 case ptSpace_WideGamutRGB_D50 :
609 case ptSpace_ProPhotoRGB_D50 :
610 DToReference = D50;
611 break;
612 default:
613 assert(0);
614 }
615
616 OutProfile = cmsCreateRGBProfile(&DToReference,
617 (cmsCIExyYTRIPLE*)&RGBPrimaries[To],
618 Gamma3);
619 if (!OutProfile) {
620 ptLogError(ptError_Profile,"Could not open OutProfile profile.");
621 cmsCloseProfile (InProfile);
622 return NULL;
623 }
624
625 cmsFreeToneCurve(Gamma);
626
627 cmsHTRANSFORM Transform;
628 Transform = cmsCreateTransform(InProfile,
629 TYPE_XYZ_16,
630 OutProfile,
631 TYPE_RGB_16,
632 Intent,
633 cmsFLAGS_NOOPTIMIZE | cmsFLAGS_BLACKPOINTCOMPENSATION);
634
635 int32_t Size = m_Width*m_Height;
636 int32_t Step = 100000;
637 #pragma omp parallel for schedule(static)
638 for (int32_t i = 0; i < Size; i+=Step) {
639 int32_t Length = (i+Step)<Size ? Step : Size - i;
640 uint16_t* Tile = &(m_Image[i][0]);
641 cmsDoTransform(Transform,Tile,Tile,Length);
642 }
643 cmsDeleteTransform(Transform);
644 cmsCloseProfile(OutProfile);
645 cmsCloseProfile(InProfile);
646
647 m_ColorSpace = To;
648
649 return this;
650 }
651
652 ////////////////////////////////////////////////////////////////////////////////
653 //
654 // RGBToLab
655 // Conversion from RGB space to Lab space.
656 //
657 ////////////////////////////////////////////////////////////////////////////////
658
659 // A lookup table used for this purpose.
660 // Initialized at the first call.
661 static double ToLABFunctionTable[0x20000];
662 static short ToLABFunctionInited = 0;
663
RGBToLab()664 ptImage* ptImage::RGBToLab() {
665
666 if (m_ColorSpace == ptSpace_Lab) return this;
667
668 assert (3 == m_Colors);
669
670 double DReference[3];
671 switch (m_ColorSpace) {
672 case ptSpace_sRGB_D65:
673 case ptSpace_AdobeRGB_D65:
674 DReference[0] = D65Reference[0];
675 DReference[1] = D65Reference[1];
676 DReference[2] = D65Reference[2];
677 break;
678 case ptSpace_WideGamutRGB_D50:
679 case ptSpace_ProPhotoRGB_D50:
680 DReference[0] = D50Reference[0];
681 DReference[1] = D50Reference[1];
682 DReference[2] = D50Reference[2];
683 break;
684 default:
685 assert(0);
686 }
687
688 // First go to XYZ
689 double xyz[3];
690 #pragma omp parallel for private(xyz)
691 for (uint32_t i=0; i < (uint32_t)m_Width*m_Height; i++) {
692 xyz[0] = MatrixRGBToXYZ[m_ColorSpace][0][0] * m_Image[i][0] +
693 MatrixRGBToXYZ[m_ColorSpace][0][1] * m_Image[i][1] +
694 MatrixRGBToXYZ[m_ColorSpace][0][2] * m_Image[i][2] ;
695 xyz[1] = MatrixRGBToXYZ[m_ColorSpace][1][0] * m_Image[i][0] +
696 MatrixRGBToXYZ[m_ColorSpace][1][1] * m_Image[i][1] +
697 MatrixRGBToXYZ[m_ColorSpace][1][2] * m_Image[i][2] ;
698 xyz[2] = MatrixRGBToXYZ[m_ColorSpace][2][0] * m_Image[i][0] +
699 MatrixRGBToXYZ[m_ColorSpace][2][1] * m_Image[i][1] +
700 MatrixRGBToXYZ[m_ColorSpace][2][2] * m_Image[i][2] ;
701
702 // Reference White
703 xyz[0] /= DReference[0];
704 xyz[1] /= DReference[1];
705 xyz[2] /= DReference[2];
706
707 // Then to Lab
708 xyz[0] = ToLABFunctionTable[ (int32_t) MAX(0.0,xyz[0]+0.5) ];
709 xyz[1] = ToLABFunctionTable[ (int32_t) MAX(0.0,xyz[1]+0.5) ];
710 xyz[2] = ToLABFunctionTable[ (int32_t) MAX(0.0,xyz[2]+0.5) ];
711
712 // L in 0 , a in 1, b in 2
713 m_Image[i][0] = CLIP((int32_t) (0xffff*(116.0/100.0 * xyz[1] - 16.0/100.0)));
714 m_Image[i][1] = CLIP((int32_t) (0x101*(128.0+500.0*(xyz[0]-xyz[1]))));
715 m_Image[i][2] = CLIP((int32_t) (0x101*(128.0+200.0*(xyz[1]-xyz[2]))));
716
717 }
718
719 // And that's it.
720 m_ColorSpace = ptSpace_Lab;
721
722 return this;
723 }
724
725 ////////////////////////////////////////////////////////////////////////////////
726 //
727 // lcmsRGBToLab
728 // Conversion from RGB space to Lab space using lcms.
729 //
730 ////////////////////////////////////////////////////////////////////////////////
731
lcmsRGBToLab(const int Intent)732 ptImage* ptImage::lcmsRGBToLab(const int Intent) {
733
734 assert (3 == m_Colors);
735 assert ((m_ColorSpace>0) && (m_ColorSpace<5));
736
737 cmsToneCurve* Gamma = cmsBuildGamma(NULL, 1.0);
738 cmsToneCurve* Gamma3[3];
739 Gamma3[0] = Gamma3[1] = Gamma3[2] = Gamma;
740
741 cmsHPROFILE InProfile = 0;
742
743 cmsCIExyY DFromReference;
744
745 switch (m_ColorSpace) {
746 case ptSpace_sRGB_D65 :
747 case ptSpace_AdobeRGB_D65 :
748 DFromReference = D65;
749 break;
750 case ptSpace_WideGamutRGB_D50 :
751 case ptSpace_ProPhotoRGB_D50 :
752 DFromReference = D50;
753 break;
754 default:
755 assert(0);
756 }
757
758 InProfile = cmsCreateRGBProfile(&DFromReference,
759 (cmsCIExyYTRIPLE*)&RGBPrimaries[m_ColorSpace],
760 Gamma3);
761
762 if (!InProfile) {
763 ptLogError(ptError_Profile,"Could not open InProfile profile.");
764 return NULL;
765 }
766
767 cmsHPROFILE OutProfile = 0;
768
769 OutProfile = cmsCreateLab4Profile(NULL);
770
771 if (!OutProfile) {
772 ptLogError(ptError_Profile,"Could not open OutProfile profile.");
773 cmsCloseProfile (InProfile);
774 return NULL;
775 }
776
777 cmsFreeToneCurve(Gamma);
778
779 cmsHTRANSFORM Transform;
780 Transform = cmsCreateTransform(InProfile,
781 TYPE_RGB_16,
782 OutProfile,
783 TYPE_Lab_16,
784 Intent,
785 cmsFLAGS_NOOPTIMIZE | cmsFLAGS_BLACKPOINTCOMPENSATION);
786
787 int32_t Size = m_Width*m_Height;
788 int32_t Step = 100000;
789 #pragma omp parallel for schedule(static)
790 for (int32_t i = 0; i < Size; i+=Step) {
791 int32_t Length = (i+Step)<Size ? Step : Size - i;
792 uint16_t* Tile = &(m_Image[i][0]);
793 cmsDoTransform(Transform,Tile,Tile,Length);
794 }
795 cmsDeleteTransform(Transform);
796 cmsCloseProfile(OutProfile);
797 cmsCloseProfile(InProfile);
798
799 // And that's it.
800 m_ColorSpace = ptSpace_Lab;
801
802 return this;
803 }
804
805 ////////////////////////////////////////////////////////////////////////////////
806 //
807 // LabToRGB
808 // Conversion from Lab space to RGB space.
809 //
810 ////////////////////////////////////////////////////////////////////////////////
811
LabToRGB(const short To)812 ptImage* ptImage::LabToRGB(const short To) {
813
814 if (m_ColorSpace == To) return this;
815
816 if (!(m_ColorSpace == ptSpace_Lab)) {
817 ptMessageBox::critical(0,"Error","Too fast! Keep cool ;-)");
818 return this;
819 }
820
821 assert (3 == m_Colors);
822 assert ((To>0) && (To<5));
823
824 double DReference[3];
825 switch (To) {
826 case ptSpace_sRGB_D65:
827 case ptSpace_AdobeRGB_D65:
828 DReference[0] = D65Reference[0];
829 DReference[1] = D65Reference[1];
830 DReference[2] = D65Reference[2];
831 break;
832 case ptSpace_WideGamutRGB_D50:
833 case ptSpace_ProPhotoRGB_D50:
834 DReference[0] = D50Reference[0];
835 DReference[1] = D50Reference[1];
836 DReference[2] = D50Reference[2];
837 break;
838 default:
839 assert(0);
840 }
841
842 // This code has been tweaked with performance in mind.
843 // Probably some of it could be detected by compiler, but overall I
844 // still gain something like going from 167 to 35 clockticks.
845 // (especially avoiding the stupid pow(x,3.0) -> x*x*x !
846
847 double xyz[3];
848 double fx,fy,fz;
849 double xr,yr,zr;
850 const double epsilon = 216.0/24389.0;
851 const double kappa = 24389.0/27.0;
852
853 const uint32_t NrPixels = m_Width*m_Height;
854 #pragma omp parallel for private(xyz,fx,fy,fz,xr,yr,zr)
855 for (uint32_t i=0; i < NrPixels ; i++) {
856
857 const double L = ToFloatTable[m_Image[i][0]]*100.0;
858 const double a = (m_Image[i][1]-0x8080) /((double)0x8080/128.0) ;
859 const double b = (m_Image[i][2]-0x8080) /((double)0x8080/128.0) ;
860
861 const double Tmp1 = (L+16.0)/116.0;
862
863 yr = (L<=kappa*epsilon)?
864 (L/kappa):(Tmp1*Tmp1*Tmp1);
865 fy = (yr<=epsilon) ? ((kappa*yr+16.0)/116.0) : Tmp1;
866 fz = fy - b/200.0;
867 fx = a/500.0 + fy;
868 const double fz3 = fz*fz*fz;
869 const double fx3 = fx*fx*fx;
870 zr = (fz3<=epsilon) ? ((116.0*fz-16.0)/kappa) : fz3;
871 xr = (fx3<=epsilon) ? ((116.0*fx-16.0)/kappa) : fx3;
872
873 xyz[0] = xr*DReference[0]*65535.0 - 0.5;
874 xyz[1] = yr*DReference[1]*65535.0 - 0.5;
875 xyz[2] = zr*DReference[2]*65535.0 - 0.5;
876
877 // And finally to RGB via the matrix.
878 for (short c=0; c<3; c++) {
879 double Value = 0;
880 Value += MatrixXYZToRGB[To][c][0] * xyz[0];
881 Value += MatrixXYZToRGB[To][c][1] * xyz[1];
882 Value += MatrixXYZToRGB[To][c][2] * xyz[2];
883 m_Image[i][c] = (uint16_t) CLIP((int32_t)(Value));
884 }
885 }
886
887 // And that's it.
888
889 m_ColorSpace = To;
890
891 return this;
892 }
893
894 ////////////////////////////////////////////////////////////////////////////////
895 //
896 // lcmsLabToRGB
897 // Conversion from Lab space to RGB space using lcms.
898 //
899 ////////////////////////////////////////////////////////////////////////////////
900
lcmsLabToRGB(const short To,const int Intent)901 ptImage* ptImage::lcmsLabToRGB(const short To,
902 const int Intent) {
903
904 assert (3 == m_Colors);
905 assert ((To>0) && (To<5));
906 assert (m_ColorSpace == ptSpace_Lab);
907
908 cmsToneCurve* Gamma = cmsBuildGamma(NULL, 1.0);
909 cmsToneCurve* Gamma3[3];
910 Gamma3[0] = Gamma3[1] = Gamma3[2] = Gamma;
911
912 cmsHPROFILE InProfile = 0;
913
914 InProfile = cmsCreateLab4Profile(NULL);
915
916 if (!InProfile) {
917 ptLogError(ptError_Profile,"Could not open InProfile profile.");
918 return NULL;
919 }
920
921 cmsHPROFILE OutProfile = 0;
922
923 cmsCIExyY DToReference;
924
925 switch (To) {
926 case ptSpace_sRGB_D65 :
927 case ptSpace_AdobeRGB_D65 :
928 DToReference = D65;
929 break;
930 case ptSpace_WideGamutRGB_D50 :
931 case ptSpace_ProPhotoRGB_D50 :
932 DToReference = D50;
933 break;
934 default:
935 assert(0);
936 }
937
938 OutProfile = cmsCreateRGBProfile(&DToReference,
939 (cmsCIExyYTRIPLE*)&RGBPrimaries[To],
940 Gamma3);
941
942 if (!OutProfile) {
943 ptLogError(ptError_Profile,"Could not open OutProfile profile.");
944 cmsCloseProfile (InProfile);
945 return NULL;
946 }
947
948 cmsFreeToneCurve(Gamma);
949
950 cmsHTRANSFORM Transform;
951 Transform = cmsCreateTransform(InProfile,
952 TYPE_Lab_16,
953 OutProfile,
954 TYPE_RGB_16,
955 Intent,
956 cmsFLAGS_NOOPTIMIZE | cmsFLAGS_BLACKPOINTCOMPENSATION);
957
958 int32_t Size = m_Width*m_Height;
959 int32_t Step = 100000;
960 #pragma omp parallel for schedule(static)
961 for (int32_t i = 0; i < Size; i+=Step) {
962 int32_t Length = (i+Step)<Size ? Step : Size - i;
963 uint16_t* Tile = &(m_Image[i][0]);
964 cmsDoTransform(Transform,Tile,Tile,Length);
965 }
966 cmsDeleteTransform(Transform);
967 cmsCloseProfile(OutProfile);
968 cmsCloseProfile(InProfile);
969
970 // And that's it.
971
972 m_ColorSpace = To;
973
974 return this;
975 }
976
977 //==============================================================================
978
RGBToLch()979 ptImage *ptImage::RGBToLch() {
980 if (m_ColorSpace == ptSpace_LCH) return this;
981
982 this->RGBToLab();
983 this->LabToLch();
984
985 return this;
986 }
987
988 //==============================================================================
989
LchToRGB(const short To)990 ptImage *ptImage::LchToRGB(const short To) {
991 if (m_ColorSpace == To) return this;
992
993 LchToLab();
994 LabToRGB(To);
995 return this;
996 }
997
998 //==============================================================================
999
ResizeLCH(size_t ASize)1000 void ptImage::ResizeLCH(size_t ASize) {
1001 m_ImageL.resize(ASize);
1002 m_ImageL.shrink_to_fit();
1003 m_ImageC.resize(ASize);
1004 m_ImageC.shrink_to_fit();
1005 m_ImageH.resize(ASize);
1006 m_ImageH.shrink_to_fit();
1007 }
1008
1009 //==============================================================================
1010
LabToLch()1011 ptImage *ptImage::LabToLch()
1012 {
1013 if (m_ColorSpace == ptSpace_LCH) return this;
1014
1015 assert (m_ColorSpace == ptSpace_Lab);
1016
1017 uint32_t hSize = (uint32_t)m_Width*m_Height;
1018
1019 ResizeLCH(hSize);
1020
1021 float hValueA = 0.0f;
1022 float hValueB = 0.0f;
1023 #pragma omp parallel for schedule(static) private(hValueA, hValueB)
1024 for (uint32_t i = 0; i < hSize; i++) {
1025 hValueA = ToFloatABNeutral[m_Image[i][1]];
1026 hValueB = ToFloatABNeutral[m_Image[i][2]];
1027 m_ImageL.at(i) = m_Image[i][0];
1028 m_ImageH.at(i) = ToHue(hValueA, hValueB);
1029 m_ImageC.at(i) = powf(hValueA*hValueA + hValueB*hValueB, 0.5f);
1030 }
1031
1032 setSize(0);
1033 m_ColorSpace = ptSpace_LCH;
1034 return this;
1035 }
1036
1037 //==============================================================================
1038
toRGB()1039 ptImage *ptImage::toRGB()
1040 {
1041 if (m_ColorSpace == ptSpace_Lab) LabToRGB(getCurrentRGB());
1042 else if (m_ColorSpace == ptSpace_XYZ) XYZToRGB(getCurrentRGB());
1043 else if (m_ColorSpace == ptSpace_Profiled) GInfo->Raise("Cannot transform color space!");
1044
1045 return this;
1046 }
1047
1048 //==============================================================================
1049
toLab()1050 ptImage *ptImage::toLab()
1051 {
1052 if (m_ColorSpace == ptSpace_Lab) return this;
1053 else if (m_ColorSpace == ptSpace_XYZ) XYZToRGB(getCurrentRGB());
1054 else if (m_ColorSpace == ptSpace_Profiled) GInfo->Raise("Cannot transform color space!");
1055
1056 RGBToLab();
1057
1058 return this;
1059 }
1060
1061 //==============================================================================
1062
LchToLab()1063 ptImage *ptImage::LchToLab() {
1064 if (m_ColorSpace == ptSpace_Lab) return this;
1065
1066 assert (m_ColorSpace == ptSpace_LCH);
1067 assert (m_Image == 0);
1068
1069 uint32_t hSize = (uint32_t)m_Width*m_Height;
1070 setSize((size_t)hSize);
1071
1072 #pragma omp parallel for schedule(static)
1073 for (uint32_t i = 0; i < hSize; i++) {
1074 m_Image[i][0] = m_ImageL.at(i);
1075 m_Image[i][1] = CLIP((int32_t)(cosf(m_ImageH.at(i))*m_ImageC.at(i) + ptWPHLab));
1076 m_Image[i][2] = CLIP((int32_t)(sinf(m_ImageH.at(i))*m_ImageC.at(i) + ptWPHLab));
1077 }
1078
1079 ResizeLCH(0);
1080
1081 m_ColorSpace = ptSpace_Lab;
1082 return this;
1083 }
1084
1085 ////////////////////////////////////////////////////////////////////////////////
1086 //
1087 // Constructor.
1088 //
1089 ////////////////////////////////////////////////////////////////////////////////
1090
ptImage()1091 ptImage::ptImage() {
1092 m_Width = 0;
1093 m_Height = 0;
1094 m_Image = nullptr;
1095 m_Colors = 0;
1096 m_ColorSpace = ptSpace_sRGB_D65;
1097 m_Data.clear();
1098 ResizeLCH(0);
1099
1100 // Initialize the lookup table for the RGB->LAB function
1101 // if this would be the first time.
1102 if (!ToLABFunctionInited) {
1103 // Remark : we extend the table well beyond r>1.0 for numerical
1104 // stability purposes. XYZ>1.0 occurs sometimes and this way
1105 // it stays stable (srgb->lab->srgb correct within 0.02%)
1106 #pragma omp parallel for
1107 for (uint32_t i=0; i<0x20000; i++) {
1108 double r = (double)(i) / ptWP;
1109 ToLABFunctionTable[i] = r > 216.0/24389.0 ? pow(r,1/3.0) : (24389.0/27.0*r + 16.0)/116.0;
1110 }
1111 ToLABFunctionInited = 1;
1112 }
1113
1114 // Some lcsm initialization.
1115 //cmsErrorAction (LCMS_ERROR_SHOW);
1116 cmsWhitePointFromTemp(&D65, 6503);
1117 cmsWhitePointFromTemp(&D50, 5003);
1118 }
1119
1120 ////////////////////////////////////////////////////////////////////////////////
1121 //
1122 // Destructor.
1123 //
1124 ////////////////////////////////////////////////////////////////////////////////
1125
~ptImage()1126 ptImage::~ptImage() {
1127 // nothing to free :-)
1128 }
1129
1130 //==============================================================================
1131
setSize(size_t Size)1132 void ptImage::setSize(size_t Size) {
1133 m_Data.resize(Size);
1134 m_Data.shrink_to_fit();
1135 m_Image = (uint16_t (*)[3]) m_Data.data();
1136 }
1137
1138 // -----------------------------------------------------------------------------
1139
clear()1140 void ptImage::clear() {
1141 this->setSize(0);
1142 m_Width = 0;
1143 m_Height = 0;
1144 }
1145
1146 ////////////////////////////////////////////////////////////////////////////////
1147 //
1148 // Set
1149 //
1150 ////////////////////////////////////////////////////////////////////////////////
1151
Set(const ptDcRaw * DcRawObject,const short TargetSpace,const char * ProfileName,const int Intent,const int ProfileGamma)1152 ptImage* ptImage::Set(const ptDcRaw* DcRawObject,
1153 const short TargetSpace,
1154 const char* ProfileName,
1155 const int Intent,
1156 const int ProfileGamma) {
1157
1158 assert(NULL != DcRawObject);
1159 assert ((TargetSpace>0) && (TargetSpace<5));
1160
1161 m_Width = DcRawObject->m_Width;
1162 m_Height = DcRawObject->m_Height;
1163
1164 // Temp image for later flip
1165 TImage16Data PreFlip;
1166 PreFlip.resize((size_t)m_Width*m_Height);
1167
1168 if (!ProfileName) {
1169 // Matrix for conversion RGB to RGB : is a multiplication via XYZ
1170 double MatrixRGBToRGB[3][3];
1171 for(short i=0;i<3;i++) {
1172 for (short j=0;j<3;j++) {
1173 MatrixRGBToRGB[i][j] = 0.0;
1174 for (short k=0;k<3;k++) {
1175 MatrixRGBToRGB[i][j] +=
1176 MatrixXYZToRGB[TargetSpace][i][k] *
1177 // Yes dcraw assumes that the result of rgb_cam is in sRGB !
1178 MatrixRGBToXYZ[ptSpace_sRGB_D65][k][j];
1179 }
1180 }
1181 }
1182
1183 // We just need to add the correct matrix calculation with rgb_cam
1184 // Be attentive : we still can have 4 channels at this moment !
1185 double Matrix[3][4];
1186 for(short i=0;i<3;i++) {
1187 for (short j=0;j<DcRawObject->m_Colors;j++) {
1188 Matrix[i][j] = 0.0;
1189 for (short k=0;k<3;k++) {
1190 Matrix[i][j] +=
1191 MatrixRGBToRGB[i][k] * DcRawObject->m_MatrixCamRGBToSRGB[k][j];
1192 }
1193 }
1194 }
1195
1196 // Convert the image.
1197 #pragma omp parallel for
1198 for (uint32_t i=0; i<(uint32_t)m_Height*m_Width; i++) {
1199 int32_t Value[3];
1200 for (short c=0; c<3; c++) {
1201 Value[c] = 0;
1202 for (short k=0; k<DcRawObject->m_Colors; k++) {
1203 Value[c] += (int32_t) (Matrix[c][k]*DcRawObject->m_Image[i][k]);
1204 }
1205 }
1206 for (short c=0; c<3; c++) {
1207 PreFlip[i][c] = (uint16_t) CLIP(Value[c]);
1208 }
1209 }
1210
1211 m_Colors = MIN((int)(DcRawObject->m_Colors),3);
1212 m_ColorSpace = TargetSpace;
1213
1214 if (m_Colors != 3) {
1215 fprintf(stderr,
1216 "At this moment anything with m_Colors != 3 was not yet "
1217 "tested sufficiently. You can remove the assertion related "
1218 "to this, but unexpected issues can pop up. Consider submitting "
1219 "a picture from this camera to the author in order to test "
1220 "further.\n");
1221 assert(0);
1222 }
1223 }// (if !ProfileName)
1224 else if (ProfileName) {
1225
1226 m_Colors = MIN((int)(DcRawObject->m_Colors),3);
1227 m_ColorSpace = TargetSpace;
1228
1229 if (m_Colors != 3) {
1230 fprintf(stderr,
1231 "At this moment anything with m_Colors != 3 was not yet "
1232 "tested sufficiently. You can remove the assertion related "
1233 "to this, but unexpected issues can pop up. Consider submitting "
1234 "a picture from this camera to the author in order to test "
1235 "further. Moreover , is this profile correct ??\n");
1236 assert(0);
1237 }
1238
1239 short NrProfiles = (ProfileGamma) ? 3:2;
1240 cmsHPROFILE Profiles[NrProfiles];
1241 short ProfileIdx = 0;
1242 if (ProfileGamma) {
1243 cmsHPROFILE PreProfile = 0;
1244 double Params[6];
1245 switch(ProfileGamma) {
1246 case ptCameraColorGamma_sRGB :
1247 // Parameters for standard (inverse) sRGB in lcms
1248 Params[0] = 1.0/2.4;
1249 Params[1] = 1.1371189;
1250 Params[2] = 0.0;
1251 Params[3] = 12.92;
1252 Params[4] = 0.0031308;
1253 Params[5] = -0.055;
1254 break;
1255 case ptCameraColorGamma_BT709 :
1256 // Parameters for standard (inverse) BT709 in lcms
1257 Params[0] = 0.45;
1258 Params[1] = 1.233405791;
1259 Params[2] = 0.0;
1260 Params[3] = 4.5;
1261 Params[4] = 0.018;
1262 Params[5] = -0.099;
1263 break;
1264 case ptCameraColorGamma_Pure22 :
1265 // Parameters for standard (inverse) 2.2 gamma in lcms
1266 Params[0] = 1.0/2.2;
1267 Params[1] = 1.0;
1268 Params[2] = 0.0;
1269 Params[3] = 0.0;
1270 Params[4] = 0.0;
1271 Params[5] = 0.0;
1272 break;
1273 default:
1274 assert(0);
1275 }
1276 cmsToneCurve* Gamma = cmsBuildParametricToneCurve(0,4,Params);
1277 cmsToneCurve* Gamma4[4];
1278 Gamma4[0] = Gamma4[1] = Gamma4[2] = Gamma4[3]= Gamma;
1279 PreProfile = cmsCreateLinearizationDeviceLink(cmsSigRgbData,Gamma4);
1280 Profiles[ProfileIdx] = PreProfile;
1281 if (!Profiles[ProfileIdx]) {
1282 ptLogError(ptError_Profile,"Could not open sRGB preprofile.");
1283 return NULL;
1284 }
1285 ProfileIdx++;
1286 cmsFreeToneCurve(Gamma);
1287 }
1288
1289 Profiles[ProfileIdx] = cmsOpenProfileFromFile(ProfileName,"r");
1290 if (!Profiles[ProfileIdx]) {
1291 ptLogError(ptError_Profile,"Could not open profile %s.",ProfileName);
1292 for (; ProfileIdx>=0; ProfileIdx--) {
1293 cmsCloseProfile(Profiles[ProfileIdx]);
1294 }
1295 return NULL;
1296 }
1297 ProfileIdx++;
1298
1299 // Calculate the output profile (with gamma 1)
1300
1301 // Linear gamma.
1302 cmsToneCurve* Gamma = cmsBuildGamma(NULL, 1.0);
1303 cmsToneCurve* Gamma3[3];
1304 Gamma3[0] = Gamma3[1] = Gamma3[2] = Gamma;
1305
1306 cmsCIExyY DToReference;
1307
1308 switch (TargetSpace) {
1309 case ptSpace_sRGB_D65 :
1310 case ptSpace_AdobeRGB_D65 :
1311 DToReference = D65;
1312 break;
1313 case ptSpace_WideGamutRGB_D50 :
1314 case ptSpace_ProPhotoRGB_D50 :
1315 DToReference = D50;
1316 break;
1317 default:
1318 assert(0);
1319 }
1320
1321 Profiles[ProfileIdx] =
1322 cmsCreateRGBProfile(&DToReference,
1323 (cmsCIExyYTRIPLE*)
1324 &RGBPrimaries[TargetSpace],
1325 Gamma3);
1326 if (!Profiles[ProfileIdx]) {
1327 ptLogError(ptError_Profile,"Could not open OutProfile profile.");
1328 for (; ProfileIdx>=0; ProfileIdx--) {
1329 cmsCloseProfile(Profiles[ProfileIdx]);
1330 }
1331 return NULL;
1332 }
1333
1334 cmsFreeToneCurve(Gamma);
1335
1336 cmsHTRANSFORM Transform;
1337 Transform = cmsCreateMultiprofileTransform(Profiles,
1338 NrProfiles,
1339 TYPE_RGBA_16,
1340 TYPE_RGB_16,
1341 Intent,
1342 0);
1343 int32_t Size = m_Width*m_Height;
1344 int32_t Step = 100000;
1345 #pragma omp parallel for schedule(static)
1346 for (int32_t i = 0; i < Size; i+=Step) {
1347 int32_t Length = (i+Step)<Size ? Step : Size - i;
1348 uint16_t* Tile1 = &(DcRawObject->m_Image[i][0]);
1349 uint16_t* Tile2 = &(PreFlip[i][0]);
1350 cmsDoTransform(Transform,Tile1,Tile2,Length);
1351 }
1352 cmsDeleteTransform(Transform);
1353 for (; ProfileIdx>=0; ProfileIdx--) {
1354 cmsCloseProfile(Profiles[ProfileIdx]);
1355 }
1356 }
1357
1358 uint16_t TargetWidth = m_Width;
1359 uint16_t TargetHeight = m_Height;
1360 // Free and allocate
1361 setSize((int32_t)TargetWidth*TargetHeight);
1362
1363 // Flip image. With m_Flip as in dcraw.
1364 // (see also flip_index() function in dcraw)
1365 if (DcRawObject->m_Flip & 4) {
1366 SWAP(TargetWidth,TargetHeight);
1367 }
1368 #pragma omp parallel for
1369 for (uint16_t TargetRow=0; TargetRow<TargetHeight; TargetRow++) {
1370 for (uint16_t TargetCol=0; TargetCol<TargetWidth; TargetCol++) {
1371 uint16_t OriginRow = TargetRow;
1372 uint16_t OriginCol = TargetCol;
1373 if (DcRawObject->m_Flip & 4) SWAP(OriginRow,OriginCol);
1374 if (DcRawObject->m_Flip & 2) OriginRow = m_Height-1-OriginRow;
1375 if (DcRawObject->m_Flip & 1) OriginCol = m_Width-1-OriginCol;
1376 for (short c=0; c<3; c++) {
1377 m_Data[TargetRow*TargetWidth+TargetCol][c] =
1378 PreFlip[OriginRow*m_Width+OriginCol][c];
1379 }
1380 }
1381 }
1382
1383 m_Height = TargetHeight;
1384 m_Width = TargetWidth;
1385
1386 return this;
1387 }
1388
1389 ////////////////////////////////////////////////////////////////////////////////
1390 //
1391 // Set
1392 // To DcRaw's ImageAfterPhase2. (ad hoc, look at 'sensor' values)
1393 //
1394 ////////////////////////////////////////////////////////////////////////////////
1395
Set(const ptDcRaw * DcRawObject,const short TargetSpace)1396 ptImage* ptImage::Set(const ptDcRaw* DcRawObject,
1397 const short TargetSpace) {
1398
1399 assert(NULL != DcRawObject);
1400 assert ((TargetSpace>0) && (TargetSpace<5));
1401
1402 m_Width = DcRawObject->m_Width;
1403 m_Height = DcRawObject->m_Height;
1404 m_ColorSpace = TargetSpace;
1405
1406 // Temp image for later flip
1407 TImage16Data PreFlip;
1408 PreFlip.resize((size_t)m_Width*m_Height);
1409
1410 // Convert the image.
1411 #pragma omp parallel for
1412 for (uint32_t i=0; i<(uint32_t)m_Height*m_Width; i++) {
1413 for (short c=0; c<3; c++) {
1414 PreFlip[i][c] = DcRawObject->m_Image_AfterPhase2[i][c];
1415 }
1416 }
1417
1418 m_Colors = MIN((int)(DcRawObject->m_Colors),3);
1419
1420 if (m_Colors != 3) {
1421 fprintf(stderr,
1422 "At this moment anything with m_Colors != 3 was not yet "
1423 "tested sufficiently. You can remove the assertion related "
1424 "to this, but unexpected issues can pop up. Consider submitting "
1425 "a picture from this camera to the author in order to test "
1426 "further.\n");
1427 assert(0);
1428 }
1429
1430 uint16_t TargetWidth = m_Width;
1431 uint16_t TargetHeight = m_Height;
1432
1433 // Free a maybe preexisting and allocate space.
1434 setSize((size_t)TargetWidth*TargetHeight);
1435
1436 if (DcRawObject->m_Flip & 4) {
1437 SWAP(TargetWidth,TargetHeight);
1438 }
1439 #pragma omp parallel for schedule(static)
1440 for (uint16_t TargetRow=0; TargetRow<TargetHeight; TargetRow++) {
1441 for (uint16_t TargetCol=0; TargetCol<TargetWidth; TargetCol++) {
1442 uint16_t OriginRow = TargetRow;
1443 uint16_t OriginCol = TargetCol;
1444 if (DcRawObject->m_Flip & 4) SWAP(OriginRow,OriginCol);
1445 if (DcRawObject->m_Flip & 2) OriginRow = m_Height-1-OriginRow;
1446 if (DcRawObject->m_Flip & 1) OriginCol = m_Width-1-OriginCol;
1447 for (short c=0; c<3; c++) {
1448 m_Data[TargetRow*TargetWidth+TargetCol][c] =
1449 PreFlip[OriginRow*m_Width+OriginCol][c];
1450 }
1451 }
1452 }
1453
1454 m_Height = TargetHeight;
1455 m_Width = TargetWidth;
1456
1457 return this;
1458 }
1459
1460 ////////////////////////////////////////////////////////////////////////////////
1461 //
1462 // Set, just allocation
1463 //
1464 ////////////////////////////////////////////////////////////////////////////////
1465
Set(const uint16_t Width,const uint16_t Height)1466 ptImage* ptImage::Set(const uint16_t Width,
1467 const uint16_t Height) {
1468
1469 m_Width = Width;
1470 m_Height = Height;
1471
1472 // Free a maybe preexisting and allocate space.
1473 setSize((size_t)m_Width*m_Height);
1474
1475 m_Colors = 3;
1476 m_ColorSpace = ptSpace_sRGB_D65;
1477
1478 return this;
1479 }
1480
1481 ////////////////////////////////////////////////////////////////////////////////
1482 //
1483 // Set
1484 //
1485 ////////////////////////////////////////////////////////////////////////////////
1486
Set(const ptImage * Origin)1487 ptImage* ptImage::Set(const ptImage *Origin) { // Always deep
1488
1489 assert(NULL != Origin);
1490
1491 m_Width = Origin->m_Width;
1492 m_Height = Origin->m_Height;
1493 m_Colors = Origin->m_Colors;
1494 m_ColorSpace = Origin->m_ColorSpace;
1495 setSize((size_t)m_Width*m_Height);
1496 m_Data = Origin->m_Data;
1497
1498 return this;
1499 }
1500
1501 ////////////////////////////////////////////////////////////////////////////////
1502 //
1503 // Set scaled
1504 //
1505 ////////////////////////////////////////////////////////////////////////////////
1506
SetScaled(const ptImage * Origin,const short ScaleFactor)1507 ptImage* ptImage::SetScaled(const ptImage *Origin,
1508 const short ScaleFactor) {
1509
1510 assert(NULL != Origin);
1511
1512 m_Width = Origin->m_Width;
1513 m_Height = Origin->m_Height;
1514 m_Colors = Origin->m_Colors;
1515 m_ColorSpace = Origin->m_ColorSpace;
1516
1517 if (ScaleFactor == 0) {
1518 setSize((int32_t)m_Width*m_Height);
1519 m_Data = Origin->m_Data;
1520 } else {
1521 m_Width >>= ScaleFactor;
1522 m_Height >>= ScaleFactor;
1523
1524 short Step = 1 << ScaleFactor;
1525 float InvAverage = 1.0/powf(2.0,2.0 * ScaleFactor);
1526
1527 // Allocate new.
1528 setSize((int32_t)m_Width*m_Height);
1529
1530 #pragma omp parallel for schedule(static)
1531 for (uint16_t Row=0; Row < m_Height; Row++) {
1532 for (uint16_t Col=0; Col < m_Width; Col++) {
1533 float PixelValue[3] = {0.0,0.0,0.0};
1534 for (uint8_t sRow=0; sRow < Step; sRow++) {
1535 for (uint8_t sCol=0; sCol < Step; sCol++) {
1536 int32_t index = (Row*Step+sRow)*Origin->m_Width+Col*Step+sCol;
1537 for (short c=0; c < 3; c++) {
1538 PixelValue[c] += Origin->m_Data[index][c];
1539 }
1540 }
1541 }
1542 for (short c=0; c < 3; c++) {
1543 m_Data[Row*m_Width+Col][c]
1544 = (int32_t) (PixelValue[c] * InvAverage);
1545 }
1546 }
1547 }
1548 }
1549 return this;
1550 }
1551
1552 ////////////////////////////////////////////////////////////////////////////////
1553 //
1554 // Get the RGB value at a given point
1555 //
1556 ////////////////////////////////////////////////////////////////////////////////
1557
GetRGB(const uint16_t x,const uint16_t y)1558 RGBValue ptImage::GetRGB(const uint16_t x, const uint16_t y) {
1559 RGBValue RGB;
1560
1561 if (m_ColorSpace != ptSpace_Lab) {
1562 RGB.R = m_Image[y*m_Width+x][0];
1563 RGB.G = m_Image[y*m_Width+x][1];
1564 RGB.B = m_Image[y*m_Width+x][2];
1565 } else {
1566 RGB.R = RGB.G = RGB.B = 0;
1567 }
1568
1569 return RGB;
1570 }
1571
1572 ////////////////////////////////////////////////////////////////////////////////
1573 //
1574 // IndicateOverExposure
1575 //
1576 ////////////////////////////////////////////////////////////////////////////////
1577
IndicateExposure(const short Over,const short Under,const uint8_t ChannelMask,const uint16_t OverExposureLevel[3],const uint16_t UnderExposureLevel[3])1578 ptImage* ptImage::IndicateExposure(const short Over,
1579 const short Under,
1580 const uint8_t ChannelMask,
1581 const uint16_t OverExposureLevel[3],
1582 const uint16_t UnderExposureLevel[3]) {
1583
1584 if (!Over && !Under) return this;
1585 if (!ChannelMask) return this;
1586 #pragma omp parallel for
1587 for (uint32_t i=0; i< (uint32_t)m_Height*m_Width; i++) {
1588 for (short Color=0; Color<m_Colors; Color++) {
1589 if (! (ChannelMask & (1<<Color) )) continue;
1590 if (m_Image[i][Color] >= OverExposureLevel[Color]) {
1591 if (Over) m_Image[i][Color] = 0;
1592 } else if (m_Image[i][Color] <= UnderExposureLevel[Color]) {
1593 if (Under) m_Image[i][Color] = 0xffff;
1594 }
1595 }
1596 }
1597 return this;
1598 }
1599
IndicateExposure(ptImage * ValueImage,const short Over,const short Under,const uint8_t ChannelMask,const uint16_t OverExposureLevel[3],const uint16_t UnderExposureLevel[3])1600 ptImage* ptImage::IndicateExposure(ptImage* ValueImage,
1601 const short Over,
1602 const short Under,
1603 const uint8_t ChannelMask,
1604 const uint16_t OverExposureLevel[3],
1605 const uint16_t UnderExposureLevel[3]) {
1606
1607 if (!Over && !Under) return this;
1608 if (!ChannelMask) return this;
1609 #pragma omp parallel for
1610 for (uint32_t i=0; i< (uint32_t)m_Height*m_Width; i++) {
1611 for (short Color=0; Color<m_Colors; Color++) {
1612 if (! (ChannelMask & (1<<Color) )) continue;
1613 if (ValueImage->m_Image[i][Color] >= OverExposureLevel[Color]) {
1614 if (Over) m_Image[i][Color] = 0;
1615 } else if (ValueImage->m_Image[i][Color] <= UnderExposureLevel[Color]) {
1616 if (Under) m_Image[i][Color] = 0xffff;
1617 }
1618 }
1619 }
1620 return this;
1621 }
1622
1623 ////////////////////////////////////////////////////////////////////////////////
1624 //
1625 // Expose
1626 //
1627 ////////////////////////////////////////////////////////////////////////////////
1628
Expose(const double Exposure,const TExposureClipMode ExposureClipMode)1629 ptImage* ptImage::Expose(const double Exposure,
1630 const TExposureClipMode ExposureClipMode) {
1631
1632 assert (m_Colors == 3);
1633 assert (m_ColorSpace != ptSpace_XYZ);
1634 assert (Exposure < (2<<15));
1635 const short NrChannels = (m_ColorSpace == ptSpace_Lab)?1:3;
1636
1637 # pragma omp parallel for schedule(static)
1638 for (uint32_t i=0; i< (uint32_t)m_Height*m_Width; i++) {
1639 uint32_t Pixel[3];
1640 uint32_t Highest = 0;
1641 for (short Color=0; Color<NrChannels; Color++) {
1642 Pixel[Color] = (uint32_t)(m_Image[i][Color]*Exposure);
1643 if (Pixel[Color] > Highest) Highest = Pixel[Color];
1644 }
1645 if (Highest<= 0XFFFF) {
1646 // No clipping.
1647 for (short Color=0; Color<NrChannels; Color++) {
1648 m_Image[i][Color] = Pixel[Color];
1649 }
1650 } else {
1651 if (ExposureClipMode == TExposureClipMode::Hard) {
1652 for (short Color=0; Color<NrChannels; Color++) {
1653 m_Image[i][Color] = CLIP((int32_t)Pixel[Color]);
1654 }
1655 } else if (ExposureClipMode == TExposureClipMode::Ratio) {
1656 for (short Color=0; Color<NrChannels; Color++) {
1657 m_Image[i][Color] = CLIP((int32_t)(Pixel[Color]*(float)0xFFFF/Highest));
1658 }
1659 } else {
1660 assert("Unexpected clip mode.");
1661 }
1662 }
1663 }
1664
1665 return this;
1666 }
1667
1668 ////////////////////////////////////////////////////////////////////////////////
1669 //
1670 // FractionLevel
1671 //
1672 // Calculates the level in the image where Fraction of pixels is above.
1673 //
1674 ////////////////////////////////////////////////////////////////////////////////
1675
CalculateFractionLevel(const double Fraction,const uint8_t ChannelMask)1676 uint16_t ptImage::CalculateFractionLevel(const double Fraction,
1677 const uint8_t ChannelMask) {
1678 // Build histogram that is valid at this point.
1679 // (we might not have one yet or it might be on altered data)
1680 uint32_t Histogram[0x10000][4];
1681 memset (Histogram, 0, sizeof Histogram);
1682 #pragma omp parallel for schedule(static)
1683 for (uint32_t i=0; i<(uint32_t)m_Height*m_Width; i++) {
1684 for (short c=0; c<m_Colors; c++) {
1685 Histogram[m_Image[i][c]][c]++;
1686 }
1687 }
1688 uint32_t ExpectedPixels = (uint32_t) (m_Width*m_Height*Fraction);
1689 uint32_t Total;
1690 uint16_t Level = 0;
1691 for (short c=0;c<3;c++){
1692 if (! (ChannelMask & (1<<c) )) continue;
1693 Total = 0;
1694 uint16_t Value = 0;
1695 for (Value=0xFFFF; Value > 0 ; Value--) {
1696 if ((Total+= Histogram[Value][c]) > ExpectedPixels)
1697 break;
1698 }
1699 if (Level < Value) Level = Value;
1700 }
1701
1702 return Level;
1703 }
1704
1705
1706 ////////////////////////////////////////////////////////////////////////////////
1707 //
1708 // ApplyCurve
1709 //
1710 ////////////////////////////////////////////////////////////////////////////////
1711
1712
ApplyCurve(const ptCurve * Curve,const uint8_t ChannelMask)1713 ptImage* ptImage::ApplyCurve(const ptCurve *Curve,
1714 const uint8_t ChannelMask) {
1715
1716 assert (NULL != Curve);
1717 assert (m_Colors == 3);
1718 assert (m_ColorSpace != ptSpace_XYZ);
1719
1720 std::vector<short> Channel;
1721 if (ChannelMask & 1) Channel.push_back(0);
1722 if (ChannelMask & 2) Channel.push_back(1);
1723 if (ChannelMask & 4) Channel.push_back(2);
1724 std::for_each (m_Data.begin(), m_Data.end(), [&](TPixel16 &Pixel) {
1725 std::for_each (Channel.begin(), Channel.end(), [&](const short &Value){
1726 Pixel[Value] = Curve->Curve[ Pixel[Value] ];
1727 });
1728 });
1729
1730 return this;
1731 }
1732
1733 ////////////////////////////////////////////////////////////////////////////////
1734 //
1735 // Apply L by Hue Curve
1736 //
1737 ////////////////////////////////////////////////////////////////////////////////
1738
ApplyLByHueCurve(const ptCurve * Curve)1739 ptImage* ptImage::ApplyLByHueCurve(const ptCurve *Curve) {
1740
1741 // Best solution would be to use the Lab <-> Lch conversion from lcms.
1742 // This should be faster without sacrificing much quality.
1743
1744 assert (m_ColorSpace == ptSpace_Lab);
1745 // neutral value for a* and b* channel
1746 const float WPH = 0x8080;
1747
1748 std::for_each (m_Data.begin(), m_Data.end(), [&](TPixel16 &Pixel) {
1749 // Factor by hue
1750 float ValueA = (float)Pixel[1]-WPH;
1751 float ValueB = (float)Pixel[2]-WPH;
1752 float Hue = 0;
1753 if (ValueA == 0.0 && ValueB == 0.0) {
1754 Hue = 0; // value for grey pixel
1755 } else {
1756 Hue = atan2f(ValueB,ValueA);
1757 }
1758 while (Hue < 0) Hue += 2.*ptPI;
1759
1760 float Factor = Curve->Curve[CLIP((int32_t)(Hue/ptPI*WPH))]/(float)0x7fff - 1.0f;
1761 if (Factor != 0.0) {
1762 float Col = powf(ValueA * ValueA + ValueB * ValueB, 0.25) / (float) 0xb5;
1763 Factor = powf(2.0f, 3.0f*Factor*Col);
1764
1765 Pixel[0] = CLIP((int32_t)(Pixel[0] * Factor));
1766 }
1767 });
1768
1769 return this;
1770 }
1771
1772 ////////////////////////////////////////////////////////////////////////////////
1773 //
1774 // Apply Hue Curve
1775 //
1776 ////////////////////////////////////////////////////////////////////////////////
1777
ApplyHueCurve(const ptCurve * Curve)1778 ptImage* ptImage::ApplyHueCurve(const ptCurve *Curve) {
1779
1780 assert (m_ColorSpace == ptSpace_Lab);
1781 // neutral value for a* and b* channel
1782 const float WPH = 0x8080;
1783 const float ScalePi = ptPI / 0x7fff;
1784 const float InvScalePi = 0x7fff / ptPI;
1785
1786 float ValueA = 0.0;
1787 float ValueB = 0.0;
1788 float Col = 0.0;
1789 float Hue = 0.0;
1790
1791 if (Curve->mask() == ptCurve::ChromaMask) { // by chroma
1792 #pragma omp parallel for schedule(static) private(ValueA, ValueB, Col, Hue)
1793 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
1794
1795 ValueA = (float)m_Image[i][1]-WPH;
1796 ValueB = (float)m_Image[i][2]-WPH;
1797
1798 if (ValueA == 0.0 && ValueB == 0.0) {
1799 Hue = 0; // value for grey pixel
1800 } else {
1801 Hue = atan2f(ValueB,ValueA);
1802 }
1803 while (Hue < 0) Hue += 2.*ptPI;
1804 Col = powf(ValueA * ValueA + ValueB * ValueB, 0.5);
1805
1806 Hue += ((float)Curve->Curve[CLIP((int32_t)(Hue*InvScalePi))]-(float)0x7fff)*ScalePi;
1807
1808 m_Image[i][1] = CLIP((int32_t)((cosf(Hue)*Col)+WPH));
1809 m_Image[i][2] = CLIP((int32_t)((sinf(Hue)*Col)+WPH));
1810 }
1811
1812 } else { // by luma
1813 #pragma omp parallel for schedule(static) private(ValueA, ValueB, Col, Hue)
1814 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
1815
1816 ValueA = (float)m_Image[i][1]-WPH;
1817 ValueB = (float)m_Image[i][2]-WPH;
1818
1819 if (ValueA == 0.0 && ValueB == 0.0) {
1820 Hue = 0; // value for grey pixel
1821 } else {
1822 Hue = atan2f(ValueB,ValueA);
1823 }
1824 Col = powf(ValueA * ValueA + ValueB * ValueB, 0.5);
1825
1826 Hue += ((float)Curve->Curve[m_Image[i][0]]-(float)0x7fff)*ScalePi;
1827
1828 m_Image[i][1] = CLIP((int32_t)((cosf(Hue)*Col)+WPH));
1829 m_Image[i][2] = CLIP((int32_t)((sinf(Hue)*Col)+WPH));
1830 }
1831 }
1832 return this;
1833 }
1834
1835 ////////////////////////////////////////////////////////////////////////////////
1836 //
1837 // Apply Saturation Curve
1838 //
1839 ////////////////////////////////////////////////////////////////////////////////
1840
ApplySaturationCurve(const ptCurve * Curve,const short Mode)1841 ptImage* ptImage::ApplySaturationCurve(const ptCurve *Curve,
1842 const short Mode) {
1843
1844 // Best solution would be to use the Lab <-> Lch conversion from lcms.
1845 // This should be faster without sacrificing much quality.
1846
1847 assert (m_ColorSpace == ptSpace_Lab);
1848 // neutral value for a* and b* channel
1849 const float WPH = 0x8080;
1850 const float InvScalePi = 0x7fff / ptPI;
1851
1852 float ValueA = 0.0;
1853 float ValueB = 0.0;
1854
1855 if (Curve->mask() == ptCurve::ChromaMask) { // by chroma
1856 #pragma omp parallel for schedule(static) private(ValueA, ValueB)
1857 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
1858 // Factor by hue
1859 ValueA = (float)m_Image[i][1]-WPH;
1860 ValueB = (float)m_Image[i][2]-WPH;
1861 float Hue = 0;
1862 if (ValueA == 0.0 && ValueB == 0.0) {
1863 Hue = 0; // value for grey pixel
1864 } else {
1865 Hue = atan2f(ValueB,ValueA);
1866 }
1867 while (Hue < 0) Hue += 2.*ptPI;
1868
1869 float Factor = Curve->Curve[CLIP((int32_t)(Hue*InvScalePi))]/(float)0x7fff;
1870 if (Factor == 1.0) continue;
1871 Factor *= Factor;
1872 float m = 0;
1873 if (Mode == 1) {
1874 float Col = powf(ValueA * ValueA + ValueB * ValueB, 0.125);
1875 Col /= 0xd; // normalizing to 0..1
1876
1877 if (Factor > 1)
1878 // work more on desaturated pixels
1879 m = Factor*(1-Col)+Col;
1880 else
1881 // work more on saturated pixels
1882 m = Factor*Col+(1-Col);
1883 } else {
1884 m = Factor;
1885 }
1886 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1.0 - m)));
1887 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1.0 - m)));
1888 }
1889
1890 } else { // by luma
1891 #pragma omp parallel for schedule(static) private(ValueA, ValueB)
1892 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
1893 // Factor by luminance
1894 float Factor = Curve->Curve[m_Image[i][0]]/(float)0x7fff;
1895 if (Factor == 1.0) continue;
1896 Factor *= Factor;
1897 float m = 0;
1898 if (Mode == 1) {
1899 ValueA = (float)m_Image[i][1]-WPH;
1900 ValueB = (float)m_Image[i][2]-WPH;
1901 float Col = powf(ValueA * ValueA + ValueB * ValueB, 0.125);
1902 Col /= 0xd; // normalizing to 0..1
1903
1904 if (Factor > 1)
1905 // work more on desaturated pixels
1906 m = Factor*(1-Col)+Col;
1907 else
1908 // work more on saturated pixels
1909 m = Factor*Col+(1-Col);
1910 } else {
1911 m = Factor;
1912 }
1913 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1.0 - m)));
1914 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1.0 - m)));
1915 }
1916 }
1917 return this;
1918 }
1919
1920 ////////////////////////////////////////////////////////////////////////////////
1921 //
1922 // Apply Texture Curve
1923 //
1924 ////////////////////////////////////////////////////////////////////////////////
1925
ApplyTextureCurve(const ptCurve * Curve,const short Scaling)1926 ptImage* ptImage::ApplyTextureCurve(const ptCurve *Curve, const short Scaling) {
1927 assert (m_ColorSpace == ptSpace_Lab);
1928
1929 const TChannelMask ChannelMask = ChMask_L;
1930
1931 const float Threshold = 10.0/pow(2,Scaling);
1932 const float Softness = 0.01;
1933 const float Opacity = 1.0;
1934 const float EdgeControl = 1.0;
1935
1936 float hValueA = 0.0;
1937 float hValueB = 0.0;
1938 float m = 0.0;
1939
1940 ptImage *ContrastLayer = new ptImage;
1941 ContrastLayer->Set(this);
1942
1943 ContrastLayer->fastBilateralChannel(Threshold, Softness, 2, ChannelMask);
1944
1945 if (Curve->mask() == ptCurve::ChromaMask) {
1946 #pragma omp parallel for schedule(static) private(hValueA, hValueB, m)
1947 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
1948 // Factor by hue
1949 hValueA = ToFloatABNeutral[m_Image[i][1]];
1950 hValueB = ToFloatABNeutral[m_Image[i][2]];
1951
1952 float hHue = ToHue(hValueA, hValueB);
1953 float hCol = powf(ptSqr(hValueA) + ptSqr(hValueB), 0.125f);
1954 hCol /= 0x7; // normalizing to 0..2
1955
1956 float Factor = Curve->Curve[CLIP((int32_t)(hHue/ptPI*ptWPHf))]/(float)0x3fff - 1.0f;
1957 m = 20.0f * Factor * hCol;
1958
1959 ContrastLayer->m_Image[i][0] = CLIP((int32_t) ((ptWPH-(int32_t)ContrastLayer->m_Image[i][0])+m_Image[i][0]));
1960 if (Factor < 0) ContrastLayer->m_Image[i][0] = ToInvertTable[ContrastLayer->m_Image[i][0]];
1961 if (fabsf(Factor*hCol)<0.1f) continue;
1962 ContrastLayer->m_Image[i][0] = Sigmoidal_4_Value(ContrastLayer->m_Image[i][0], m);
1963 }
1964
1965 } else { // by luma
1966 #pragma omp parallel for schedule(static) private(m)
1967 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
1968 // Factor by luminance
1969 float Factor = Curve->Curve[m_Image[i][0]]/(float)0x3fff - 1.0;
1970 m = 20.0f * Factor;
1971
1972 ContrastLayer->m_Image[i][0] = CLIP((int32_t) ((ptWPH-(int32_t)ContrastLayer->m_Image[i][0])+m_Image[i][0]));
1973 if (Factor < 0) ContrastLayer->m_Image[i][0] = ToInvertTable[ContrastLayer->m_Image[i][0]];
1974 if (fabsf(Factor)<0.1f) continue;
1975 ContrastLayer->m_Image[i][0] = Sigmoidal_4_Value(ContrastLayer->m_Image[i][0], m);
1976 }
1977 }
1978
1979 if (EdgeControl)
1980 ContrastLayer->WaveletDenoise(ChannelMask, EdgeControl, 0.2, 0);
1981
1982 Overlay(ContrastLayer->m_Image,Opacity,NULL);
1983
1984 delete ContrastLayer;
1985
1986 return this;
1987 }
1988
1989 ////////////////////////////////////////////////////////////////////////////////
1990 //
1991 // Sigmoidal Contrast
1992 // This can be done with a curve, but this should be a bit more effective.
1993 //
1994 ////////////////////////////////////////////////////////////////////////////////
1995
SigmoidalContrast(const double Contrast,const double Threshold,const short ChannelMask)1996 ptImage* ptImage::SigmoidalContrast(const double Contrast,
1997 const double Threshold,
1998 const short ChannelMask) {
1999
2000 int Channels = 0;
2001 int Channel[3] = {0,1,2};
2002 if (ChannelMask & 1) {Channel[Channels] = 0; Channels++;}
2003 if (ChannelMask & 2) {Channel[Channels] = 1; Channels++;}
2004 if (ChannelMask & 4) {Channel[Channels] = 2; Channels++;}
2005
2006 uint16_t ContrastTable[0x10000];
2007 SigmoidalTable(ContrastTable, Contrast, Threshold);
2008
2009 #pragma omp parallel for default(shared)
2010 for (uint32_t i=0; i < (uint32_t)m_Height*m_Width; i++) {
2011 for (int c = 0; c<Channels; c++)
2012 m_Image[i][Channel[c]] = ContrastTable[ m_Image[i][Channel[c]] ];
2013 }
2014
2015 return this;
2016 }
2017
2018 ////////////////////////////////////////////////////////////////////////////////
2019 //
2020 // Crop
2021 //
2022 ////////////////////////////////////////////////////////////////////////////////
2023
Crop(const uint16_t X,const uint16_t Y,const uint16_t W,const uint16_t H)2024 ptImage* ptImage::Crop(const uint16_t X,
2025 const uint16_t Y,
2026 const uint16_t W,
2027 const uint16_t H) {
2028
2029 assert(m_Colors ==3);
2030 assert( (X+W) <= m_Width);
2031 assert( (Y+H) <= m_Height);
2032
2033 ptImage hCroppedImage;
2034 hCroppedImage.setSize((size_t)W*H);
2035
2036 #pragma omp parallel for
2037 for (uint16_t Row=0;Row<H;Row++) {
2038 for (uint16_t Column=0;Column<W;Column++) {
2039 hCroppedImage.m_Data[Row*W+Column] = m_Data[(Y+Row)*m_Width+X+Column];
2040 }
2041 }
2042
2043 setSize((size_t) W*H);
2044
2045 m_Data = hCroppedImage.m_Data;
2046 m_Width = W;
2047 m_Height = H;
2048
2049 return this;
2050 }
2051
2052 ////////////////////////////////////////////////////////////////////////////////
2053 //
2054 // Overlay
2055 //
2056 ////////////////////////////////////////////////////////////////////////////////
2057
2058 // macros are bad, but for lack of a better idea...
2059 #define Value_4_Amount(AValue) CLIP((int32_t) ((AValue)*Amount + Source*CompAmount))
2060 #define Value_4_Amount_Mask(AValue) CLIP((int32_t) (((AValue)*Mask[i] + Source*(1.0f - Mask[i]))*Amount + Source*CompAmount))
2061 #define LoopBody(ABlock, AResult) \
2062 if (!Mask) { \
2063 for (short Ch=0; Ch<3; Ch++) { \
2064 if (! (ChannelMask & (1<<Ch))) continue; \
2065 _Pragma("omp parallel for default(shared) private(Source, Blend, Multiply, Screen, Overlay, Temp)") \
2066 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) { \
2067 Source = SourceImage[i][Ch]; \
2068 Blend = BlendImage[i][Ch]; \
2069 { ABlock } \
2070 m_Image[i][Ch] = Value_4_Amount(AResult); \
2071 } \
2072 } \
2073 } else { \
2074 for (short Ch=0; Ch<3; Ch++) { \
2075 if (! (ChannelMask & (1<<Ch))) continue; \
2076 _Pragma("omp parallel for default(shared) private(Source, Blend, Multiply, Screen, Overlay, Temp)") \
2077 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) { \
2078 if (Mask[i] == 0.0f) continue; \
2079 Source = SourceImage[i][Ch]; \
2080 Blend = BlendImage[i][Ch]; \
2081 { ABlock } \
2082 m_Image[i][Ch] = Value_4_Amount_Mask(AResult); \
2083 } \
2084 } \
2085 }
2086
Overlay(uint16_t (* OverlayImage)[3],const float Amount,const float * Mask,const TOverlayMode Mode,const bool Swap)2087 ptImage* ptImage::Overlay(
2088 uint16_t (*OverlayImage)[3],
2089 const float Amount,
2090 const float *Mask,
2091 const TOverlayMode Mode,
2092 const bool Swap)
2093 {
2094 const short ChannelMask = (m_ColorSpace == ptSpace_Lab)?1:7;
2095 float Multiply = 0;
2096 float Screen = 0;
2097 float Overlay = 0;
2098 uint16_t Source = 0;
2099 uint16_t Blend = 0;
2100 float Temp = 0;
2101 float CompAmount = 1.0 - Amount;
2102 uint16_t (*SourceImage)[3];
2103 uint16_t (*BlendImage)[3];
2104 if (!Swap) {
2105 SourceImage = m_Image;
2106 BlendImage = OverlayImage;
2107 } else {
2108 BlendImage = m_Image;
2109 SourceImage = OverlayImage;
2110 }
2111
2112 switch (Mode) {
2113 case TOverlayMode::Disabled: // just for completeness
2114 break;
2115
2116 case TOverlayMode::Softlight:
2117 LoopBody({
2118 Multiply = (float)Source*Blend*ptInvWP;
2119 Screen = ptWPf - (float)ToInvertTable[Source]*ToInvertTable[Blend]*ptInvWP;
2120 Overlay = (ToInvertTable[Source]*Multiply + Source*Screen)*ptInvWP;
2121 }, Overlay)
2122 break;
2123
2124 case TOverlayMode::Multiply:
2125 LoopBody({
2126 Multiply = (float)Source*Blend*ptInvWP;
2127 }, Multiply)
2128 break;
2129
2130 case TOverlayMode::Screen:
2131 LoopBody({
2132 Screen = ptWPf - (float)ToInvertTable[Source]*ToInvertTable[Blend]*ptInvWP;
2133 }, Screen)
2134 break;
2135
2136 case TOverlayMode::GammaDark:
2137 LoopBody({
2138 if (Blend == 0) Multiply = 0;
2139 else Multiply = ptWPf*powf(Source*ptInvWP,ptWPf/Blend);
2140 }, Multiply)
2141 break;
2142
2143 case TOverlayMode::GammaBright:
2144 LoopBody({
2145 if (Blend == ptWP) Multiply = ptWPf;
2146 else Multiply = ptWPf-ptWPf*powf(ToInvertTable[Source]*ptInvWP,ptWPf/(float)ToInvertTable[Blend]);
2147 }, Multiply)
2148 break;
2149
2150 case TOverlayMode::Normal:
2151 LoopBody({
2152 }, Blend)
2153 break;
2154
2155 case TOverlayMode::Lighten:
2156 if (!Mask) {
2157 for (short Ch=0; Ch<3; Ch++) {
2158 // Is it a channel we are supposed to handle ?
2159 if (! (ChannelMask & (1<<Ch))) continue;
2160 #pragma omp parallel for default(shared) private(Source, Blend, Multiply, Screen, Overlay, Temp)
2161 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2162 Source = SourceImage[i][Ch];
2163 Blend = BlendImage[i][Ch];
2164 m_Image[i][Ch] = CLIP((int32_t) (ptMax((uint16_t)(Blend*Amount), Source)));
2165 }
2166 }
2167 } else {
2168 for (short Ch=0; Ch<3; Ch++) {
2169 // Is it a channel we are supposed to handle ?
2170 if (! (ChannelMask & (1<<Ch))) continue;
2171 #pragma omp parallel for default(shared) private(Source, Blend, Multiply, Screen, Overlay, Temp)
2172 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2173 if (Mask[i] == 0.0f) continue;
2174 Source = SourceImage[i][Ch];
2175 Blend = BlendImage[i][Ch];
2176 m_Image[i][Ch] = CLIP((int32_t) (ptMax((uint16_t)(Blend*Mask[i] + Source*(1.0f-Mask[i])*Amount),Source)));
2177 }
2178 }
2179 }
2180 break;
2181
2182 case TOverlayMode::Darken:
2183 if (!Mask) {
2184 for (short Ch=0; Ch<3; Ch++) {
2185 // Is it a channel we are supposed to handle ?
2186 if (! (ChannelMask & (1<<Ch))) continue;
2187 #pragma omp parallel for default(shared) private(Source, Blend, Multiply, Screen, Overlay, Temp)
2188 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2189 Source = SourceImage[i][Ch];
2190 Blend = BlendImage[i][Ch];
2191 m_Image[i][Ch] = CLIP((int32_t) (ptMin((uint16_t)(Blend*Amount), Source)));
2192 }
2193 }
2194 } else {
2195 for (short Ch=0; Ch<3; Ch++) {
2196 // Is it a channel we are supposed to handle ?
2197 if (! (ChannelMask & (1<<Ch))) continue;
2198 #pragma omp parallel for default(shared) private(Source, Blend, Multiply, Screen, Overlay, Temp)
2199 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2200 if (Mask[i] == 0.0f) continue;
2201 Source = SourceImage[i][Ch];
2202 Blend = BlendImage[i][Ch];
2203 m_Image[i][Ch] = CLIP((int32_t) (ptMin((uint16_t)(Blend*Mask[i]+Source*(1.0f-Mask[i])*Amount),Source)));
2204 }
2205 }
2206 }
2207 break;
2208
2209 case TOverlayMode::Overlay:
2210 LoopBody({
2211 if (Source <= ptWPH) Overlay = Source*Blend*ptInvWP;
2212 else Overlay = ptWPf - ToInvertTable[Source]*ToInvertTable[Blend]*ptInvWP;
2213 }, Overlay)
2214 break;
2215
2216 case TOverlayMode::GrainMerge:
2217 LoopBody({
2218 }, (float)Blend + Source - ptWPHf)
2219 break;
2220
2221 case TOverlayMode::ColorDodge: // a/(1-b)
2222 LoopBody({
2223 if (Source == 0) Temp = 0;
2224 else {
2225 if (Blend == ptWP) Temp = ptWPf;
2226 else Temp = (float)Source / (1.0f - (float)Blend*ptInvWP);
2227 }
2228 }, Temp)
2229 break;
2230
2231 case TOverlayMode::ColorBurn: // 1-(1-a)/b
2232 LoopBody({
2233 if (Source == ptWP) Temp = ptWPf;
2234 else {
2235 if (Blend == 0) Temp = 0;
2236 else Temp = ptWPf - ( ToInvertTable[Source] / (Blend*ptInvWP));
2237 }
2238 }, Temp)
2239 break;
2240
2241 case TOverlayMode::ShowMask:
2242 if (Mask) {
2243 for (short Ch=0; Ch<3; Ch++) {
2244 // Is it a channel we are supposed to handle ?
2245 if (! (ChannelMask & (1<<Ch))) continue;
2246 #pragma omp parallel for default(shared) private(Source, Blend, Multiply, Screen, Overlay, Temp)
2247 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2248 m_Image[i][Ch] = CLIP((int32_t) (Mask[i]*ptWPf));
2249 }
2250 }
2251 }
2252 break;
2253
2254 case TOverlayMode::Replace: // Replace, just for testing
2255 for (short Ch=0; Ch<3; Ch++) {
2256 // Is it a channel we are supposed to handle ?
2257 if (! (ChannelMask & (1<<Ch))) continue;
2258 #pragma omp parallel for default(shared) private(Source, Blend, Multiply, Screen, Overlay, Temp)
2259 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2260 Blend = BlendImage[i][Ch];
2261 m_Image[i][Ch] = CLIP((int32_t) Blend);
2262 }
2263 }
2264 break;
2265
2266 }
2267 return this;
2268 }
2269 #undef Value_4_Amount
2270 #undef Value_4_Amount_Mask
2271 #undef LoopBody
2272
2273 ////////////////////////////////////////////////////////////////////////////////
2274 //
2275 // Flip
2276 //
2277 ////////////////////////////////////////////////////////////////////////////////
2278
Flip(const short FlipMode)2279 ptImage* ptImage::Flip(const short FlipMode) {
2280
2281 uint16_t Width = m_Width;
2282 uint16_t Height = m_Height;
2283 uint16_t Cache = 0;
2284 if (FlipMode == ptFlipMode_Vertical) {
2285 #pragma omp parallel for default(shared) private(Cache)
2286 for (uint16_t Row=0;Row<Height/2;Row++) {
2287 for (uint16_t Column=0;Column<Width;Column++) {
2288 Cache = m_Image[Row*Width+Column][0];
2289 m_Image[Row*Width+Column][0] = m_Image[(Height-1-Row)*Width+Column][0];
2290 m_Image[(Height-1-Row)*Width+Column][0] = Cache;
2291 Cache = m_Image[Row*Width+Column][1];
2292 m_Image[Row*Width+Column][1] = m_Image[(Height-1-Row)*Width+Column][1];
2293 m_Image[(Height-1-Row)*Width+Column][1] = Cache;
2294 Cache = m_Image[Row*Width+Column][2];
2295 m_Image[Row*Width+Column][2] = m_Image[(Height-1-Row)*Width+Column][2];
2296 m_Image[(Height-1-Row)*Width+Column][2] = Cache;
2297 }
2298 }
2299 } else if (FlipMode == ptFlipMode_Horizontal) {
2300 #pragma omp parallel for default(shared) private(Cache)
2301 for (uint16_t Row=0;Row<Height;Row++) {
2302 for (uint16_t Column=0;Column<Width/2;Column++) {
2303 Cache = m_Image[Row*Width+Column][0];
2304 m_Image[Row*Width+Column][0] = m_Image[Row*Width+(Width-1-Column)][0];
2305 m_Image[Row*Width+(Width-1-Column)][0] = Cache;
2306 Cache = m_Image[Row*Width+Column][1];
2307 m_Image[Row*Width+Column][1] = m_Image[Row*Width+(Width-1-Column)][1];
2308 m_Image[Row*Width+(Width-1-Column)][1] = Cache;
2309 Cache = m_Image[Row*Width+Column][2];
2310 m_Image[Row*Width+Column][2] = m_Image[Row*Width+(Width-1-Column)][2];
2311 m_Image[Row*Width+(Width-1-Column)][2] = Cache;
2312 }
2313 }
2314 }
2315 return this;
2316 }
2317
2318
2319 ////////////////////////////////////////////////////////////////////////////////
2320 //
2321 // Levels
2322 //
2323 ////////////////////////////////////////////////////////////////////////////////
2324
Levels(const float BlackPoint,const float WhitePoint)2325 ptImage* ptImage::Levels(const float BlackPoint,
2326 const float WhitePoint) {
2327
2328 const float WP = 0xffff;
2329 const short NrChannels = (m_ColorSpace == ptSpace_Lab)?1:3;
2330
2331 if (fabs(BlackPoint-WhitePoint)>0.001) {
2332 float m = 1.0/(WhitePoint-BlackPoint);
2333 float t = -BlackPoint/(WhitePoint-BlackPoint)*WP;
2334 #pragma omp parallel for
2335 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2336 for (short Ch=0; Ch<NrChannels; Ch++) {
2337 m_Image[i][Ch] = CLIP((int32_t)(m_Image[i][Ch]*m+t));
2338 }
2339 }
2340 }
2341 return this;
2342 }
2343
2344 ////////////////////////////////////////////////////////////////////////////////
2345 //
2346 // DeFringe
2347 // Original implementation by Emil Martinec for RawTherapee
2348 //
2349 ////////////////////////////////////////////////////////////////////////////////
2350
DeFringe(const double Radius,const short Threshold,const int Flags,const double Shift)2351 ptImage* ptImage::DeFringe(const double Radius,
2352 const short Threshold,
2353 const int Flags,
2354 const double Shift) {
2355
2356 assert (m_ColorSpace == ptSpace_Lab);
2357
2358 int Neighborhood = ceil(2*Radius)+1;
2359
2360 float (*ChromaDiff) = (float(*)) CALLOC(m_Width*m_Height,sizeof(*ChromaDiff));
2361 ptMemoryError(ChromaDiff,__FILE__,__LINE__);
2362
2363 ptImage *SaveLayer = new ptImage;
2364 SaveLayer->Set(this);
2365
2366 ptCIBlur(Radius, 6);
2367
2368 uint32_t Size = m_Width * m_Height;
2369
2370 float Average = 0.0f;
2371 #pragma omp parallel for schedule(static) reduction (+:Average)
2372 for (uint32_t i = 0; i < Size; i++) {
2373 ChromaDiff[i] = SQR((float)m_Image[i][1]-(float)SaveLayer->m_Image[i][1]) +
2374 SQR((float)m_Image[i][2]-(float)SaveLayer->m_Image[i][2]);
2375 Average += ChromaDiff[i];
2376 }
2377 Average /= Size;
2378
2379 float NewThreshold = Threshold*Average/33.0f;
2380 float HueShift = Shift * ptPI/6;
2381 float Val1 = MAX(0.f,HueShift);
2382 float Val2 = ptPI/3+HueShift;
2383 float Val3 = 2*ptPI/3+HueShift;
2384 float Val4 = 3*ptPI/3+HueShift;
2385 float Val5 = 4*ptPI/3+HueShift;
2386 float Val6 = 5*ptPI/3+HueShift;
2387 float Val7 = 6*ptPI/3+HueShift;
2388
2389 #pragma omp parallel for schedule(dynamic)
2390 for (uint16_t Row = 0; Row < m_Height; Row++) {
2391 for (uint16_t Col = 0; Col < m_Width; Col++) {
2392 short CorrectPixel = 0;
2393 uint32_t Index = Row*m_Width+Col;
2394 if (ChromaDiff[Index] > NewThreshold) {
2395 // Calculate hue.
2396 float ValueA = (float)SaveLayer->m_Image[Index][1]-0x8080;
2397 float ValueB = (float)SaveLayer->m_Image[Index][2]-0x8080;
2398 float Hue = 0;
2399 if (ValueA == 0.0 && ValueB == 0.0) {
2400 Hue = 0; // value for grey pixel
2401 } else {
2402 Hue = atan2f(ValueB,ValueA);
2403 }
2404 while (Hue < 0) Hue += 2.*ptPI;
2405 // Check if we want to treat that hue
2406 if (Flags & 1) // red
2407 if ((Hue >= Val1 && Hue < Val2) || Hue >= Val7) CorrectPixel = 1;
2408 if (Flags & 2) // yellow
2409 if (Hue >= Val2 && Hue < Val3) CorrectPixel = 1;
2410 if (Flags & 4) // green
2411 if (Hue >= Val3 && Hue < Val4) CorrectPixel = 1;
2412 if (Flags & 8) // cyan
2413 if (Hue >= Val4 && Hue < Val5) CorrectPixel = 1;
2414 if (Flags & 16) // blue
2415 if (Hue >= Val5 && Hue < Val6) CorrectPixel = 1;
2416 if (Flags & 32) // purple
2417 if ((Hue >= Val6 && Hue < Val7) || Hue < Val1) CorrectPixel = 1;
2418 }
2419 if (CorrectPixel == 1 ) {
2420 float TotalA=0;
2421 float TotalB=0;
2422 float Total=0;
2423 float Weight;
2424 for (int i1 = MAX(0,Row-Neighborhood+1); i1 < MIN((int)m_Height,Row+Neighborhood); i1++)
2425 for (int j1 = MAX(0,Col-Neighborhood+1); j1 < MIN((int)m_Width,Col+Neighborhood); j1++) {
2426 // Neighborhood average of pixels weighted by chrominance
2427 uint32_t Index2 = i1*m_Width+j1;
2428 Weight = 1/(ChromaDiff[Index2]+Average);
2429 TotalA += Weight*SaveLayer->m_Image[Index2][1];
2430 TotalB += Weight*SaveLayer->m_Image[Index2][2];
2431 Total += Weight;
2432 }
2433 m_Image[Index][1] = CLIP((int32_t)(TotalA/Total));
2434 m_Image[Index][2] = CLIP((int32_t)(TotalB/Total));
2435 } else {
2436 m_Image[Index][1] = SaveLayer->m_Image[Index][1];
2437 m_Image[Index][2] = SaveLayer->m_Image[Index][2];
2438 }
2439 }
2440 }
2441
2442 delete SaveLayer;
2443 FREE(ChromaDiff);
2444
2445 return this;
2446 }
2447
2448 ////////////////////////////////////////////////////////////////////////////////
2449 //
2450 // Impluse noise reduction
2451 // Original implementation by Emil Martinec for RawTherapee
2452 //
2453 ////////////////////////////////////////////////////////////////////////////////
2454
DenoiseImpulse(const double ThresholdL,const double ThresholdAB)2455 ptImage* ptImage::DenoiseImpulse(const double ThresholdL,
2456 const double ThresholdAB) {
2457
2458 assert (m_ColorSpace == ptSpace_Lab);
2459
2460 float hpfabs, hfnbrave;
2461 float hpfabs1, hfnbrave1, hpfabs2, hfnbrave2;
2462
2463 ptImage *LowPass = new ptImage;
2464 LowPass->Set(this);
2465
2466 short ChannelMask = 0;
2467 if (ThresholdL != 0.0) ChannelMask += 1;
2468 if (ThresholdAB != 0.0) ChannelMask += 6;
2469 LowPass->ptCIBlur(2.0f, ChannelMask);
2470
2471 short (*Impulse)[3] = (short (*)[3]) CALLOC(m_Width*m_Height,sizeof(*Impulse));
2472 ptMemoryError(Impulse,__FILE__,__LINE__);
2473
2474 static float eps = 1.0;
2475 float wtdsum, dirwt, norm;
2476
2477 uint32_t Index = 0;
2478 uint32_t Index1 = 0;
2479
2480 if (ThresholdL != 0.0) {
2481 #pragma omp parallel for schedule(static) private(hpfabs, hfnbrave, Index, Index1)
2482 for (uint16_t Row = 0; Row < m_Height; Row++) {
2483 for (uint16_t Col = 0; Col < m_Width; Col++) {
2484 Index = Row*m_Width+Col;
2485 hpfabs = fabs((int32_t)m_Image[Index][0] - (int32_t)LowPass->m_Image[Index][0]);
2486 hfnbrave = 0;
2487 //block average of high pass data
2488 for (uint16_t i1 = MAX(0,Row-2); i1 <= MIN(Row+2,m_Height-1); i1++ ) {
2489 for (uint16_t j1 = MAX(0,Col-2); j1 <= MIN(Col+2,m_Width-1); j1++ ) {
2490 Index1 = i1*m_Width+j1;
2491 hfnbrave += fabs((int32_t)m_Image[Index1][0] - (int32_t)LowPass->m_Image[Index1][0]);
2492 }
2493 }
2494 hfnbrave = (hfnbrave - hpfabs) / 24.0f;
2495 hpfabs > (hfnbrave*(5.5-ThresholdL)) ? Impulse[Index][0]=1 : Impulse[Index][0]=0;
2496 }//now impulsive values have been identified
2497 }
2498 }
2499
2500 if (ThresholdAB != 0.0) {
2501 #pragma omp parallel for schedule(static) private(hpfabs1, hfnbrave1, hpfabs2, hfnbrave2, Index, Index1)
2502 for (uint16_t Row = 0; Row < m_Height; Row++) {
2503 for (uint16_t Col = 0; Col < m_Width; Col++) {
2504 Index = Row*m_Width+Col;
2505 hpfabs1 = fabs((int32_t)m_Image[Index][1] - (int32_t)LowPass->m_Image[Index][1]);
2506 hpfabs2 = fabs((int32_t)m_Image[Index][2] - (int32_t)LowPass->m_Image[Index][2]);
2507 hfnbrave1 = 0;
2508 hfnbrave2 = 0;
2509 //block average of high pass data
2510 for (uint16_t i1 = MAX(0,Row-2); i1 <= MIN(Row+2,m_Height-1); i1++ ) {
2511 for (uint16_t j1 = MAX(0,Col-2); j1 <= MIN(Col+2,m_Width-1); j1++ ) {
2512 Index1 = i1*m_Width+j1;
2513 hfnbrave1 += fabs((int32_t)m_Image[Index1][1] - (int32_t)LowPass->m_Image[Index1][1]);
2514 hfnbrave2 += fabs((int32_t)m_Image[Index1][2] - (int32_t)LowPass->m_Image[Index1][2]);
2515 }
2516 }
2517 hfnbrave1 = (hfnbrave1 - hpfabs1) / 24.0f;
2518 hfnbrave2 = (hfnbrave2 - hpfabs2) / 24.0f;
2519 hpfabs1 > (hfnbrave1*(5.5-ThresholdAB)) ? Impulse[Index][1]=1 : Impulse[Index][1]=0;
2520 hpfabs2 > (hfnbrave2*(5.5-ThresholdAB)) ? Impulse[Index][2]=1 : Impulse[Index][2]=0;
2521 }//now impulsive values have been identified
2522 }
2523 }
2524
2525 for (short Channel=0;Channel<3;Channel++) {
2526 // Is it a channel we are supposed to handle ?
2527 if (! (ChannelMask & (1<<Channel))) continue;
2528 #pragma omp parallel for schedule(static) private(norm, wtdsum, dirwt, Index, Index1)
2529 for (uint16_t Row = 0; Row < m_Height; Row++) {
2530 for (uint16_t Col = 0; Col < m_Width; Col++) {
2531 Index = Row*m_Width+Col;
2532 if (!Impulse[Index][Channel]) continue;
2533 norm=0.0;
2534 wtdsum=0.0;
2535 for (uint16_t i1 = MAX(0,Row-2); i1 <= MIN(Row+2,m_Height-1); i1++ ) {
2536 for (uint16_t j1 = MAX(0,Col-2); j1 <= MIN(Col+2,m_Width-1); j1++ ) {
2537 if (i1==Row && j1==Col) continue;
2538 Index1 = i1*m_Width+j1;
2539 if (Impulse[Index1][Channel]) continue;
2540 float Temp = (float)m_Image[Index1][Channel]-(float)m_Image[Index][Channel];
2541 dirwt = 1/(SQR(Temp)+eps);//use more sophisticated rangefn???
2542 wtdsum += dirwt*m_Image[Index1][Channel];
2543 norm += dirwt;
2544 }
2545 }
2546 //wtdsum /= norm;
2547 if (norm) {
2548 m_Image[Index][Channel] = CLIP((int32_t) (wtdsum/norm));//low pass filter
2549 }
2550 }
2551 }//now impulsive values have been corrected
2552 }
2553
2554
2555 FREE(Impulse);
2556 delete LowPass;
2557
2558 return this;
2559 }
2560
2561 ////////////////////////////////////////////////////////////////////////////////
2562 //
2563 // Reinhard 05 tone mapping operator
2564 // Adapted from PFSTMO package, original Copyright (C) 2007 Grzegorz Krawczyk
2565 //
2566 ////////////////////////////////////////////////////////////////////////////////
2567
Reinhard05(const float Brightness,const float Chromatic,const float Light)2568 ptImage* ptImage::Reinhard05(const float Brightness,
2569 const float Chromatic,
2570 const float Light)
2571 {
2572 assert (m_ColorSpace != ptSpace_Lab);
2573 // OpenMP is done for RGB
2574 // const short NrChannels = (m_ColorSpace == ptSpace_Lab)?1:3;
2575 const short NrChannels = 3;
2576
2577 const float DChromatic = 1 - Chromatic;
2578 const float DLight = 1 - Light;
2579
2580 uint32_t m_Size = m_Width * m_Height;
2581 float (*Temp)[3] = (float (*)[3]) CALLOC(m_Size,sizeof(*Temp));
2582 ptMemoryError(Temp,__FILE__,__LINE__);
2583
2584 // luminance
2585 float (*Y) = (float (*)) CALLOC(m_Width*m_Height,sizeof(*Y));
2586 ptMemoryError(Y,__FILE__,__LINE__);
2587
2588 if (NrChannels == 1)
2589 #pragma omp parallel for schedule(static)
2590 for (uint32_t i=0; i < m_Size; i++) {
2591 Y[i] = ToFloatTable[m_Image[i][0]];
2592 }
2593 else
2594 #pragma omp parallel for schedule(static)
2595 for (uint32_t i=0; i < m_Size; i++) {
2596 Y[i] = ToFloatTable[RGB_2_L(m_Image[i])];
2597 }
2598
2599 float max_lum = 0.0f;
2600 float min_lum = 0.0f;
2601 float world_lum = 0.0f;
2602 float Cav[] = {0.0f, 0.0f, 0.0f};
2603 float Cav1 = 0.0f;
2604 float Cav2 = 0.0f;
2605 float Cav3 = 0.0f;
2606 float Lav = 0.0f;
2607
2608 #pragma omp parallel
2609 {
2610 float thread_max = 0.0;
2611 float thread_min = 0.0;
2612 #pragma omp for reduction(+:world_lum,Cav1,Cav2,Cav3,Lav)
2613 for (uint32_t i=0; i < m_Size; i++) {
2614 float lum = Y[i];
2615 thread_max = (thread_max > lum) ? thread_max : lum;
2616 thread_min = (thread_min < lum) ? thread_min : lum;
2617 world_lum += logf(2.3e-5+lum);
2618 Cav1 += m_Image[i][0];
2619 Cav2 += m_Image[i][1];
2620 Cav3 += m_Image[i][2];
2621 Lav += lum;
2622 }
2623 #pragma omp critical
2624 if (thread_max > max_lum) {
2625 max_lum = thread_max;
2626 }
2627 #pragma omp critical
2628 if (thread_min < min_lum) {
2629 min_lum = thread_min;
2630 }
2631 }
2632 Cav[0] = Cav1;
2633 Cav[1] = Cav2;
2634 Cav[2] = Cav3;
2635 world_lum /= (float)m_Size;
2636 for (int c = 0; c < NrChannels; c++)
2637 Cav[c] = Cav[c] / ((float)m_Size * (float)0xffff);
2638 Lav /= (float)m_Size;
2639
2640 //--- tone map image
2641 max_lum = logf( max_lum );
2642 min_lum = logf( min_lum );
2643
2644 // image key
2645 float k = (max_lum - world_lum) / (max_lum - min_lum);
2646 // image contrast based on key value
2647 float m = 0.3f+0.7f*powf(k,1.4f);
2648 // image brightness
2649 float f = expf(-Brightness);
2650
2651 float max_col = 0.0f;
2652 float min_col = 1.0f;
2653
2654 #pragma omp parallel
2655 {
2656 float thread_max = 0.0f;
2657 float thread_min = 1.0f;
2658 #pragma omp for
2659 for (uint32_t i=0; i < m_Size; i++) {
2660 float l = Y[i];
2661 float col;
2662 if( l != 0.0f ) {
2663 for (int c = 0; c < NrChannels; c++) {
2664 col = ToFloatTable[m_Image[i][c]];
2665
2666 if(col != 0.0f) {
2667 // local light adaptation
2668 float Il = Chromatic * col + DChromatic * l;
2669 // global light adaptation
2670 float Ig = Chromatic * Cav[c] + DChromatic * Lav;
2671 // interpolated light adaptation
2672 float Ia = Light*Il + DLight * Ig;
2673 // photoreceptor equation
2674 col /= col + powf(f*Ia, m);
2675 }
2676
2677 thread_max = (col>thread_max) ? col : thread_max;
2678 thread_min = (col<thread_min) ? col : thread_min;
2679
2680 Temp[i][c]=col;
2681 }
2682 }
2683 }
2684 #pragma omp critical
2685 if (thread_max > max_col) {
2686 max_col = thread_max;
2687 }
2688 #pragma omp critical
2689 if (thread_min < min_col) {
2690 min_col = thread_min;
2691 }
2692 }
2693
2694 //--- normalize intensities
2695 #pragma omp parallel for schedule(static)
2696 for (uint32_t i=0; i < m_Size; i++) {
2697 for (int c = 0; c < NrChannels; c++)
2698 m_Image[i][c] = CLIP((int32_t) ((Temp[i][c]-min_col)/(max_col-min_col)*0xffff));
2699 }
2700
2701 FREE (Temp);
2702 FREE (Y);
2703
2704 return this;
2705 }
2706
2707 //------------------------------------------------------------------------------
2708
ColorIntensity(int AVibrance,int ARed,int AGreen,int ABlue)2709 ptImage* ptImage::ColorIntensity(int AVibrance, int ARed, int AGreen, int ABlue) {
2710 TChannelMatrix hMixer;
2711
2712 if (AVibrance != 0) {
2713 hMixer[0][0] = 1.0 + (AVibrance/150.0);
2714 hMixer[0][1] = -(AVibrance/300.0);
2715 hMixer[0][2] = hMixer[0][1];
2716 hMixer[1][0] = hMixer[0][1];
2717 hMixer[1][1] = hMixer[0][0];
2718 hMixer[1][2] = hMixer[0][1];
2719 hMixer[2][0] = hMixer[0][1];
2720 hMixer[2][1] = hMixer[0][1];
2721 hMixer[2][2] = hMixer[0][0];
2722
2723 this->mixChannels(hMixer);
2724 }
2725
2726 if ((ARed != 0) || (AGreen != 0) || (ABlue != 0)) {
2727 hMixer[0][0] = 1.0 + (ARed/150.0);
2728 hMixer[0][1] = -(ARed/300.0);
2729 hMixer[0][2] = hMixer[0][1];
2730 hMixer[1][0] = -(AGreen/300.0);
2731 hMixer[1][1] = 1.0+(AGreen/150.0);;
2732 hMixer[1][2] = hMixer[1][0];
2733 hMixer[2][0] = -(ABlue/300.0);
2734 hMixer[2][1] = hMixer[2][0];
2735 hMixer[2][2] = 1.0+(ABlue/150.0);
2736
2737 this->mixChannels(hMixer);
2738 }
2739
2740 return this;
2741 }
2742
2743 ////////////////////////////////////////////////////////////////////////////////
2744 //
2745 // Color Boost
2746 //
2747 ////////////////////////////////////////////////////////////////////////////////
2748
ColorBoost(const double ValueA,const double ValueB)2749 ptImage* ptImage::ColorBoost(const double ValueA,
2750 const double ValueB) {
2751
2752 assert ((m_ColorSpace == ptSpace_Lab));
2753
2754 const double WPH = 0x8080; // Neutral in Lab A and B
2755
2756 double t1 = (1-ValueA)*WPH;
2757 double t2 = (1-ValueB)*WPH;
2758 if (ValueA!=1.0) {
2759 #pragma omp parallel for
2760 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2761 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1]*ValueA+t1));
2762 }
2763 }
2764 if (ValueB!=1.0) {
2765 #pragma omp parallel for
2766 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2767 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2]*ValueB+t2));
2768 }
2769 }
2770 return this;
2771 }
2772
2773 ////////////////////////////////////////////////////////////////////////////////
2774 //
2775 // LCH conversion and L adjustment
2776 //
2777 ////////////////////////////////////////////////////////////////////////////////
2778
LumaAdjust(const double LC1,const double LC2,const double LC3,const double LC4,const double LC5,const double LC6,const double LC7,const double LC8)2779 ptImage* ptImage::LumaAdjust(const double LC1, // 8 colors for L
2780 const double LC2,
2781 const double LC3,
2782 const double LC4,
2783 const double LC5,
2784 const double LC6,
2785 const double LC7,
2786 const double LC8)
2787 {
2788 assert (m_ColorSpace == ptSpace_Lab);
2789 float WPH = 0x7fff;
2790 float IQPI = 4/ptPI;
2791
2792 #pragma omp parallel for schedule(static)
2793 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
2794 float Col = powf(((float)m_Image[i][1]-WPH)*((float)m_Image[i][1]-WPH) +
2795 ((float)m_Image[i][2]-WPH)*((float)m_Image[i][2]-WPH), 0.25);
2796 Col /= 0xb5; // normalizing to 0..1, sqrt(0x7fff)
2797 float Hue = 0;
2798 if (m_Image[i][1] == WPH && m_Image[i][2] == WPH) {
2799 Hue = 0; // value for grey pixel
2800 } else {
2801 Hue = atan2f((float)m_Image[i][2]-WPH,
2802 (float)m_Image[i][1]-WPH);
2803 }
2804 while (Hue < 0) Hue += 2.*ptPI;
2805
2806 if ( LC1 != 0 && Hue > -.1 && Hue < ptPI/4)
2807 m_Image[i][0] = CLIP((int32_t)(m_Image[i][0] * powf(2,(1.-fabsf(Hue-0)*IQPI)*LC1*Col)));
2808 if ( LC2 != 0 && Hue > 0 && Hue < ptPI/2)
2809 m_Image[i][0] = CLIP((int32_t)(m_Image[i][0] * powf(2,(1.-fabsf(Hue-ptPI/4)*IQPI)*LC2*Col)));
2810 if ( LC3 != 0 && Hue > ptPI/4 && Hue < ptPI*3/4)
2811 m_Image[i][0] = CLIP((int32_t)(m_Image[i][0] * powf(2,(1.-fabsf(Hue-ptPI/2)*IQPI)*LC3*Col)));
2812 if ( LC4 != 0 && Hue > ptPI/2 && Hue < ptPI)
2813 m_Image[i][0] = CLIP((int32_t)(m_Image[i][0] * powf(2,(1.-fabsf(Hue-ptPI*3/4)*IQPI)*LC4*Col)));
2814 if ( LC5 != 0 && Hue > ptPI*3/4 && Hue < ptPI*5/4)
2815 m_Image[i][0] = CLIP((int32_t)(m_Image[i][0] * powf(2,(1.-fabsf(Hue-ptPI)*IQPI)*LC5*Col)));
2816 if ( LC6 != 0 && Hue > ptPI && Hue < ptPI*6/4)
2817 m_Image[i][0] = CLIP((int32_t)(m_Image[i][0] * powf(2,(1.-fabsf(Hue-ptPI*5/4)*IQPI)*LC6*Col)));
2818 if ( LC7 != 0 && Hue > ptPI*5/4 && Hue < ptPI*7/4)
2819 m_Image[i][0] = CLIP((int32_t)(m_Image[i][0] * powf(2,(1.-fabsf(Hue-ptPI*6/4)*IQPI)*LC7*Col)));
2820 if ( LC8 != 0 && Hue > ptPI*6/4 && Hue < ptPI*8/4)
2821 m_Image[i][0] = CLIP((int32_t)(m_Image[i][0] * powf(2,(1.-fabsf(Hue-ptPI*7/4)*IQPI)*LC8*Col)));
2822 if ( LC1 != 0 && Hue > ptPI*7/4 && Hue < ptPI*2.1)
2823 m_Image[i][0] = CLIP((int32_t)(m_Image[i][0] * powf(2,(1.-fabsf(Hue-ptPI*2)*IQPI)*LC1*Col)));
2824 }
2825
2826 return this;
2827 }
2828
2829 //==============================================================================
2830
SatAdjust(const double SC1,const double SC2,const double SC3,const double SC4,const double SC5,const double SC6,const double SC7,const double SC8)2831 ptImage* ptImage::SatAdjust(const double SC1, // 8 colors for saturation
2832 const double SC2,
2833 const double SC3,
2834 const double SC4,
2835 const double SC5,
2836 const double SC6,
2837 const double SC7,
2838 const double SC8)
2839 {
2840 assert (m_ColorSpace == ptSpace_Lab);
2841 float WPH = 0x7fff;
2842 float IQPI = 4/ptPI;
2843
2844 #pragma omp parallel for schedule(static)
2845 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
2846 float Col = powf(((float)m_Image[i][1]-WPH)*((float)m_Image[i][1]-WPH) +
2847 ((float)m_Image[i][2]-WPH)*((float)m_Image[i][2]-WPH), 0.25);
2848 Col /= 0xb5; // normalizing to 0..1, sqrt(0x7fff)
2849 float Hue = 0;
2850 if (m_Image[i][1] == WPH && m_Image[i][2] == WPH) {
2851 Hue = 0; // value for grey pixel
2852 } else {
2853 Hue = atan2f((float)m_Image[i][2]-WPH,
2854 (float)m_Image[i][1]-WPH);
2855 }
2856 while (Hue < 0) Hue += 2.*ptPI;
2857
2858 float m = 0;
2859 if ( SC1 != 0 && Hue > -.1 && Hue < ptPI/4) {
2860 m = powf(8,(1.-fabsf(Hue-0)*IQPI)*SC1*Col);
2861 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1. - m)));
2862 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1. - m)));
2863 }
2864 if ( SC2 != 0 && Hue > 0 && Hue < ptPI/2) {
2865 m = powf(8,(1.-fabsf(Hue-ptPI/4)*IQPI)*SC2*Col);
2866 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1. - m)));
2867 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1. - m)));
2868 }
2869 if ( SC3 != 0 && Hue > ptPI/4 && Hue < ptPI*3/4) {
2870 m = powf(8,(1.-fabsf(Hue-ptPI/2)*IQPI)*SC3*Col);
2871 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1. - m)));
2872 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1. - m)));
2873 }
2874 if ( SC4 != 0 && Hue > ptPI/2 && Hue < ptPI) {
2875 m = powf(8,(1.-fabsf(Hue-ptPI*3/4)*IQPI)*SC4*Col);
2876 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1. - m)));
2877 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1. - m)));
2878 }
2879 if ( SC5 != 0 && Hue > ptPI*3/4 && Hue < ptPI*5/4) {
2880 m = powf(8,(1.-fabsf(Hue-ptPI)*IQPI)*SC5*Col);
2881 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1. - m)));
2882 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1. - m)));
2883 }
2884 if ( SC6 != 0 && Hue > ptPI && Hue < ptPI*6/4) {
2885 m = powf(8,(1.-fabsf(Hue-ptPI*5/4)*IQPI)*SC6*Col);
2886 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1. - m)));
2887 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1. - m)));
2888 }
2889 if ( SC7 != 0 && Hue > ptPI*5/4 && Hue < ptPI*7/4) {
2890 m = powf(8,(1.-fabsf(Hue-ptPI*6/4)*IQPI)*SC7*Col);
2891 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1. - m)));
2892 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1. - m)));
2893 }
2894 if ( SC8 != 0 && Hue > ptPI*6/4 && Hue < ptPI*8/4) {
2895 m = powf(8,(1.-fabsf(Hue-ptPI*7/4)*IQPI)*SC8*Col);
2896 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1. - m)));
2897 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1. - m)));
2898 }
2899 if ( SC1 != 0 && Hue > ptPI*7/4 && Hue < ptPI*2.1) {
2900 m = powf(8,(1.-fabsf(Hue-ptPI*2)*IQPI)*SC1*Col);
2901 m_Image[i][1] = CLIP((int32_t)(m_Image[i][1] * m + WPH * (1. - m)));
2902 m_Image[i][2] = CLIP((int32_t)(m_Image[i][2] * m + WPH * (1. - m)));
2903 }
2904 }
2905
2906 return this;
2907 }
2908
2909 ////////////////////////////////////////////////////////////////////////////////
2910 //
2911 // Outline
2912 //
2913 ////////////////////////////////////////////////////////////////////////////////
2914
Outline(const TOverlayMode Mode,const short GradientMode,const ptCurve * Curve,const double Weight,const double Radius,const bool SwitchLayer)2915 ptImage* ptImage::Outline(
2916 const TOverlayMode Mode,
2917 const short GradientMode,
2918 const ptCurve *Curve,
2919 const double Weight,
2920 const double Radius,
2921 const bool SwitchLayer)
2922 {
2923 assert (m_ColorSpace == ptSpace_Lab);
2924
2925 if (Mode == TOverlayMode::Disabled) return this;
2926
2927 ptImage *Gradient = new ptImage;
2928
2929 Gradient->Set(this);
2930
2931 ptCimgEdgeDetectionSum(Gradient, Weight, GradientMode);
2932
2933 Gradient->ptCIBlur(Radius, 1);
2934
2935 Gradient->ApplyCurve(Curve, 1);
2936
2937 if (Mode != TOverlayMode::Replace)
2938 Overlay(Gradient->m_Image, 1.0f, NULL, Mode, SwitchLayer);
2939 else
2940 Overlay(Gradient->m_Image, 1.0f, NULL, Mode);
2941
2942 delete Gradient;
2943
2944 return this;
2945 }
2946
2947 ////////////////////////////////////////////////////////////////////////////////
2948 //
2949 // Color Enhance
2950 // http://docs.google.com/View?id=dsgjq79_829f9wv8ncd
2951 //
2952 ////////////////////////////////////////////////////////////////////////////////
2953
ColorEnhance(const float AShadows,const float AHighlights)2954 ptImage* ptImage::ColorEnhance(const float AShadows,
2955 const float AHighlights)
2956 {
2957 assert (m_ColorSpace != ptSpace_Lab);
2958
2959 if (AShadows) {
2960 ptImage *ShadowsLayer = new ptImage;
2961 ShadowsLayer->Set(this);
2962
2963 // Invert and greyscale
2964 #pragma omp parallel for default(shared) schedule(static)
2965 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2966 ShadowsLayer->m_Image[i][0] = ToInvertTable[RGB_2_L(ShadowsLayer->m_Image[i])];
2967 ShadowsLayer->m_Image[i][1] = ShadowsLayer->m_Image[i][2] = ShadowsLayer->m_Image[i][0];
2968 }
2969
2970 ShadowsLayer->Overlay(m_Image, 1.0, NULL, TOverlayMode::ColorDodge, 1 /*Swap */);
2971 Overlay(ShadowsLayer->m_Image, AShadows, NULL, TOverlayMode::ColorBurn);
2972 delete ShadowsLayer;
2973 }
2974 // I trade processing time for memory, so invert and greyscale will be
2975 // recalculated to save another parallel memory instance
2976
2977 if (AHighlights) {
2978 ptImage *HighlightsLayer = new ptImage;
2979 HighlightsLayer->Set(this);
2980
2981 // Invert and greyscale
2982 #pragma omp parallel for default(shared) schedule(static)
2983 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
2984 HighlightsLayer->m_Image[i][0] = ToInvertTable[RGB_2_L(HighlightsLayer->m_Image[i])];
2985 HighlightsLayer->m_Image[i][1] = HighlightsLayer->m_Image[i][2] = HighlightsLayer->m_Image[i][0];
2986 }
2987
2988 HighlightsLayer->Overlay(m_Image, 1.0, NULL, TOverlayMode::ColorBurn, 1 /*Swap */);
2989 Overlay(HighlightsLayer->m_Image, AHighlights, NULL, TOverlayMode::ColorDodge);
2990 delete HighlightsLayer;
2991 }
2992 return this;
2993 }
2994
2995 ////////////////////////////////////////////////////////////////////////////////
2996 //
2997 // LMHLightRecovery
2998 //
2999 ////////////////////////////////////////////////////////////////////////////////
3000
LMHRecovery(const TMaskType MaskType,const float Amount,const float LowerLimit,const float UpperLimit,const float Softness)3001 ptImage* ptImage::LMHRecovery(const TMaskType MaskType,
3002 const float Amount,
3003 const float LowerLimit,
3004 const float UpperLimit,
3005 const float Softness)
3006 {
3007 const float ExposureFactor = pow(2,Amount);
3008 const float InverseExposureFactor = 1/ExposureFactor;
3009
3010 // Precalculated table for the transform of the original.
3011 // The transform is an exposure (>1) or a gamma driven darkening.
3012 // (Table generates less math except for images with < 20K pixels)
3013 uint16_t TransformTable[0x10000];
3014 #pragma omp parallel for
3015 for (uint32_t i=0; i<0x10000; i++) {
3016 if (ExposureFactor<1.0) {
3017 TransformTable[i] = CLIP((int32_t)(powf(i*ptInvWP,InverseExposureFactor)*ptWPf));
3018 } else {
3019 TransformTable[i] = CLIP((int32_t)(i*ExposureFactor+0.5f));
3020 }
3021 }
3022
3023 const double Soft = pow(2,Softness);
3024
3025 // Precalculated table for softening the mask.
3026 double SoftTable[0x100]; // Assuming a 256 table is fine grained enough.
3027 for (int16_t i=0; i<0x100; i++) {
3028 if (Soft>1.0) {
3029 SoftTable[i] = ptBound((float)(pow(i/(float)0xff,Soft)), 0.0f, 1.0f);
3030 } else {
3031 SoftTable[i] = ptBound((float)(i/(float)0xff/Soft),0.0f,1.0f);
3032 }
3033 }
3034
3035 const short NrChannels = (m_ColorSpace == ptSpace_Lab)?1:3;
3036
3037 const float ReciprocalRange = 1.0f/MAX(UpperLimit-LowerLimit,0.001f);
3038 const float ReciprocalLowerLimit = 1.0f/MAX(LowerLimit,0.001f);
3039 const float ReciprocalUpperMargin = 1.0f/MAX(1.0f-UpperLimit,0.001f);
3040 #pragma omp parallel for
3041 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3042 // Init Mask with luminance
3043 float Mask = ((m_ColorSpace == ptSpace_Lab) ?
3044 m_Image[i][0] :
3045 // Remark this classic 30/59/11 should be in fact colour space
3046 // dependent. TODO
3047 RGB_2_L(m_Image[i]))*ptInvWP;
3048 switch(MaskType) {
3049 case TMaskType::Shadows:
3050 // Mask is an inverted luminance mask, normalized 0..1 and
3051 // shifted over the limits such that it clips beyond the limits.
3052 // The Mask varies from 1 at LowerLimit to 0 at UpperLimit
3053 // Meaning that deep shadows will be pulled up a lot and
3054 // approximating the upperlimit we will take more of the original
3055 // image.
3056 Mask = 1.0f-LIM((Mask-LowerLimit)*ReciprocalRange,0.0f,1.0f);
3057 break;
3058 case TMaskType::Midtones:
3059 // Not fully understood but generates a useful and nice
3060 // midtone luminance mask.
3061 Mask = 1.0f -
3062 LIM((LowerLimit-Mask)*ReciprocalLowerLimit,0.0f,0.1f) -
3063 LIM((Mask-UpperLimit)*ReciprocalUpperMargin,0.0f,1.0f);
3064 Mask = LIM(Mask,0.0f,1.0f);
3065 break;
3066 case TMaskType::Highlights:
3067 // Mask is a luminance mask, normalized 0..1 and
3068 // shifted over the limits such that it clips beyond the limits.
3069 // The Mask varies from 0 at LowerLimit to 1 at UpperLimit
3070 // Meaning that as from the LowerLimit on , we will take more and
3071 // more of the darkened image.
3072 Mask = LIM((Mask-LowerLimit)*ReciprocalRange,0.0f,1.0f);
3073 break;
3074 case TMaskType::All:
3075 Mask = 1.0f;
3076 break;
3077
3078 default:
3079 GInfo->Raise(QString("Unknown mask type: ") + QString::number(static_cast<int>(MaskType)), AT);
3080 }
3081
3082 // Softening the mask
3083 Mask = SoftTable[(uint8_t)(Mask*0xff+0.5)];
3084
3085 // Blend transformed and original according to mask.
3086 for (short Ch=0; Ch<NrChannels; Ch++) {
3087 uint16_t PixelValue = m_Image[i][Ch];
3088 m_Image[i][Ch] = CLIP((int32_t)
3089 (TransformTable[PixelValue]*Mask + PixelValue*(1.0f-Mask)));
3090 // Uncomment me to 'see' the mask.
3091 // m_Image[i][Ch] = Mask*WP;
3092 }
3093 }
3094 return this;
3095 }
3096
3097
3098 ////////////////////////////////////////////////////////////////////////////////
3099 //
3100 // Highpass
3101 //
3102 ////////////////////////////////////////////////////////////////////////////////
3103
Highpass(const double Radius,const double Amount,const double HaloControl,const double Denoise)3104 ptImage* ptImage::Highpass(const double Radius,
3105 const double Amount,
3106 const double HaloControl,
3107 const double Denoise)
3108 {
3109 double LowerLimit = 0.1;
3110 double UpperLimit = 1 - LowerLimit;
3111 double Softness = 0;
3112
3113 const double WPH = 0x7fff; // WPH=WP/2
3114 const short NrChannels = (m_ColorSpace == ptSpace_Lab)?1:3;
3115 const TChannelMask ChannelMask = (m_ColorSpace == ptSpace_Lab) ? ChMask_L : ChMask_RGB;
3116
3117 ptImage *HighpassLayer = new ptImage;
3118 HighpassLayer->Set(this);
3119 HighpassLayer->ptCIBlur(Radius, ChannelMask);
3120
3121 // also calculates the curve
3122 auto AmpCurve = new ptCurve(createAmpAnchors(Amount, HaloControl));
3123
3124 #pragma omp parallel for default(shared)
3125 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3126 for (short Ch=0; Ch<NrChannels; Ch++) {
3127 HighpassLayer->m_Image[i][Ch] = CLIP((int32_t) ((WPH-(int32_t)HighpassLayer->m_Image[i][Ch])+m_Image[i][Ch]));
3128 }
3129 }
3130
3131 HighpassLayer->ApplyCurve(AmpCurve,ChannelMask);
3132 delete AmpCurve;
3133
3134 if (Denoise) {
3135 HighpassLayer->fastBilateralChannel(4.0, Denoise/3.0, 1, ChMask_L);
3136 }
3137
3138 float (*Mask);
3139 Mask = (m_ColorSpace == ptSpace_Lab)?
3140 GetMask(TMaskType::Midtones, LowerLimit, UpperLimit, Softness,1,0,0):
3141 GetMask(TMaskType::Midtones, LowerLimit, UpperLimit, Softness);
3142
3143 Overlay(HighpassLayer->m_Image,0.5,Mask);
3144 delete HighpassLayer;
3145 FREE(Mask);
3146
3147 return this;
3148 }
3149
3150 ////////////////////////////////////////////////////////////////////////////////
3151 //
3152 // Gradient Sharpen
3153 //
3154 ////////////////////////////////////////////////////////////////////////////////
3155
3156 // To the extent possible under law, Manuel Llorens <manuelllorens@gmail.com>
3157 // has waived all copyright and related or neighboring rights to this work.
3158 // This code is licensed under CC0 v1.0, see license information at
3159 // http://creativecommons.org/publicdomain/zero/1.0/
3160
GradientSharpen(const short Passes,const double Strength)3161 ptImage* ptImage::GradientSharpen(const short Passes,
3162 const double Strength) {
3163
3164 assert (m_ColorSpace == ptSpace_Lab);
3165
3166 int32_t offset,c,i,j,p,width2;
3167 float lumH,lumV,lumD1,lumD2,v,contrast,s;
3168 float difL,difR,difT,difB,difLT,difRB,difLB,difRT,wH,wV,wD1,wD2,chmax[3];
3169 float f1,f2,f3,f4;
3170
3171 uint16_t width = m_Width;
3172 uint16_t height = m_Height;
3173
3174 width2=2*m_Width;
3175
3176 float (*L) = (float (*)) CALLOC(m_Width*m_Height,sizeof(*L));
3177 ptMemoryError(L,__FILE__,__LINE__);
3178
3179 chmax[0]=0.08;
3180 chmax[1]=3.0;
3181 chmax[2]=3.0;
3182 c = 0;
3183
3184 #pragma omp parallel for private(offset) schedule(static)
3185 for(offset=0;offset<m_Width*m_Height;offset++)
3186 L[offset]=ToFloatTable[m_Image[offset][c]];
3187
3188 //for(c=0;c<=channels;c++)
3189 for(p=0;p<Passes;p++){
3190 #pragma omp parallel for private(j,i,offset,wH,wV,wD1,wD2,s,lumH,lumV,lumD1,lumD2,v,contrast,f1,f2,f3,f4,difT,difB,difL,difR,difLT,difLB,difRT,difRB) schedule(static)
3191 for(j=2;j<height-2;j++)
3192 for(i=2,offset=j*width+i;i<width-2;i++,offset++){
3193 // weight functions
3194 wH=fabsf(L[offset+1]-L[offset-1]);
3195 wV=fabsf(L[offset+width]-L[offset-width]);
3196
3197 s=1.0+fabs(wH-wV)/2.0;
3198 wD1=fabsf(L[offset+width+1]-L[offset-width-1])/s;
3199 wD2=fabsf(L[offset+width-1]-L[offset-width+1])/s;
3200 s=wD1;
3201 wD1/=wD2;
3202 wD2/=wD1;
3203
3204 // initial values
3205 lumH=lumV=lumD1=lumD2=v=ToFloatTable[m_Image[offset][c]];
3206
3207 // contrast detection
3208 contrast=sqrtf(fabsf(L[offset+1]-L[offset-1])*fabsf(L[offset+1]-L[offset-1])+fabsf(L[offset+width]-L[offset-width])*fabsf(L[offset+width]-L[offset-width]))/chmax[c];
3209 if(contrast>1.0) contrast=1.0;
3210
3211 // new possible values
3212 if(((L[offset]<L[offset-1])&&(L[offset]>L[offset+1])) ||
3213 ((L[offset]>L[offset-1])&&(L[offset]<L[offset+1]))){
3214 f1=fabsf(L[offset-2]-L[offset-1]);
3215 f2=fabsf(L[offset-1]-L[offset]);
3216 f3=fabsf(L[offset-1]-L[offset-width])*fabsf(L[offset-1]-L[offset+width]);
3217 f4=sqrtf(fabsf(L[offset-1]-L[offset-width2])*fabsf(L[offset-1]-L[offset+width2]));
3218 difL=f1*f2*f2*f3*f3*f4;
3219 f1=fabsf(L[offset+2]-L[offset+1]);
3220 f2=fabsf(L[offset+1]-L[offset]);
3221 f3=fabsf(L[offset+1]-L[offset-width])*fabsf(L[offset+1]-L[offset+width]);
3222 f4=sqrtf(fabs(L[offset+1]-L[offset-width2])*fabsf(L[offset+1]-L[offset+width2]));
3223 difR=f1*f2*f2*f3*f3*f4;
3224 if((difR!=0)&&(difL!=0)){
3225 lumH=(L[offset-1]*difR+L[offset+1]*difL)/(difL+difR);
3226 lumH=v*(1-contrast)+lumH*contrast;
3227 }
3228 }
3229
3230 if(((L[offset]<L[offset-width])&&(L[offset]>L[offset+width])) ||
3231 ((L[offset]>L[offset-width])&&(L[offset]<L[offset+width]))){
3232 f1=fabsf(L[offset-width2]-L[offset-width]);
3233 f2=fabsf(L[offset-width]-L[offset]);
3234 f3=fabsf(L[offset-width]-L[offset-1])*fabsf(L[offset-width]-L[offset+1]);
3235 f4=sqrtf(fabsf(L[offset-width]-L[offset-2])*fabsf(L[offset-width]-L[offset+2]));
3236 difT=f1*f2*f2*f3*f3*f4;
3237 f1=fabsf(L[offset+width2]-L[offset+width]);
3238 f2=fabsf(L[offset+width]-L[offset]);
3239 f3=fabsf(L[offset+width]-L[offset-1])*fabsf(L[offset+width]-L[offset+1]);
3240 f4=sqrtf(fabsf(L[offset+width]-L[offset-2])*fabsf(L[offset+width]-L[offset+2]));
3241 difB=f1*f2*f2*f3*f3*f4;
3242 if((difB!=0)&&(difT!=0)){
3243 lumV=(L[offset-width]*difB+L[offset+width]*difT)/(difT+difB);
3244 lumV=v*(1-contrast)+lumV*contrast;
3245 }
3246 }
3247
3248 if(((L[offset]<L[offset-1-width])&&(L[offset]>L[offset+1+width])) ||
3249 ((L[offset]>L[offset-1-width])&&(L[offset]<L[offset+1+width]))){
3250 f1=fabsf(L[offset-2-width2]-L[offset-1-width]);
3251 f2=fabsf(L[offset-1-width]-L[offset]);
3252 f3=fabsf(L[offset-1-width]-L[offset-width+1])*fabsf(L[offset-1-width]-L[offset+width-1]);
3253 f4=sqrtf(fabsf(L[offset-1-width]-L[offset-width2+2])*fabsf(L[offset-1-width]-L[offset+width2-2]));
3254 difLT=f1*f2*f2*f3*f3*f4;
3255 f1=fabsf(L[offset+2+width2]-L[offset+1+width]);
3256 f2=fabsf(L[offset+1+width]-L[offset]);
3257 f3=fabsf(L[offset+1+width]-L[offset-width+1])*fabsf(L[offset+1+width]-L[offset+width-1]);
3258 f4=sqrtf(fabsf(L[offset+1+width]-L[offset-width2+2])*fabsf(L[offset+1+width]-L[offset+width2-2]));
3259 difRB=f1*f2*f2*f3*f3*f4;
3260 if((difLT!=0)&&(difRB!=0)){
3261 lumD1=(L[offset-1-width]*difRB+L[offset+1+width]*difLT)/(difLT+difRB);
3262 lumD1=v*(1-contrast)+lumD1*contrast;
3263 }
3264 }
3265
3266 if(((L[offset]<L[offset+1-width])&&(L[offset]>L[offset-1+width])) ||
3267 ((L[offset]>L[offset+1-width])&&(L[offset]<L[offset-1+width]))){
3268 f1=fabsf(L[offset-2+width2]-L[offset-1+width]);
3269 f2=fabsf(L[offset-1+width]-L[offset]);
3270 f3=fabsf(L[offset-1+width]-L[offset-width-1])*fabsf(L[offset-1+width]-L[offset+width+1]);
3271 f4=sqrtf(fabsf(L[offset-1+width]-L[offset-width2-2])*fabsf(L[offset-1+width]-L[offset+width2+2]));
3272 difLB=f1*f2*f2*f3*f3*f4;
3273 f1=fabsf(L[offset+2-width2]-L[offset+1-width]);
3274 f2=fabsf(L[offset+1-width]-L[offset])*fabsf(L[offset+1-width]-L[offset]);
3275 f3=fabsf(L[offset+1-width]-L[offset+width+1])*fabsf(L[offset+1-width]-L[offset-width-1]);
3276 f4=sqrtf(fabsf(L[offset+1-width]-L[offset+width2+2])*fabsf(L[offset+1-width]-L[offset-width2-2]));
3277 difRT=f1*f2*f2*f3*f3*f4;
3278 if((difLB!=0)&&(difRT!=0)){
3279 lumD2=(L[offset+1-width]*difLB+L[offset-1+width]*difRT)/(difLB+difRT);
3280 lumD2=v*(1-contrast)+lumD2*contrast;
3281 }
3282 }
3283
3284 s=Strength;
3285
3286 // avoid sharpening diagonals too much
3287 if(((fabsf(wH/wV)<0.45)&&(fabsf(wH/wV)>0.05))||((fabsf(wV/wH)<0.45)&&(fabsf(wV/wH)>0.05)))
3288 s=Strength/3.0;
3289
3290 // final mix
3291 if((wH!=0)&&(wV!=0)&&(wD1!=0)&&(wD2!=0))
3292 m_Image[offset][c]= CLIP((int32_t) ((v*(1-s)+(lumH*wH+lumV*wV+lumD1*wD1+lumD2*wD2)/(wH+wV+wD1+wD2)*s)*0xffff));
3293 }
3294 }
3295
3296 FREE(L);
3297 return this;
3298 }
3299
3300 // To the extent possible under law, Manuel Llorens <manuelllorens@gmail.com>
3301 // has waived all copyright and related or neighboring rights to this work.
3302 // This code is licensed under CC0 v1.0, see license information at
3303 // http://creativecommons.org/publicdomain/zero/1.0/
3304
MLMicroContrast(const double Strength,const double Scaling,const double Weight,const ptCurve * Curve)3305 ptImage* ptImage::MLMicroContrast(const double Strength,
3306 const double Scaling,
3307 const double Weight,
3308 const ptCurve *Curve)
3309 {
3310 assert (m_ColorSpace == ptSpace_Lab);
3311
3312 int32_t offset,offset2,c,i,j,col,row,n;
3313 float v,s,contrast,temp;
3314 float CompWeight = 1.0f - Weight;
3315 const float WPH = 0x8080;
3316
3317 float ValueA = 0.0;
3318 float ValueB = 0.0;
3319
3320 uint16_t width = m_Width;
3321 uint16_t height = m_Height;
3322
3323 float (*L) = (float (*)) CALLOC(m_Width*m_Height,sizeof(*L));
3324 ptMemoryError(L,__FILE__,__LINE__);
3325
3326 int signs[9];
3327
3328 float chmax = 8.0f/Scaling;
3329
3330 c=0;
3331
3332 #pragma omp parallel for private(offset) schedule(static)
3333 for(offset=0;offset<m_Width*m_Height;offset++)
3334 L[offset]=ToFloatTable[m_Image[offset][c]];
3335
3336 #pragma omp parallel for private(j,i,offset,s,signs,v,n,row,col,offset2,contrast,temp, ValueA, ValueB) schedule(static)
3337 for(j=1;j<height-1;j++)
3338 for(i=1,offset=j*width+i;i<width-1;i++,offset++){
3339 if (Curve == NULL) s=Strength;
3340 else {
3341 // set s according to the curve
3342 if (Curve->mask() == ptCurve::ChromaMask) {
3343 ValueA = (float)m_Image[offset][1]-WPH;
3344 ValueB = (float)m_Image[offset][2]-WPH;
3345 float Hue = 0;
3346 if (ValueA == 0.0 && ValueB == 0.0) {
3347 Hue = 0; // value for grey pixel
3348 } else {
3349 Hue = atan2f(ValueB,ValueA);
3350 }
3351 while (Hue < 0) Hue += 2.*ptPI;
3352
3353 float Col = powf(ValueA * ValueA + ValueB * ValueB, 0.125);
3354 Col /= 0x7; // normalizing to 0..2
3355
3356 float Factor = Curve->Curve[CLIP((int32_t)(Hue/ptPI*WPH))]/(float)0x3333 - 1.0;
3357 s = Strength * Factor * Col;
3358
3359 } else { //by luma
3360 float Factor = Curve->Curve[m_Image[offset][0]]/(float)0x3333 - 1.0;
3361 s = Strength * Factor;
3362 }
3363 }
3364 v=L[offset];
3365
3366 n=0;
3367 for(row=j-1;row<=j+1;row++) {
3368 for(col=i-1,offset2=row*width+col;col<=i+1;col++,offset2++){
3369 signs[n]=0;
3370 if(v<L[offset2]) signs[n]=-1;
3371 if(v>L[offset2]) signs[n]=1;
3372 n++;
3373 }
3374 }
3375
3376 contrast=sqrtf(fabsf(L[offset+1]-L[offset-1])*fabsf(L[offset+1]-L[offset-1])+fabsf(L[offset+width]-L[offset-width])*fabsf(L[offset+width]-L[offset-width]))/chmax;
3377 if(contrast>1.0) contrast=1.0;
3378 temp = ToFloatTable[m_Image[offset][c]];
3379 temp +=(v-L[offset-width-1])*sqrtf(2)*s;
3380 temp +=(v-L[offset-width])*s;
3381 temp +=(v-L[offset-width+1])*sqrtf(2)*s;
3382
3383 temp +=(v-L[offset-1])*s;
3384 temp +=(v-L[offset+1])*s;
3385
3386 temp +=(v-L[offset+width-1])*sqrtf(2)*s;
3387 temp +=(v-L[offset+width])*s;
3388 temp +=(v-L[offset+width+1])*sqrtf(2)*s;
3389
3390 temp = MAX(0.0f,temp);
3391
3392 // Reduce halo looking artifacs
3393 v=temp;
3394 n=0;
3395 for(row=j-1;row<=j+1;row++)
3396 for(col=i-1,offset2=row*width+col;col<=i+1;col++,offset2++){
3397 if(((v<L[offset2])&&(signs[n]>0))||((v>L[offset2])&&(signs[n]<0)))
3398 temp=v*Weight+L[offset2]*CompWeight;
3399 n++;
3400 }
3401 m_Image[offset][c]=CLIP((int32_t) ((temp*(1-contrast)+L[offset]*contrast)*0xffff));
3402 }
3403
3404 FREE(L);
3405 return this;
3406 }
3407
3408 ////////////////////////////////////////////////////////////////////////////////
3409 //
3410 // Hotpixel
3411 //
3412 ////////////////////////////////////////////////////////////////////////////////
3413
HotpixelReduction(const double Threshold)3414 ptImage* ptImage::HotpixelReduction(const double Threshold) {
3415 uint16_t Thres = (int32_t) ((1.0-Threshold)*0x2fff);
3416 uint16_t Width = m_Width;
3417 uint16_t Height = m_Height;
3418
3419 #pragma omp parallel for schedule(static) default(shared)
3420 for (uint16_t Row=0; Row<Height; Row++) {
3421 for (uint16_t Col=0; Col<Width; Col++) {
3422 uint16_t TempValue = 0;
3423 for (int c=0; c<3; c++) {
3424 // bright pixels
3425 if (Row > 1) {
3426 TempValue = MAX(m_Image[(Row-1)*Width+Col][c],TempValue);
3427 if (Col > 1) TempValue = MAX(m_Image[(Row-1)*Width+Col-1][c],TempValue);
3428 if (Col < Width-1) TempValue = MAX(m_Image[(Row-1)*Width+Col+1][c],TempValue);
3429 }
3430 if (Row < Height-1) {
3431 TempValue = MAX(m_Image[(Row+1)*Width+Col][c],TempValue);
3432 if (Col > 1) TempValue = MAX(m_Image[(Row+1)*Width+Col-1][c],TempValue);
3433 if (Col < Width-1) TempValue = MAX(m_Image[(Row+1)*Width+Col+1][c],TempValue);
3434 }
3435 if (Col > 1) TempValue = MAX(m_Image[Row*Width+Col-1][c],TempValue);
3436 if (Col < Width-1) TempValue = MAX(m_Image[Row*Width+Col+1][c],TempValue);
3437 if (TempValue+Thres<m_Image[Row*Width+Col][c])
3438 m_Image[Row*Width+Col][c] = TempValue;
3439
3440 // dark pixels
3441 TempValue = 0xffff;
3442 if (Row > 1) {
3443 TempValue = MIN(m_Image[(Row-1)*Width+Col][c],TempValue);
3444 if (Col > 1) TempValue = MIN(m_Image[(Row-1)*Width+Col-1][c],TempValue);
3445 if (Col < Width-1) TempValue = MIN(m_Image[(Row-1)*Width+Col+1][c],TempValue);
3446 }
3447 if (Row < Height-1) {
3448 TempValue = MIN(m_Image[(Row+1)*Width+Col][c],TempValue);
3449 if (Col > 1) TempValue = MIN(m_Image[(Row+1)*Width+Col-1][c],TempValue);
3450 if (Col < Width-1) TempValue = MIN(m_Image[(Row+1)*Width+Col+1][c],TempValue);
3451 }
3452 if (Col > 1) TempValue = MIN(m_Image[Row*Width+Col-1][c],TempValue);
3453 if (Col < Width-1) TempValue = MIN(m_Image[Row*Width+Col+1][c],TempValue);
3454 if (TempValue-Thres>m_Image[Row*Width+Col][c])
3455 m_Image[Row*Width+Col][c] = TempValue;
3456 }
3457 }
3458 }
3459 return this;
3460 }
3461
3462 ////////////////////////////////////////////////////////////////////////////////
3463 //
3464 // Refined Shadows and Highlights
3465 //
3466 ////////////////////////////////////////////////////////////////////////////////
3467
ShadowsHighlights(const ptCurve * Curve,const double Radius,const double AmountCoarse,const double AmountFine)3468 ptImage* ptImage::ShadowsHighlights(const ptCurve *Curve,
3469 const double Radius,
3470 const double AmountCoarse,
3471 const double AmountFine) {
3472
3473 assert (m_ColorSpace == ptSpace_Lab);
3474
3475 uint16_t (*CoarseLayer) = (uint16_t (*)) CALLOC(m_Width*m_Height,sizeof(*CoarseLayer));
3476 ptMemoryError(CoarseLayer,__FILE__,__LINE__);
3477 uint16_t (*FineLayer) = (uint16_t (*)) CALLOC(m_Width*m_Height,sizeof(*FineLayer));
3478 ptMemoryError(FineLayer,__FILE__,__LINE__);
3479
3480 #pragma omp parallel for default(shared)
3481 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3482 CoarseLayer[i] = m_Image[i][0];
3483 }
3484
3485 this->fastBilateralChannel(Radius, 0.14, 2, ChMask_L);
3486
3487 #pragma omp parallel for default(shared)
3488 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3489 FineLayer[i] = CLIP((int32_t) ((ptWPH - (int32_t)m_Image[i][0]) + CoarseLayer[i]));
3490 CoarseLayer[i] = m_Image[i][0];
3491 }
3492
3493 this->fastBilateralChannel(4*Radius, 0.14, 2, ChMask_L);
3494
3495 #pragma omp parallel for default(shared)
3496 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3497 CoarseLayer[i] = CLIP((int32_t) ((ptWPH - (int32_t)m_Image[i][0]) + CoarseLayer[i]));
3498 }
3499
3500 //Curve on residual
3501 if (Curve != NULL) {
3502 ApplyCurve(Curve,1);
3503 }
3504
3505 //Sigmoidal Contrast
3506 float Threshold = 0.5;
3507 float Contrast = AmountCoarse;
3508 uint16_t ContrastTable[0x10000];
3509
3510 if (AmountCoarse != 0) {
3511 SigmoidalTable(ContrastTable, Contrast, Threshold);
3512 #pragma omp parallel for default(shared)
3513 for (uint32_t i=0; i < (uint32_t)m_Height*m_Width; i++) {
3514 CoarseLayer[i] = ContrastTable[ CoarseLayer[i] ];
3515 }
3516 }
3517
3518 Contrast = AmountFine;
3519
3520 if (AmountFine != 0) {
3521 SigmoidalTable(ContrastTable, Contrast, Threshold);
3522 #pragma omp parallel for default(shared)
3523 for (uint32_t i=0; i < (uint32_t)m_Height*m_Width; i++) {
3524 FineLayer[i] = ContrastTable[ FineLayer[i] ];
3525 }
3526 }
3527
3528 #pragma omp parallel for default(shared)
3529 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3530 m_Image[i][0] = CLIP((int32_t) (((int32_t)CoarseLayer[i] - ptWPH) + m_Image[i][0]));
3531 m_Image[i][0] = CLIP((int32_t) (((int32_t)FineLayer[i] - ptWPH) + m_Image[i][0]));
3532 }
3533
3534 FREE(CoarseLayer);
3535 FREE(FineLayer);
3536
3537 return this;
3538 }
3539
3540 ////////////////////////////////////////////////////////////////////////////////
3541 //
3542 // Microcontrast
3543 //
3544 ////////////////////////////////////////////////////////////////////////////////
3545
Microcontrast(const double Radius,const double Amount,const double Opacity,const double HaloControl,const TMaskType MaskType,const double LowerLimit,const double UpperLimit,const double Softness)3546 ptImage* ptImage::Microcontrast(const double Radius,
3547 const double Amount,
3548 const double Opacity,
3549 const double HaloControl,
3550 const TMaskType MaskType,
3551 const double LowerLimit,
3552 const double UpperLimit,
3553 const double Softness) {
3554
3555 const double WPH = 0x7fff; // WPH=WP/2
3556 const short NrChannels = (m_ColorSpace == ptSpace_Lab)?1:3;
3557 const short ChannelMask = (m_ColorSpace == ptSpace_Lab)?1:7;
3558
3559 ptImage *MicrocontrastLayer = new ptImage;
3560 MicrocontrastLayer->Set(this);
3561 MicrocontrastLayer->ptCIBlur(Radius, ChannelMask);
3562
3563 auto AmpCurve = new ptCurve(createAmpAnchors(Amount, HaloControl));
3564
3565 #pragma omp parallel for default(shared)
3566 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3567 for (short Ch=0; Ch<NrChannels; Ch++) {
3568 MicrocontrastLayer->m_Image[i][Ch] = CLIP((int32_t) ((WPH-(int32_t)MicrocontrastLayer->m_Image[i][Ch])+m_Image[i][Ch]));
3569 }
3570 }
3571
3572 MicrocontrastLayer->ApplyCurve(AmpCurve,ChannelMask);
3573 delete AmpCurve;
3574
3575 float (*Mask);
3576 Mask = (m_ColorSpace == ptSpace_Lab)?
3577 GetMask( MaskType, LowerLimit, UpperLimit, Softness,1,0,0):
3578 GetMask( MaskType, LowerLimit, UpperLimit, Softness);
3579
3580 Overlay(MicrocontrastLayer->m_Image,Opacity,Mask);
3581 delete MicrocontrastLayer;
3582
3583 FREE(Mask);
3584
3585 return this;
3586 }
3587
3588 ////////////////////////////////////////////////////////////////////////////////
3589 //
3590 // Colorcontrast
3591 //
3592 ////////////////////////////////////////////////////////////////////////////////
3593
Colorcontrast(const double Radius,const double Amount,const double Opacity,const double HaloControl)3594 ptImage* ptImage::Colorcontrast(const double Radius,
3595 const double Amount,
3596 const double Opacity,
3597 const double HaloControl)
3598 {
3599 const double WP = 0xffff;
3600 const double WPH = 0x7fff; // WPH=WP/2
3601 const short NrChannels = 3;
3602 const short ChannelMask = 6;
3603
3604 ptImage *MicrocontrastLayer = new ptImage;
3605 MicrocontrastLayer->Set(this);
3606 MicrocontrastLayer->ptCIBlur(Radius, ChannelMask);
3607
3608 auto AmpCurve = new ptCurve(createAmpAnchors(Amount, HaloControl));
3609
3610 #pragma omp parallel for default(shared)
3611 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3612 for (short Ch=1; Ch<NrChannels; Ch++) {
3613 MicrocontrastLayer->m_Image[i][Ch] = CLIP((int32_t) ((WPH-(int32_t)MicrocontrastLayer->m_Image[i][Ch])+m_Image[i][Ch]));
3614 }
3615 }
3616
3617 MicrocontrastLayer->ApplyCurve(AmpCurve,ChannelMask);
3618 delete AmpCurve;
3619
3620 float Multiply = 0;
3621 float Screen = 0;
3622 float Overlay = 0;
3623 float Source = 0;
3624 float Blend = 0;
3625
3626 for (short Ch=1; Ch<3; Ch++) {
3627 // Is it a channel we are supposed to handle ?
3628 if (! (ChannelMask & (1<<Ch))) continue;
3629 #pragma omp parallel for default(shared) private(Source, Blend, Multiply, Screen, Overlay)
3630 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3631 Source = m_Image[i][Ch];
3632 Blend = MicrocontrastLayer->m_Image[i][Ch];
3633 Multiply = CLIP((int32_t)(Source*Blend/WP));
3634 Screen = CLIP((int32_t)(WP-(WP-Source)*(WP-Blend)/WP));
3635 Overlay = CLIP((int32_t)((((WP-Source)*Multiply+Source*Screen)/WP)));
3636 m_Image[i][Ch] = CLIP((int32_t) (-WP*MIN(Opacity, 0.0)+Overlay*Opacity+Source*(1-fabs(Opacity))));
3637 }
3638 }
3639
3640 delete MicrocontrastLayer;
3641
3642 return this;
3643 }
3644
3645 ////////////////////////////////////////////////////////////////////////////////
3646 //
3647 // Bilateral Denoise
3648 //
3649 ////////////////////////////////////////////////////////////////////////////////
3650
BilateralDenoise(const double Threshold,const double Softness,const double Opacity,const double UseMask)3651 ptImage* ptImage::BilateralDenoise(const double Threshold,
3652 const double Softness,
3653 const double Opacity,
3654 const double UseMask /* = 0*/) {
3655
3656 ptImage *DenoiseLayer = new ptImage;
3657 DenoiseLayer->Set(this);
3658 const short NrChannels = (m_ColorSpace == ptSpace_Lab)?1:3;
3659 const TChannelMask ChannelMask = (m_ColorSpace == ptSpace_Lab) ? ChMask_L : ChMask_RGB;
3660 const double WPH = 0x7fff;
3661
3662 DenoiseLayer->fastBilateralChannel(Threshold, Softness, 2, ChannelMask);
3663
3664 if (UseMask){
3665
3666 double m = 10.0;
3667 double t = (1.0 - m)*WPH;
3668
3669 uint16_t Table[0x10000];
3670 #pragma omp parallel for schedule(static)
3671 for (uint32_t i=0; i<0x10000; i++) {
3672 Table[i] = CLIP((int32_t)(m*i+t));
3673 }
3674
3675 ptImage *MaskLayer = new ptImage;
3676 MaskLayer->m_Colors = 3;
3677 MaskLayer->Set(DenoiseLayer);
3678
3679 #pragma omp parallel for default(shared) schedule(static)
3680 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3681 for (short Ch=0; Ch<NrChannels; Ch++) {
3682 MaskLayer->m_Image[i][Ch] = CLIP((int32_t) ((WPH-(int32_t)DenoiseLayer->m_Image[i][Ch])+m_Image[i][Ch]));
3683 MaskLayer->m_Image[i][Ch] = Table[MaskLayer->m_Image[i][Ch]];
3684 }
3685 }
3686
3687 if (ChannelMask == 7) {
3688 #pragma omp parallel for default(shared) schedule(static)
3689 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++)
3690 MaskLayer->m_Image[i][0] = RGB_2_L(MaskLayer->m_Image[i]);
3691 }
3692
3693 ptCurve* Curve = new ptCurve({TAnchor(0.0, 1.0),
3694 TAnchor(0.4, 0.3),
3695 TAnchor(0.5, 0.0),
3696 TAnchor(0.6, 0.3),
3697 TAnchor(1.0, 1.0)});
3698
3699 MaskLayer->ApplyCurve(Curve,1);
3700 MaskLayer->ptCIBlur(UseMask, 1);
3701
3702 Curve->setFromAnchors({TAnchor(0.0, 0.0),
3703 TAnchor(0.6, 0.4),
3704 TAnchor(1.0, 0.8)});
3705
3706 MaskLayer->ApplyCurve(Curve,1);
3707 delete Curve;
3708
3709 float (*Mask);
3710 Mask = MaskLayer->GetMask(TMaskType::Shadows, 0.0, 1.0, 0.0, 1,0,0);
3711 Overlay(DenoiseLayer->m_Image,Opacity,Mask,TOverlayMode::Normal);
3712 FREE(Mask);
3713 delete MaskLayer;
3714 } else {
3715 Overlay(DenoiseLayer->m_Image,Opacity,NULL,TOverlayMode::Normal);
3716 }
3717 delete DenoiseLayer;
3718 return this;
3719 }
3720
3721 ////////////////////////////////////////////////////////////////////////////////
3722 //
3723 // Denoise curve
3724 //
3725 ////////////////////////////////////////////////////////////////////////////////
3726
3727
ApplyDenoiseCurve(const double Threshold,const double Softness,const ptCurve * MaskCurve)3728 ptImage* ptImage::ApplyDenoiseCurve(const double Threshold,
3729 const double Softness,
3730 const ptCurve *MaskCurve) {
3731
3732 assert (m_ColorSpace == ptSpace_Lab);
3733
3734 ptImage *DenoiseLayer = new ptImage;
3735 DenoiseLayer->Set(this);
3736 const short NrChannels = (m_ColorSpace == ptSpace_Lab)?1:3;
3737 const TChannelMask ChannelMask = (m_ColorSpace == ptSpace_Lab) ? ChMask_L : ChMask_RGB;
3738 float WPH = 0x7fff;
3739
3740 DenoiseLayer->fastBilateralChannel(Threshold, Softness, 2, ChannelMask);
3741
3742 double UseMask = 50.0;
3743
3744 float m = 10.0;
3745 float t = (1.0 - m)*WPH;
3746
3747 uint16_t Table[0x10000];
3748 #pragma omp parallel for schedule(static)
3749 for (uint32_t i=0; i<0x10000; i++) {
3750 Table[i] = CLIP((int32_t)(m*i+t));
3751 }
3752
3753 ptImage *MaskLayer = new ptImage;
3754 MaskLayer->m_Colors = 3;
3755 MaskLayer->Set(DenoiseLayer);
3756
3757 uint16_t (*Temp) = (uint16_t(*)) CALLOC(m_Width*m_Height,sizeof(*Temp));
3758 ptMemoryError(Temp,__FILE__,__LINE__);
3759
3760 #pragma omp parallel for default(shared) schedule(static)
3761 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3762 for (short Ch=0; Ch<NrChannels; Ch++) {
3763 MaskLayer->m_Image[i][Ch] = CLIP((int32_t) ((WPH-(int32_t)DenoiseLayer->m_Image[i][Ch])+m_Image[i][Ch]));
3764 MaskLayer->m_Image[i][Ch] = Table[MaskLayer->m_Image[i][Ch]];
3765 }
3766 Temp[i] = m_Image[i][0];
3767 }
3768
3769 ptCurve* Curve = new ptCurve({TAnchor(0.0, 1.0),
3770 TAnchor(0.4, 0.3),
3771 TAnchor(0.5, 0.0),
3772 TAnchor(0.6, 0.3),
3773 TAnchor(1.0, 1.0)}); // also triggers curve calc
3774 MaskLayer->ApplyCurve(Curve, ChMask_L);
3775
3776 MaskLayer->ptCIBlur(UseMask, ChMask_L);
3777
3778 Curve->setFromAnchors({TAnchor(0.0, 0.0),
3779 TAnchor(0.6, 0.4),
3780 TAnchor(1.0, 0.8)});
3781 Curve->calcCurve();
3782 MaskLayer->ApplyCurve(Curve, ChMask_L);
3783 delete Curve;
3784
3785 float (*Mask);
3786 Mask = MaskLayer->GetMask(TMaskType::Shadows, 0.0, 1.0, 0.0, 1,0,0);
3787 delete MaskLayer;
3788 Overlay(DenoiseLayer->m_Image,1.0f,Mask,TOverlayMode::Normal);
3789 FREE(Mask);
3790 delete DenoiseLayer;
3791
3792 // at this point Temp contains the unaltered L channel and m_Image
3793 // contains the denoised layer.
3794
3795 // neutral value for a* and b* channel
3796 WPH = 0x8080;
3797
3798 float ValueA = 0.0;
3799 float ValueB = 0.0;
3800 float Hue = 0.0;
3801
3802 if (MaskCurve->mask() == ptCurve::ChromaMask) {
3803 #pragma omp parallel for schedule(static) private(ValueA, ValueB, Hue)
3804 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
3805 // Factor by hue
3806 ValueA = (float)m_Image[i][1]-WPH;
3807 ValueB = (float)m_Image[i][2]-WPH;
3808
3809 if (ValueA == 0.0 && ValueB == 0.0) {
3810 Hue = 0; // value for grey pixel
3811 } else {
3812 Hue = atan2f(ValueB,ValueA);
3813 }
3814 while (Hue < 0) Hue += 2.*ptPI;
3815
3816 float Factor = MaskCurve->Curve[CLIP((int32_t)(Hue/ptPI*WPH))]-(float)0x7fff;
3817 Factor /= (float)0x7fff;
3818
3819 m_Image[i][0]=CLIP((int32_t)(Temp[i]*(1.0f - Factor)+m_Image[i][0]*Factor ));
3820 }
3821
3822 } else { // by luma
3823 #pragma omp parallel for schedule(static) private(ValueA, ValueB)
3824 for(uint32_t i = 0; i < (uint32_t) m_Width*m_Height; i++) {
3825 // Factor by luminance
3826 float Factor = MaskCurve->Curve[Temp[i]]-(float)0x7fff;
3827 Factor /= (float)0x7fff;
3828
3829 m_Image[i][0]=CLIP((int32_t)(Temp[i]*(1.0f - Factor)+m_Image[i][0]*Factor ));
3830 }
3831 }
3832
3833 FREE(Temp);
3834
3835 return this;
3836 }
3837
3838
3839 ////////////////////////////////////////////////////////////////////////////////
3840 //
3841 // Texturecontrast
3842 //
3843 ////////////////////////////////////////////////////////////////////////////////
3844
TextureContrast(const double Threshold,const double Softness,const double Amount,const double Opacity,const double EdgeControl,const double Masking)3845 ptImage* ptImage::TextureContrast(const double Threshold,
3846 const double Softness,
3847 const double Amount,
3848 const double Opacity,
3849 const double EdgeControl,
3850 const double Masking) {
3851
3852 const int32_t WPH = 0x7fff; // WPH=WP/2
3853 const short NrChannels = (m_ColorSpace == ptSpace_Lab)?1:3;
3854 const TChannelMask ChannelMask = (m_ColorSpace == ptSpace_Lab) ? ChMask_L : ChMask_RGB;
3855
3856 ptImage *ContrastLayer = new ptImage;
3857 ContrastLayer->Set(this);
3858
3859 ContrastLayer->fastBilateralChannel(Threshold, Softness, 2, ChannelMask);
3860
3861 #pragma omp parallel for default(shared) schedule(static)
3862 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
3863 for (short Ch=0; Ch<NrChannels; Ch++) {
3864 ContrastLayer->m_Image[i][Ch] = CLIP((int32_t) ((WPH-(int32_t)ContrastLayer->m_Image[i][Ch])+m_Image[i][Ch]));
3865 if (Amount < 0) ContrastLayer->m_Image[i][Ch] = 0xffff-ContrastLayer->m_Image[i][Ch];
3866 }
3867 }
3868 ContrastLayer->SigmoidalContrast(fabs(Amount), 0.5, ChannelMask);
3869
3870 if (EdgeControl)
3871 ContrastLayer->WaveletDenoise(ChannelMask, EdgeControl, 0.2, 0);
3872
3873 if (Masking) {
3874 ptImage *MaskLayer = new ptImage;
3875 MaskLayer->Set(ContrastLayer);
3876
3877 if (ChannelMask == 7) {
3878 #pragma omp parallel for default(shared) schedule(static)
3879 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++)
3880 MaskLayer->m_Image[i][0] = RGB_2_L(MaskLayer->m_Image[i]);
3881 }
3882
3883 ptCurve* Curve = new ptCurve({TAnchor(0.0, 1.0),
3884 TAnchor(0.4, 0.3),
3885 TAnchor(0.5, 0.0),
3886 TAnchor(0.6, 0.3),
3887 TAnchor(1.0, 1.0)});
3888
3889 MaskLayer->ApplyCurve(Curve,1);
3890 MaskLayer->ptCIBlur(Masking, 1);
3891
3892 Curve->setFromAnchors({TAnchor(0.0, 0.0),
3893 TAnchor(0.4, 0.6),
3894 TAnchor(0.8, 1.0)});
3895
3896 MaskLayer->ApplyCurve(Curve,1);
3897 delete Curve;
3898
3899 float (*Mask);
3900 Mask = MaskLayer->GetMask(TMaskType::Highlights, 0.0, 1.0, 0.0, 1,0,0);
3901 Overlay(ContrastLayer->m_Image,Opacity,Mask);
3902 FREE(Mask);
3903 delete MaskLayer;
3904 } else {
3905 Overlay(ContrastLayer->m_Image,Opacity,NULL);
3906 }
3907
3908 delete ContrastLayer;
3909
3910 return this;
3911 }
3912
3913 ////////////////////////////////////////////////////////////////////////////////
3914 //
3915 // Localcontrast
3916 //
3917 ////////////////////////////////////////////////////////////////////////////////
3918
Localcontrast(const int Radius1,const double Opacity,const double m,const double Feather,const short Method)3919 ptImage* ptImage::Localcontrast(const int Radius1, const double Opacity, const double m, const double Feather, const short Method) {
3920
3921 uint16_t Width = m_Width;
3922 uint16_t Height = m_Height;
3923 int32_t Size = Width*Height;
3924
3925 int MaxExponent=MAX((int)floor(log((float)Radius1)/log(2.))+1,1);
3926 uint16_t MaxRadius=(int)pow(2.,(float)(MaxExponent+1));
3927 int Step = 0;
3928 int Blocksize = 0;
3929 int SquareStep = 0;
3930 int Shift = 0;
3931 int32_t x1 = 0;
3932 int32_t x2 = 0;
3933 int32_t y1 = 0;
3934 int32_t y2 = 0;
3935 int32_t pos1x = 0;
3936 int32_t pos1y = 0;
3937 int32_t pos2x = 0;
3938 int32_t pos2y = 0;
3939 int32_t Index1 = 0;
3940 int32_t Index2 = 0;
3941 float value1 = 0;
3942 float value2 = 0;
3943 float ker = 0;
3944
3945 float **Kernel;
3946 Kernel = (float **)CALLOC(MaxRadius,sizeof(float*));
3947 for(int i = 0; i < MaxRadius; i++) Kernel[i] = (float*)CALLOC(MaxRadius,sizeof(float));
3948
3949 float r=0;
3950 for (uint16_t x=0; x<MaxRadius; x++) {
3951 for (uint16_t y=0; y<MaxRadius; y++) {
3952 r = powf((float)x*(float)x+(float)y*(float)y,0.5);
3953 if (r <= Radius1) {
3954 Kernel[x][y] = 1;
3955 } else {
3956 Kernel[x][y] = 0;
3957 }
3958 }
3959 }
3960 int16_t (*MinVector)[2] = (int16_t (*)[2]) CALLOC(Size,sizeof(*MinVector));
3961 ptMemoryError(MinVector,__FILE__,__LINE__);
3962 int16_t (*MaxVector)[2] = (int16_t (*)[2]) CALLOC(Size,sizeof(*MaxVector));
3963 ptMemoryError(MaxVector,__FILE__,__LINE__);
3964
3965 memset(MinVector,0,Size*sizeof(*MinVector));
3966 memset(MaxVector,0,Size*sizeof(*MaxVector));
3967
3968 for (int runcounter=0; runcounter<1;runcounter++){
3969 // runcounter loop should be superfluous now.
3970 // clean up after some more testing.
3971 for (int k=1; k<=MaxExponent; k++) {
3972 Step=(int) powf(2.,(float)(k));
3973 SquareStep = (int) powf(2.,(float)(k-1));
3974 Blocksize = SquareStep;
3975 Shift = Blocksize;
3976 #pragma omp parallel for private(x1, x2, y1, y2, pos1x, pos1y, pos2x, pos2y, value1, value2, ker) schedule(static)
3977 for (int32_t Row=-(runcounter&1)*Shift; Row < Height; Row += Step) {
3978 for (int32_t Col=-(runcounter&1)*Shift; Col < Width; Col += Step) {
3979 for (uint16_t incx=0; incx < Blocksize; incx++) {
3980 for (uint16_t incy=0; incy < Blocksize; incy++) {
3981 for (int fx=0; fx < 2; fx++) {
3982 for (int fy=0; fy < 2; fy++) {
3983 x1 = MAX(MIN(Col+incx+fx*SquareStep,Width-1),0);
3984 y1 = MAX(MIN(Row+incy+fy*SquareStep,Height-1),0);
3985 for (int lx=0; lx < 2; lx++) {
3986 for (int ly=0; ly < 2; ly++) {
3987 x2 = MAX(MIN(Col+incx+lx*SquareStep,Width-1),0);
3988 y2 = MAX(MIN(Row+incy+ly*SquareStep,Height-1),0);
3989 Index1 = y1*Width+x1;
3990 Index2 = y2*Width+x2;
3991 // Calculate Min
3992 pos1x = MinVector[Index1][0] + x1;
3993 pos1y = MinVector[Index1][1] + y1;
3994 pos2x = MinVector[Index2][0] + x2;
3995 pos2y = MinVector[Index2][1] + y2;
3996 value1 = (float)m_Image[pos1y*Width+pos1x][0];
3997 value2 = (float)m_Image[pos2y*Width+pos2x][0];
3998 ker = Kernel[(int32_t)fabs(pos2x - x1)][(int32_t)fabs(pos2y - y1)];
3999 if (ker && value2 <= value1) {
4000 MinVector[Index1][0] = pos2x - x1;
4001 MinVector[Index1][1] = pos2y - y1;
4002 }
4003 // Calculate Max
4004 pos1x = MaxVector[Index1][0] + x1;
4005 pos1y = MaxVector[Index1][1] + y1;
4006 pos2x = MaxVector[Index2][0] + x2;
4007 pos2y = MaxVector[Index2][1] + y2;
4008 value1 = (float)m_Image[pos1y*Width+pos1x][0];
4009 value2 = (float)m_Image[pos2y*Width+pos2x][0];
4010 ker = Kernel[(int32_t)fabs(pos2x - x1)][(int32_t)fabs(pos2y - y1)];
4011 if (ker && value2 >= value1) {
4012 MaxVector[Index1][0] = pos2x - x1;
4013 MaxVector[Index1][1] = pos2y - y1;
4014 }
4015 }
4016 }
4017 }
4018 }
4019 }
4020 }
4021 }
4022 }
4023 }
4024 }
4025
4026 uint16_t (*MinLayer) = (uint16_t (*)) CALLOC(Size,sizeof(*MinLayer));
4027 ptMemoryError(MinLayer,__FILE__,__LINE__);
4028 uint16_t (*MaxLayer) = (uint16_t (*)) CALLOC(Size,sizeof(*MaxLayer));
4029 ptMemoryError(MaxLayer,__FILE__,__LINE__);
4030
4031 #pragma omp parallel for private(pos1x, pos1y, value1) schedule(static)
4032 for (uint16_t Row=0; Row < Height; Row++) {
4033 for (uint16_t Col=0; Col < Width; Col++) {
4034 int32_t Index = Row*Width+Col;
4035 pos1x = MinVector[Index][0] + Col;
4036 pos1y = MinVector[Index][1] + Row;
4037 value1 = (float)m_Image[pos1y*Width+pos1x][0];
4038 MinLayer[Index] = CLIP((int32_t) value1);
4039 pos1x = MaxVector[Index][0] + Col;
4040 pos1y = MaxVector[Index][1] + Row;
4041 value1 = (float)m_Image[pos1y*Width+pos1x][0];
4042 MaxLayer[Index] = CLIP((int32_t) value1);
4043 }
4044 }
4045
4046 double BlurRadius = Radius1*pow(4,Feather);
4047 ptCimgBlurLayer(MinLayer, Width, Height, BlurRadius);
4048 ptCimgBlurLayer(MaxLayer, Width, Height, BlurRadius);
4049
4050 float Z = 0;
4051 float N = 0;
4052 float Mask = 0;
4053 if (Method == 1) {
4054 #pragma omp parallel for private(Z,N,Mask) schedule(static)
4055 for (int32_t i = 0; i < Size; i++) {
4056 Z = m_Image[i][0] - MinLayer[i];
4057 N = MaxLayer[i] - MinLayer[i];
4058 if (m>0) {
4059 Mask = m*N/(float)0xffff+(1.-m);
4060 } else {
4061 Mask = -m*(1-(N/(float)0xffff))+(1.+m);
4062 }
4063 m_Image[i][0]=CLIP((int32_t) ((Z/N*0xFFFF*Mask+(1-Mask)*m_Image[i][0])*Opacity+(1.-Opacity)*m_Image[i][0]));
4064 }
4065 } else if (Method == 2) {
4066 #pragma omp parallel for private(Z,N,Mask) schedule(static)
4067 for (int32_t i = 0; i < Size; i++) {
4068 m_Image[i][0] = MinLayer[i];
4069 }
4070 } else {
4071 #pragma omp parallel for private(Z,N,Mask) schedule(static)
4072 for (int32_t i = 0; i < Size; i++) {
4073 m_Image[i][0] = MaxLayer[i];
4074 }
4075 }
4076
4077 FREE(MinVector);
4078 FREE(MinLayer);
4079 FREE(MaxVector);
4080 FREE(MaxLayer);
4081 for(int i = 0; i < MaxRadius; i++) FREE(Kernel[i]);
4082 FREE(Kernel);
4083
4084 return this;
4085 }
4086
4087 ////////////////////////////////////////////////////////////////////////////////
4088 //
4089 // Grain
4090 //
4091 ////////////////////////////////////////////////////////////////////////////////
4092
Grain(const double Sigma,const TGrainType NoiseType,const double Radius,const double Opacity,const TMaskType MaskType,const double LowerLimit,const double UpperLimit,const short ScaleFactor)4093 ptImage* ptImage::Grain(const double Sigma, // 0-1
4094 const TGrainType NoiseType, // 0-5, Gaussian, uniform, salt&pepper
4095 const double Radius, // 0-20
4096 const double Opacity,
4097 const TMaskType MaskType,
4098 const double LowerLimit,
4099 const double UpperLimit,
4100 const short ScaleFactor) { // 0, 1 or 2 depending on pipe size
4101
4102 assert (m_ColorSpace == ptSpace_Lab);
4103
4104 ptImage *NoiseLayer = new ptImage;
4105 NoiseLayer->Set(this); // allocation of free layer faster? TODO!
4106 float (*Mask);
4107 short Noise = LIM(static_cast<int>(NoiseType),0,5);
4108 Noise = (Noise > 2) ? (Noise - 3) : Noise;
4109 short ScaledRadius = Radius/powf(2.0,(float)ScaleFactor);
4110
4111 ptCimgNoise(NoiseLayer, Sigma*10000, Noise, ScaledRadius);
4112
4113 const float WPH = 0x7fff;
4114
4115 // adaption to get the same optical impression when rescaled
4116 if (ScaleFactor != 2) {
4117 float m = 4/powf(2.0,(float)ScaleFactor);
4118 float t = (1-m)*WPH;
4119 #pragma omp parallel for schedule(static)
4120 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
4121 NoiseLayer->m_Image[i][0] = CLIP((int32_t)(NoiseLayer->m_Image[i][0]*m+t));
4122 }
4123 }
4124
4125 Mask = GetMask(MaskType, LowerLimit, UpperLimit, 0.0);
4126 if (Noise < 3) {
4127 Overlay(NoiseLayer->m_Image,Opacity,Mask,TOverlayMode::Softlight);
4128 } else {
4129 Overlay(NoiseLayer->m_Image,Opacity,Mask,TOverlayMode::GrainMerge);
4130 }
4131
4132 delete NoiseLayer;
4133 FREE(Mask);
4134 return this;
4135 }
4136
4137 ////////////////////////////////////////////////////////////////////////////////
4138 //
4139 // Black&White Styler
4140 //
4141 ////////////////////////////////////////////////////////////////////////////////
4142
BWStyler(const TBWFilmType FilmType,const TBWColorFilter ColorFilterType,const double MultR,const double MultG,const double MultB,const double Opacity)4143 ptImage* ptImage::BWStyler(
4144 const TBWFilmType FilmType,
4145 const TBWColorFilter ColorFilterType,
4146 const double MultR,
4147 const double MultG,
4148 const double MultB,
4149 const double Opacity)
4150 {
4151 double R = 0,G = 0,B = 0;
4152 double FR = 0, FG = 0, FB = 0;
4153
4154 switch (FilmType) {
4155 case TBWFilmType::LowSensitivity:
4156 R = 0.27;
4157 G = 0.27;
4158 B = 0.46;
4159 break;
4160
4161 case TBWFilmType::HighSensitivity:
4162 R = 0.3;
4163 G = 0.28;
4164 B = 0.42;
4165 break;
4166
4167 case TBWFilmType::Hyperpanchromatic:
4168 R = 0.41;
4169 G = 0.25;
4170 B = 0.34;
4171 break;
4172
4173 case TBWFilmType::Orthochromatic:
4174 R = 0.0;
4175 G = 0.42;
4176 B = 0.58;
4177 break;
4178
4179 //Following values corresponding to http://photographynotes.pbworks.com/bwrecipe
4180 case TBWFilmType::NormalContrast:
4181 R = 0.43;
4182 G = 0.33;
4183 B = 0.3;
4184 break;
4185
4186 case TBWFilmType::HighContrast:
4187 R = 0.4;
4188 G = 0.34;
4189 B = 0.60;
4190 break;
4191
4192 case TBWFilmType::Luminance:
4193 R = 0.3;
4194 G = 0.59;
4195 B = 0.11;
4196 break;
4197
4198 case TBWFilmType::Landscape:
4199 R = 0.66;
4200 G = 0.24;
4201 B = 0.10;
4202 break;
4203
4204 case TBWFilmType::FaceInterior:
4205 R = 0.54;
4206 G = 0.44;
4207 B = 0.12;
4208 break;
4209
4210 case TBWFilmType::ChannelMixer:
4211 R = MultR/(MultR+MultG+MultB);
4212 G = MultG/(MultR+MultG+MultB);
4213 B = MultB/(MultR+MultG+MultB);
4214 break;
4215 }
4216
4217 switch (ColorFilterType) {
4218 //Dr. Herbert Kribben and http://epaperpress.com/psphoto/bawFilters.html
4219 case TBWColorFilter::None:
4220 FR = 1.0;
4221 FG = 1.0;
4222 FB = 1.0;
4223 break;
4224
4225 case TBWColorFilter::Red:
4226 FR = 1.0;
4227 FG = 0.2;
4228 FB = 0.0;
4229 break;
4230
4231 case TBWColorFilter::Orange:
4232 FR = 1.0;
4233 FG = 0.6;
4234 FB = 0.0;
4235 break;
4236
4237 case TBWColorFilter::Yellow:
4238 FR = 1.0;
4239 FG = 1.0;
4240 FB = 0.1;
4241 break;
4242
4243 case TBWColorFilter::Lime:
4244 FR = 0.6;
4245 FG = 1.0;
4246 FB = 0.3;
4247 break;
4248
4249 case TBWColorFilter::Green:
4250 FR = 0.2;
4251 FG = 1.0;
4252 FB = 0.3;
4253 break;
4254
4255 case TBWColorFilter::Blue:
4256 FR = 0.0;
4257 FG = 0.2;
4258 FB = 1.0;
4259 break;
4260
4261 case TBWColorFilter::FakeIR:
4262 FR = 0.4;
4263 FG = 1.4;
4264 FB = -0.8;
4265 break;
4266 }
4267
4268 R = R * FR;
4269 G = G * FG;
4270 B = B * FB;
4271 R = R / (R+G+B);
4272 G = G / (R+G+B);
4273 B = B / (R+G+B);
4274
4275 TChannelMatrix Mixer;
4276 Mixer[0][0] = R;
4277 Mixer[1][0] = R;
4278 Mixer[2][0] = R;
4279 Mixer[0][1] = G;
4280 Mixer[1][1] = G;
4281 Mixer[2][1] = G;
4282 Mixer[0][2] = B;
4283 Mixer[1][2] = B;
4284 Mixer[2][2] = B;
4285
4286 if (qFuzzyCompare(Opacity, 1.0)) {
4287 mixChannels(Mixer);
4288 } else {
4289 auto BWLayer = make_unique<ptImage>();
4290 BWLayer->Set(this);
4291 BWLayer->mixChannels(Mixer);
4292 Overlay(BWLayer->m_Image, Opacity, nullptr, TOverlayMode::Normal);
4293 }
4294
4295 return this;
4296 }
4297
4298 ////////////////////////////////////////////////////////////////////////////////
4299 //
4300 // Simple toneing
4301 //
4302 ////////////////////////////////////////////////////////////////////////////////
4303
SimpleTone(const double R,const double G,const double B)4304 ptImage* ptImage::SimpleTone(const double R,
4305 const double G,
4306 const double B) {
4307
4308 assert (m_ColorSpace != ptSpace_Lab);
4309
4310 if (R) {
4311 ptCurve* RedCurve = new ptCurve({TAnchor(0.0, 0.0),
4312 TAnchor(0.5-0.2*R, 0.5+0.2*R),
4313 TAnchor(1.0, 1.0),});
4314 ApplyCurve(RedCurve,1);
4315 delete RedCurve;
4316 }
4317
4318 if (G) {
4319 ptCurve* GreenCurve = new ptCurve({TAnchor(0.0, 0.0),
4320 TAnchor(0.5-0.2*G, 0.5+0.2*G),
4321 TAnchor(1.0, 1.0),});
4322 ApplyCurve(GreenCurve,2);
4323 delete GreenCurve;
4324 }
4325
4326 if (B) {
4327 ptCurve* BlueCurve = new ptCurve({TAnchor(0.0, 0.0),
4328 TAnchor(0.5-0.2*B, 0.5+0.2*B),
4329 TAnchor(1.0, 1.0),});
4330 ApplyCurve(BlueCurve,4);
4331 delete BlueCurve;
4332 }
4333
4334 return this;
4335 }
4336
4337 ////////////////////////////////////////////////////////////////////////////////
4338 //
4339 // Temperature
4340 //
4341 ////////////////////////////////////////////////////////////////////////////////
4342
Temperature(const double Temperature,const double Tint)4343 ptImage* ptImage::Temperature(const double Temperature, const double Tint) {
4344 assert (m_ColorSpace == ptSpace_Lab);
4345
4346 ptCurve* Temp1Curve = new ptCurve(
4347 {TAnchor(0.0, 0.0),
4348 TAnchor(0.5-0.05*Temperature+0.1*Tint, 0.5+0.05*Temperature-0.1*Tint),
4349 TAnchor(1.0, 1.0)} );
4350 ApplyCurve(Temp1Curve,2);
4351 delete Temp1Curve;
4352
4353 ptCurve* Temp2Curve = new ptCurve(
4354 {TAnchor(0.0, 0.0),
4355 TAnchor(0.5-0.1*Temperature-0.05*Tint, 0.5+0.1*Temperature+0.05*Tint),
4356 TAnchor(1.0, 1.0)} );
4357 ApplyCurve(Temp2Curve,4);
4358 delete Temp2Curve;
4359
4360 return this;
4361 }
4362
4363 ////////////////////////////////////////////////////////////////////////////////
4364 //
4365 // LAB Transform
4366 //
4367 ////////////////////////////////////////////////////////////////////////////////
4368
LABTransform(const short Mode)4369 ptImage* ptImage::LABTransform(const short Mode) {
4370
4371 assert (m_ColorSpace != ptSpace_Lab);
4372
4373 double R = 0;
4374 double G = 0;
4375 double B = 0;
4376
4377 switch (Mode) {
4378 case ptLABTransform_R:
4379 R = 1.;
4380 break;
4381 case ptLABTransform_G:
4382 G = 1.;
4383 break;
4384 case ptLABTransform_B:
4385 B = 1.;
4386 break;
4387 default:
4388 break;
4389 }
4390
4391 ptImage *TempLayer = new ptImage;
4392 TempLayer->Set(this);
4393 ptCurve* GammaCurve = new ptCurve();
4394 GammaCurve->setFromFunc(ptCurve::GammaTool,0.45,0.0);
4395 TempLayer->ApplyCurve(GammaCurve,7);
4396 delete GammaCurve;
4397 RGBToLab();
4398
4399 #pragma omp parallel for
4400 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++)
4401 m_Image[i][0] = CLIP((int32_t) (R*(double)TempLayer->m_Image[i][0] +
4402 G*(double)TempLayer->m_Image[i][1] +
4403 B*(double)TempLayer->m_Image[i][2]));
4404 delete TempLayer;
4405 return this;
4406 }
4407
4408 ////////////////////////////////////////////////////////////////////////////////
4409 //
4410 // LAB Tone
4411 //
4412 ////////////////////////////////////////////////////////////////////////////////
4413
LABTone(const double Amount,const double Hue,const double Saturation,const TMaskType MaskType,const short ManualMask,const double LowerLevel,const double UpperLevel,const double Softness)4414 ptImage* ptImage::LABTone(const double Amount,
4415 const double Hue,
4416 const double Saturation, /* 0 */
4417 const TMaskType MaskType, /* TMaskType::All */
4418 const short ManualMask, /* 0 */
4419 const double LowerLevel, /* 0 */
4420 const double UpperLevel, /* 1 */
4421 const double Softness /* 0 */) {
4422
4423 assert (m_ColorSpace == ptSpace_Lab);
4424
4425 double a = Amount * cos(Hue/180.*ptPI);
4426 double b = Amount * sin(Hue/180.*ptPI);
4427
4428 const double WPH = 0x8080; // Neutral in Lab A and B
4429
4430 double t = (1.-Saturation)*WPH;
4431
4432 ptCurve* Temp1Curve = new ptCurve({TAnchor(0.0, 0.0),
4433 TAnchor(0.5-0.1*a, 0.5+0.1*a),
4434 TAnchor(1.0, 1.0)} );
4435
4436 ptCurve* Temp2Curve = new ptCurve({TAnchor(0.0, 0.0),
4437 TAnchor(0.5-0.1*b, 0.5+0.1*b),
4438 TAnchor(1.0, 1.0)} );
4439
4440 if (MaskType == TMaskType::All) {
4441 if (Saturation != 1) ColorBoost(Saturation, Saturation);
4442 if (Amount) {
4443 ApplyCurve(Temp1Curve,2);
4444 ApplyCurve(Temp2Curve,4);
4445 }
4446 } else {
4447 ptImage *ToneLayer = new ptImage;
4448 if (Amount) {
4449 ToneLayer->Set(this);
4450 if (Saturation != 1) ToneLayer->ColorBoost(Saturation, Saturation);
4451 ToneLayer->ApplyCurve(Temp1Curve,2);
4452 ToneLayer->ApplyCurve(Temp2Curve,4);
4453 }
4454 float (*Mask);
4455
4456 if (MaskType == TMaskType::Shadows)
4457 if (ManualMask)
4458 Mask = GetMask(TMaskType::Shadows, LowerLevel, UpperLevel, Softness,1,0,0);
4459 else
4460 Mask = GetMask(TMaskType::Shadows, 0,0.5,0,1,0,0);
4461 else if (MaskType == TMaskType::Midtones)
4462 if (ManualMask)
4463 Mask = GetMask(TMaskType::Midtones, LowerLevel, UpperLevel, Softness,1,0,0);
4464 else
4465 Mask = GetMask(TMaskType::Midtones, 0.5,0.5,0,1,0,0);
4466 else
4467 if (ManualMask)
4468 Mask = GetMask(TMaskType::Highlights, LowerLevel, UpperLevel, Softness,1,0,0);
4469 else
4470 Mask = GetMask(TMaskType::Highlights, 0.5,1,0,1,0,0);
4471
4472 #pragma omp parallel for
4473 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++)
4474 for (int Ch = 1; Ch < 3; Ch++) {
4475 if (Saturation != 1)
4476 m_Image[i][Ch] = CLIP((int32_t)((m_Image[i][Ch]*Saturation+t)*Mask[i] +
4477 m_Image[i][Ch]*(1-Mask[i])));
4478 if (Amount)
4479 m_Image[i][Ch] = CLIP((int32_t)((ToneLayer->m_Image[i][Ch]*Mask[i] +
4480 m_Image[i][Ch]*(1-Mask[i]))*Amount+m_Image[i][Ch]*(1-Amount)));
4481 }
4482 FREE(Mask);
4483 if (Amount) delete ToneLayer;
4484 }
4485 delete Temp1Curve;
4486 delete Temp2Curve;
4487
4488 return this;
4489 }
4490
4491 ////////////////////////////////////////////////////////////////////////////////
4492 //
4493 // Tone
4494 //
4495 ////////////////////////////////////////////////////////////////////////////////
4496
Tone(uint16_t Red,uint16_t Green,uint16_t Blue,const double Amount,const TMaskType MaskType,double LowerLimit,double UpperLimit,const double Softness,const double InputPowerFactor)4497 ptImage* ptImage::Tone(
4498 uint16_t Red,
4499 uint16_t Green,
4500 uint16_t Blue,
4501 const double Amount,
4502 const TMaskType MaskType,
4503 double LowerLimit,
4504 double UpperLimit,
4505 const double Softness,
4506 const double InputPowerFactor)
4507 {
4508 assert (m_ColorSpace != ptSpace_Lab);
4509
4510 Red = static_cast<uint16_t>(0xffff * pow(Red/255.0, InputPowerFactor));
4511 Green = static_cast<uint16_t>(0xffff * pow(Green/255.0, InputPowerFactor));
4512 Blue = static_cast<uint16_t>(0xffff * pow(Blue/255.0, InputPowerFactor));
4513
4514 LowerLimit = pow(LowerLimit, InputPowerFactor);
4515 UpperLimit = pow(UpperLimit, InputPowerFactor);
4516
4517 uint16_t (*ToneImage)[3] =
4518 (uint16_t (*)[3]) CALLOC(m_Width*m_Height,sizeof(*ToneImage));
4519 ptMemoryError(ToneImage,__FILE__,__LINE__);
4520
4521 float (*Mask);
4522 if (MaskType <= TMaskType::All)
4523 Mask=GetMask(MaskType, LowerLimit, UpperLimit, Softness);
4524 else
4525 Mask=GetMask(TMaskType::Midtones, LowerLimit, UpperLimit, Softness);
4526 #pragma omp parallel for
4527 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
4528 ToneImage[i][0] = Red;
4529 ToneImage[i][1] = Green;
4530 ToneImage[i][2] = Blue;
4531 }
4532
4533 if (MaskType <= TMaskType::All)
4534 Overlay(ToneImage, Amount, Mask);
4535 else if (MaskType == TMaskType::Screen)
4536 Overlay(ToneImage, Amount, Mask, TOverlayMode::Screen);
4537 else if (MaskType == TMaskType::Multiply)
4538 Overlay(ToneImage, Amount, Mask, TOverlayMode::Multiply);
4539 else if (MaskType == TMaskType::GammaDark)
4540 Overlay(ToneImage, Amount, Mask, TOverlayMode::GammaDark);
4541 else if (MaskType == TMaskType::GammaBright)
4542 Overlay(ToneImage, Amount, Mask, TOverlayMode::GammaBright);
4543
4544 FREE(Mask);
4545 FREE(ToneImage);
4546
4547 return this;
4548 }
4549
4550 ////////////////////////////////////////////////////////////////////////////////
4551 //
4552 // Crossprocessing
4553 //
4554 ////////////////////////////////////////////////////////////////////////////////
4555
Crossprocess(const TCrossProcessMode mode,const double color1Intensity,const double color2Intensity)4556 ptImage* ptImage::Crossprocess(
4557 const TCrossProcessMode mode,
4558 const double color1Intensity,
4559 const double color2Intensity)
4560 {
4561 assert (m_ColorSpace != ptSpace_Lab);
4562
4563 ptCurve* RedCurve = new ptCurve({TAnchor(0.0, 0.0),
4564 TAnchor(0.05, 0.05-0.05*color1Intensity),
4565 TAnchor(0.15, 0.15),
4566 TAnchor(1.0, 1.0)});
4567
4568 ptCurve* GreenCurve = new ptCurve({TAnchor(0.0, 0.0),
4569 TAnchor(0.03, 0.025),
4570 TAnchor(0.2, 0.2+0.3*color1Intensity),
4571 TAnchor(1.0, 1.0)});
4572
4573 ptCurve* BlueCurve = new ptCurve({TAnchor(0.0, 0.0),
4574 TAnchor(0.3, 0.3-0.2*color2Intensity),
4575 TAnchor(1.0, 1.0)});
4576
4577 ptImage *ColorLayer = new ptImage;
4578 ColorLayer->Set(this);
4579
4580 switch (mode) {
4581 case TCrossProcessMode::GreenYellow:
4582 ColorLayer->ApplyCurve(RedCurve,1);
4583 ColorLayer->ApplyCurve(GreenCurve,2);
4584 ColorLayer->ApplyCurve(BlueCurve,4);
4585 break;
4586 case TCrossProcessMode::GreenCyan:
4587 ColorLayer->ApplyCurve(RedCurve,4);
4588 ColorLayer->ApplyCurve(GreenCurve,2);
4589 ColorLayer->ApplyCurve(BlueCurve,1);
4590 break;
4591 case TCrossProcessMode::RedYellow:
4592 ColorLayer->ApplyCurve(RedCurve,2);
4593 ColorLayer->ApplyCurve(GreenCurve,1);
4594 ColorLayer->ApplyCurve(BlueCurve,4);
4595 break;
4596 case TCrossProcessMode::RedMagenta:
4597 ColorLayer->ApplyCurve(RedCurve,4);
4598 ColorLayer->ApplyCurve(GreenCurve,1);
4599 ColorLayer->ApplyCurve(BlueCurve,2);
4600 break;
4601 case TCrossProcessMode::BlueCyan:
4602 ColorLayer->ApplyCurve(RedCurve,2);
4603 ColorLayer->ApplyCurve(GreenCurve,4);
4604 ColorLayer->ApplyCurve(BlueCurve,1);
4605 break;
4606 case TCrossProcessMode::BlueMagenta:
4607 ColorLayer->ApplyCurve(RedCurve,1);
4608 ColorLayer->ApplyCurve(GreenCurve,4);
4609 ColorLayer->ApplyCurve(BlueCurve,2);
4610 break;
4611 default:
4612 assert("Unknown TCrossProcessMode");
4613 }
4614
4615 Overlay(ColorLayer->m_Image,0.7,NULL,TOverlayMode::Normal);
4616 delete ColorLayer;
4617
4618 delete RedCurve;
4619 delete GreenCurve;
4620 delete BlueCurve;
4621
4622 return this;
4623 }
4624
4625 ////////////////////////////////////////////////////////////////////////////////
4626 //
4627 // Gradual Overlay Mask
4628 //
4629 ////////////////////////////////////////////////////////////////////////////////
4630
GetGradualMask(const double Angle,const double LowerLevel,const double UpperLevel,const double Softness)4631 pt::c_unique_ptr<float> ptImage::GetGradualMask(
4632 const double Angle,
4633 const double LowerLevel,
4634 const double UpperLevel,
4635 const double Softness)
4636 {
4637 pt::c_unique_ptr<float> GradualMask((float (*)) CALLOC(m_Width*m_Height,sizeof(float)));
4638 ptMemoryError(GradualMask.get(),__FILE__,__LINE__);
4639
4640 float Length = 0;
4641 if (fabs(Angle) == 0 || fabs(Angle) == 180 ) {
4642 Length = m_Height;
4643 } else if (fabs(Angle) == 90) {
4644 Length = m_Width;
4645 } else if (fabs(Angle) < 90) {
4646 Length = (((float)m_Width) + ((float)m_Height)/tan(fabs(Angle)/180*ptPI))*sin(fabs(Angle)/180*ptPI);
4647 } else {
4648 Length = (((float)m_Width) + ((float)m_Height)/tan((180.0-fabs(Angle))/180.0*ptPI))*sin((180.0-fabs(Angle))/180.0*ptPI);
4649 }
4650
4651 bool Switch = UpperLevel < LowerLevel;
4652
4653 float Eps = 0.0001f;
4654
4655 float LL = Length*(Switch?UpperLevel:LowerLevel);
4656 float UL = Length*(Switch?LowerLevel:UpperLevel);
4657 float Black = Switch?1.0f:0.0f;
4658 float White = Switch?0.0f:1.0f;
4659 float Denom = 1.0f/MAX((UL-LL),Eps);
4660
4661 float coordinate = 0;
4662 float Value = 0;
4663
4664 float dist = 0;
4665
4666 float Factor1 = 0;
4667 float Factor2 = 0;
4668 if (Angle >= 0.0 && Angle < 90.0) {
4669 Factor1 = 1.0/MAX(tanf(Angle/180*ptPI),Eps);
4670 Factor2 = sinf(Angle/180*ptPI);
4671 } else if (Angle >= 90.0 && Angle < 180.0) {
4672 Factor1 = 1.0/MAX(tanf((180.0-Angle)/180*ptPI),Eps);
4673 Factor2 = sinf((180.0-Angle)/180*ptPI);
4674 } else if (Angle >= -90.0 && Angle < 0.0) {
4675 Factor1 = 1.0/MAX(tanf(fabs(Angle)/180*ptPI),Eps);
4676 Factor2 = sinf(fabs(Angle)/180*ptPI);
4677 } else if (Angle >= -180.0 && Angle < -90.0) {
4678 Factor1 = 1.0/MAX(tanf((180.0-fabs(Angle))/180*ptPI),Eps);
4679 Factor2 = sinf((180.0-fabs(Angle))/180*ptPI);
4680 }
4681
4682 #pragma omp parallel for schedule(static) firstprivate(dist, Value, coordinate)
4683 for (uint16_t Row=0; Row<m_Height; Row++) {
4684 uint32_t Idx = Row*m_Width;
4685 for (uint16_t Col=0; Col<m_Width; Col++) {
4686 if (fabs(Angle) == 0.0)
4687 dist = m_Height-Row;
4688 else if (fabs(Angle) == 180.0)
4689 dist = Row;
4690 else if (Angle == 90.0)
4691 dist = Col;
4692 else if (Angle == -90.0)
4693 dist = m_Width-Col;
4694 else if (Angle > 0.0 && Angle < 90.0)
4695 dist = (Col + (float)(m_Height-Row)*Factor1)*Factor2;
4696 else if (Angle > 90.0 && Angle < 180.0)
4697 dist = Length-((m_Width-Col) + (float)(m_Height-Row)*Factor1)*Factor2;
4698 else if (Angle > -90.0 && Angle < 0.0)
4699 dist = ((m_Width-Col) + (float)(m_Height-Row)*Factor1)*Factor2;
4700 else if (Angle > -180.0 && Angle < -90.0)
4701 dist = Length-((float)Col + (float)(m_Height-Row)*Factor1)*Factor2;
4702
4703 if (dist <= LL)
4704 GradualMask.get()[Idx+Col] = Black;
4705 else if (dist >= UL)
4706 GradualMask.get()[Idx+Col] = White;
4707 else {
4708 coordinate = 1.0f - (UL-dist)*Denom;
4709 Value = (1.0f-powf(cosf(coordinate*ptPI/2.0f),50.0f*Softness))
4710 * powf(coordinate,0.07f*Softness);
4711 GradualMask.get()[Idx+Col] = LIM(Value*White, 0.0f, 1.0f);
4712 }
4713 }
4714 }
4715
4716 return GradualMask;
4717 }
4718
4719 ////////////////////////////////////////////////////////////////////////////////
4720 //
4721 // Gradual Overlay
4722 //
4723 ////////////////////////////////////////////////////////////////////////////////
4724
GradualOverlay(const uint16_t R,const uint16_t G,const uint16_t B,const TOverlayMode Mode,const double Amount,const double Angle,const double LowerLevel,const double UpperLevel,const double Softness)4725 ptImage* ptImage::GradualOverlay(
4726 const uint16_t R,
4727 const uint16_t G,
4728 const uint16_t B,
4729 const TOverlayMode Mode,
4730 const double Amount,
4731 const double Angle,
4732 const double LowerLevel,
4733 const double UpperLevel,
4734 const double Softness)
4735 {
4736 auto GradualMask = GetGradualMask(Angle, LowerLevel, UpperLevel, Softness);
4737
4738 uint16_t (*ToneImage)[3] = (uint16_t (*)[3]) CALLOC(m_Width*m_Height,sizeof(*ToneImage));
4739 ptMemoryError(ToneImage,__FILE__,__LINE__);
4740
4741 # pragma omp for schedule(static)
4742 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
4743 ToneImage[i][0] = R;
4744 ToneImage[i][1] = G;
4745 ToneImage[i][2] = B;
4746 }
4747
4748 Overlay(ToneImage, Amount, GradualMask.get(), Mode);
4749
4750 FREE(ToneImage);
4751 return this;
4752 }
4753
4754 ////////////////////////////////////////////////////////////////////////////////
4755 //
4756 // Vignette
4757 //
4758 ////////////////////////////////////////////////////////////////////////////////
4759
Vignette(const TVignetteMask VignetteMode,const TVignetteShape Shape,const double Amount,const double InnerRadius,const double OuterRadius,const double Roundness,const double CenterX,const double CenterY,const double Softness)4760 ptImage* ptImage::Vignette(const TVignetteMask VignetteMode,
4761 const TVignetteShape Shape,
4762 const double Amount,
4763 const double InnerRadius,
4764 const double OuterRadius,
4765 const double Roundness,
4766 const double CenterX,
4767 const double CenterY,
4768 const double Softness)
4769 {
4770 auto VignetteMask = GetVignetteMask(
4771 false,
4772 Shape,
4773 InnerRadius,
4774 OuterRadius,
4775 Roundness,
4776 CenterX,
4777 CenterY,
4778 Softness);
4779
4780 switch (VignetteMode) {
4781 case TVignetteMask::Disabled:
4782 assert("Should not be called with disabled vignette mask!");
4783 break;
4784
4785 case TVignetteMask::Soft:
4786 case TVignetteMask::Hard:
4787 {
4788 uint16_t (*ToneImage)[3] = (uint16_t (*)[3]) CALLOC(m_Width*m_Height,sizeof(*ToneImage));
4789 ptMemoryError(ToneImage,__FILE__,__LINE__);
4790 uint16_t hColor = 0;
4791 auto hMode = TOverlayMode::Softlight;
4792 if (Amount > 0) {
4793 hColor = 0;
4794 if (VignetteMode == TVignetteMask::Hard) hMode = TOverlayMode::Multiply;
4795 } else {
4796 hColor = 0xffff;
4797 if (VignetteMode == TVignetteMask::Hard) hMode = TOverlayMode::Screen;
4798 }
4799
4800 #pragma omp parallel for schedule(static) default(shared)
4801 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
4802 ToneImage[i][0] = ToneImage[i][1] = ToneImage[i][2] = hColor;
4803 }
4804
4805 Overlay(ToneImage, fabs(Amount), VignetteMask.get(), hMode);
4806 FREE(ToneImage);
4807 }
4808 break;
4809
4810 case TVignetteMask::Fancy:
4811 {
4812 ptImage *VignetteLayer = new ptImage;
4813 VignetteLayer->Set(this);
4814 VignetteLayer->Expose(pow(2,-Amount*5), TExposureClipMode::Hard);
4815 ptCurve* VignetteContrastCurve = new ptCurve();
4816 VignetteContrastCurve->setFromFunc(ptCurve::Sigmoidal,0.5,fabs(Amount)*10);
4817 VignetteLayer->ApplyCurve(VignetteContrastCurve, (m_ColorSpace == ptSpace_Lab) ? 1 : 7);
4818 delete VignetteContrastCurve;
4819 Overlay(VignetteLayer->m_Image, 1, VignetteMask.get(), TOverlayMode::Normal);
4820 delete VignetteLayer;
4821 }
4822 break;
4823
4824 case TVignetteMask::ShowMask:
4825 if (m_ColorSpace == ptSpace_Lab) {
4826 # pragma omp parallel for schedule(static) default(shared)
4827 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
4828 m_Image[i][0] = CLIP((int32_t) (VignetteMask.get()[i]*0xffff));
4829 m_Image[i][1] = m_Image[i][2] = 0x8080;
4830 }
4831 } else {
4832 # pragma omp parallel for schedule(static) default(shared)
4833 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
4834 m_Image[i][0] = m_Image[i][1] = m_Image[i][2] =
4835 CLIP((int32_t) (VignetteMask.get()[i]*0xffff));
4836 }
4837 }
4838 break;
4839 }
4840
4841 return this;
4842 }
4843
4844 ////////////////////////////////////////////////////////////////////////////////
4845 //
4846 // Softglow
4847 //
4848 ////////////////////////////////////////////////////////////////////////////////
4849
Softglow(const TSoftglowMode SoftglowMode,const double Radius,const double Amount,const uint8_t ChannelMask,const double Contrast,const int Saturation)4850 ptImage* ptImage::Softglow(
4851 const TSoftglowMode SoftglowMode,
4852 const double Radius,
4853 const double Amount,
4854 const uint8_t ChannelMask,
4855 const double Contrast, /* 5 */
4856 const int Saturation /* -50 */)
4857 {
4858 // Workflow for Orton from
4859 // http://www.flickr.com/groups/gimpusers/discuss/72157621381832321/
4860 ptImage *BlurLayer = new ptImage;
4861 BlurLayer->Set(this);
4862 // for Orton
4863 if (SoftglowMode==TSoftglowMode::OrtonScreen || SoftglowMode==TSoftglowMode::OrtonSoftlight) {
4864 BlurLayer->Overlay(BlurLayer->m_Image,1.0, nullptr, TOverlayMode::Screen);
4865 }
4866
4867 // Blur
4868 if (Radius != 0)
4869 BlurLayer->ptCIBlur(Radius, ChannelMask);
4870 // Contrast
4871 if (Contrast != 0) {
4872 ptCurve* MyContrastCurve = new ptCurve();
4873 MyContrastCurve->setFromFunc(ptCurve::Sigmoidal,0.5,Contrast);
4874 BlurLayer->ApplyCurve(MyContrastCurve,7);
4875 delete MyContrastCurve;
4876 }
4877 // Desaturate
4878 if (Saturation != 0) {
4879 int Value = Saturation;
4880 TChannelMatrix VibranceMixer;
4881 VibranceMixer[0][0] = 1.0+(Value/150.0);
4882 VibranceMixer[0][1] = -(Value/300.0);
4883 VibranceMixer[0][2] = VibranceMixer[0][1];
4884 VibranceMixer[1][0] = VibranceMixer[0][1];
4885 VibranceMixer[1][1] = VibranceMixer[0][0];
4886 VibranceMixer[1][2] = VibranceMixer[0][1];
4887 VibranceMixer[2][0] = VibranceMixer[0][1];
4888 VibranceMixer[2][1] = VibranceMixer[0][1];
4889 VibranceMixer[2][2] = VibranceMixer[0][0];
4890 BlurLayer->mixChannels(VibranceMixer);
4891 }
4892 // Overlay
4893 switch (SoftglowMode) {
4894 case TSoftglowMode::Lighten:
4895 Overlay(BlurLayer->m_Image,Amount,nullptr,TOverlayMode::Lighten);
4896 break;
4897 case TSoftglowMode::Screen:
4898 Overlay(BlurLayer->m_Image,Amount,nullptr,TOverlayMode::Screen);
4899 break;
4900 case TSoftglowMode::Softlight:
4901 Overlay(BlurLayer->m_Image,Amount,nullptr,TOverlayMode::Softlight);
4902 break;
4903 case TSoftglowMode::Normal:
4904 if (Amount != 0)
4905 Overlay(BlurLayer->m_Image,Amount,nullptr,TOverlayMode::Normal);
4906 break;
4907 case TSoftglowMode::OrtonScreen:
4908 Overlay(BlurLayer->m_Image,Amount,nullptr,TOverlayMode::Screen);
4909 break;
4910 case TSoftglowMode::OrtonSoftlight:
4911 Overlay(BlurLayer->m_Image,Amount,nullptr,TOverlayMode::Softlight);
4912 break;
4913 default:
4914 break;
4915 }
4916 delete BlurLayer;
4917 return this;
4918
4919 }
4920
4921 ////////////////////////////////////////////////////////////////////////////////
4922 //
4923 // GetMask
4924 //
4925 ////////////////////////////////////////////////////////////////////////////////
4926
GetMask(const TMaskType MaskType,const double LowerLimit,const double UpperLimit,const double Softness,const double FactorR,const double FactorG,const double FactorB)4927 float *ptImage::GetMask(const TMaskType MaskType,
4928 const double LowerLimit,
4929 const double UpperLimit,
4930 const double Softness,
4931 const double FactorR,
4932 const double FactorG,
4933 const double FactorB) {
4934
4935 float (*dMask) = (float (*)) CALLOC(m_Width*m_Height,sizeof(*dMask));
4936 ptMemoryError(dMask,__FILE__,__LINE__);
4937
4938 const float WP = 0xffff;
4939 float m = 1.0/(UpperLimit-LowerLimit);
4940 float t = -LowerLimit/(UpperLimit-LowerLimit)*WP;
4941 float m1 = -1.0/MAX(LowerLimit,0.001);
4942 float m2 = 1.0/MAX(0.001,(1.0-UpperLimit));
4943 float t2 = -(UpperLimit)/MAX(0.001,(1.0-UpperLimit))*WP;
4944 float Soft = pow(2,Softness);
4945
4946 // Precalculated table for the mask.
4947 float MaskTable[0x10000];
4948 float FactorRTable[0x10000];
4949 float FactorGTable[0x10000];
4950 float FactorBTable[0x10000];
4951 #pragma omp parallel
4952 { // begin OpenMP
4953 #pragma omp for schedule(static)
4954 for (int32_t i=0; i<0x10000; i++) {
4955 switch(MaskType) {
4956 case TMaskType::All: // All values
4957 MaskTable[i] = WP;
4958 break;
4959 case TMaskType::Shadows: // Shadows
4960 MaskTable[i] = WP - CLIP((int32_t)(i*m+t));
4961 break;
4962 case TMaskType::Midtones: // Midtones
4963 MaskTable[i] = WP - CLIP((int32_t)(i*m1+WP)) - CLIP((int32_t)(i*m2+t2));
4964 break;
4965 case TMaskType::Highlights: // Highlights
4966 MaskTable[i] = CLIP((int32_t)(i*m+t));
4967 break;
4968 default:
4969 assert("Unexpected mask type");
4970 }
4971 if (Soft>1.0) {
4972 MaskTable[i] = ptBound((float)(pow(MaskTable[i]/0xffff,Soft)), 0.0f, 1.0f);
4973 } else {
4974 MaskTable[i] = ptBound((float)(MaskTable[i]/0xffff/Soft), 0.0f, 1.0f);
4975 }
4976 FactorRTable[i] = i*FactorR;
4977 FactorGTable[i] = i*FactorG;
4978 FactorBTable[i] = i*FactorB;
4979 }
4980
4981 #pragma omp for schedule(static)
4982 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
4983 dMask[i] = MaskTable[CLIP((int32_t)
4984 (FactorRTable[m_Image[i][0]] +
4985 FactorGTable[m_Image[i][1]] +
4986 FactorBTable[m_Image[i][2]]))];
4987 }
4988 } // end OpenMP
4989
4990 return dMask;
4991 }
4992
4993 //==============================================================================
4994
4995 // InMask should only be called once per pixel!
fill4stack(float * AMask,const uint16_t APointX,const uint16_t APointY,const uint16_t AWidth,const uint16_t AHeight,std::function<float (uint16_t,uint16_t)> InMask)4996 void fill4stack(float *AMask,
4997 const uint16_t APointX,
4998 const uint16_t APointY,
4999 const uint16_t AWidth,
5000 const uint16_t AHeight,
5001 std::function<float (uint16_t, uint16_t)> InMask){
5002
5003 std::stack<uint16_t> hStack;
5004 hStack.push(APointX);
5005 hStack.push(APointY);
5006
5007 uint16_t hX = 0,
5008 hY = 0;
5009 float hValue = 0.0f;
5010
5011 while (!hStack.empty()) {
5012 hY = hStack.top();
5013 hStack.pop();
5014 hX = hStack.top();
5015 hStack.pop();
5016
5017 // array bounds
5018 // hX, hY >= 0 by design
5019 if (hX >= AWidth || hY >= AHeight ) continue;
5020
5021 // already processed?
5022 if (AMask[hY*AWidth + hX] != 0.0f) continue;
5023
5024 hValue = InMask(hX, hY);
5025 // pixel suitable for mask?
5026 if (hValue > 0.0f) {
5027 AMask[hY*AWidth + hX] = hValue;
5028
5029 hStack.push(hX);
5030 hStack.push(hY + 1);
5031
5032 if (hY > 0) {
5033 hStack.push(hX);
5034 hStack.push(hY - 1);
5035 }
5036
5037 hStack.push(hX + 1);
5038 hStack.push(hY);
5039
5040 if (hX > 0) {
5041 hStack.push(hX - 1);
5042 hStack.push(hY);
5043 }
5044 }
5045 }
5046 }
5047
5048 //==============================================================================
5049
FillMask(const uint16_t APointX,const uint16_t APointY,const float AThreshold,const float AColorWeight,const uint16_t AMaxRadius,const bool AUseMaxRadius)5050 float *ptImage::FillMask(const uint16_t APointX,
5051 const uint16_t APointY,
5052 const float AThreshold,
5053 const float AColorWeight,
5054 const uint16_t AMaxRadius,
5055 const bool AUseMaxRadius)
5056 {
5057 assert (m_ColorSpace == ptSpace_LCH);
5058
5059 float (*FillMask) = (float (*)) CALLOC(m_Width*m_Height,sizeof(*FillMask));
5060 ptMemoryError(FillMask,__FILE__,__LINE__);
5061
5062 memset(FillMask, 0, m_Width*m_Height*sizeof(*FillMask));
5063
5064 float hThresholdHalf = AThreshold*0x2AAA;
5065 float hThreshold = AThreshold*0x5555;
5066 float hRadiusOut = ptSqr((float)AMaxRadius);
5067 float hLumaWeight = 1.0f - AColorWeight;
5068 float hColorWeight = AColorWeight*(float)0x7FFF;
5069
5070 float hValueL = 0.0f,
5071 hValueC = 0.0f,
5072 hValueH = 0.0f;
5073 float hDiff = 0;
5074 int32_t hIdx = 0;
5075 float hRad = 0;
5076
5077 // Calculate mean around the sample point
5078 short hSample = 1;
5079 int32_t hCnt = 0;
5080
5081 for (int32_t lX = APointX-hSample; lX <= APointX+hSample; lX++) {
5082 if (lX < 0 || lX >= m_Width) continue;
5083
5084 for (int32_t lY = APointY-hSample; lY <= APointY+hSample; lY++) {
5085 if (lY < 0 || lY >= m_Height) continue;
5086
5087 hIdx = (int32_t)lY*m_Width + lX;
5088 hValueL += m_ImageL[hIdx];
5089 hValueC += m_ImageC[hIdx];
5090 hValueH += m_ImageH[hIdx];
5091 hCnt++;
5092 }
5093 }
5094 hValueL /= hCnt;
5095 hValueC /= hCnt;
5096 hValueH /= hCnt;
5097
5098 // fill mask with radius and value threshold
5099 fill4stack(FillMask, APointX, APointY, m_Width, m_Height,
5100 [&](uint16_t X, uint16_t Y) -> float {
5101
5102 float hResult = 1.0f;
5103 if (AUseMaxRadius) {
5104 hRad = ptSqr(std::abs((int32_t)X-APointX)) + ptSqr(std::abs((int32_t)Y-APointY));
5105 if (hRad < hRadiusOut) hResult = ptSqr((hRadiusOut - hRad)/hRadiusOut);
5106 else hResult = 0.0f;
5107 }
5108
5109 if (hResult == 0.0f) return hResult;
5110
5111 hIdx = (int32_t)Y*m_Width + X;
5112
5113 hDiff = std::abs((float)m_ImageC[hIdx] - hValueC) +
5114 std::abs((float)m_ImageL[hIdx] - hValueL)*hLumaWeight +
5115 std::abs((float)m_ImageH[hIdx] - hValueH)*hColorWeight*(float)m_ImageC[hIdx]/(float)0x1fff;
5116
5117 if (hDiff > hThresholdHalf) {
5118 if (hDiff < hThreshold) hResult = hResult*powf((hThreshold - hDiff)/hThresholdHalf, 4.0f);
5119 else hResult = 0.0f;
5120 }
5121
5122 return hResult;
5123 });
5124
5125 return FillMask;
5126 }
5127
5128 //==============================================================================
5129
MaskedColorAdjust(const int Ax,const int Ay,const float AThreshold,const float AChromaWeight,const int AMaxRadius,const bool AHasMaxRadius,const ptCurve * ACurve,const bool ASatAdaptive,float ASaturation,const float AHueShift)5130 ptImage *ptImage::MaskedColorAdjust(const int Ax,
5131 const int Ay,
5132 const float AThreshold,
5133 const float AChromaWeight,
5134 const int AMaxRadius,
5135 const bool AHasMaxRadius,
5136 const ptCurve *ACurve,
5137 const bool ASatAdaptive,
5138 float ASaturation,
5139 const float AHueShift)
5140 {
5141 assert (m_ColorSpace == ptSpace_LCH);
5142
5143 float *hMask = FillMask(Ax, Ay, AThreshold*2.0f, AChromaWeight, AMaxRadius, AHasMaxRadius);
5144 const bool hSatAdjust = ASaturation != 0.0f;
5145 const bool hHueAdjust = AHueShift != 0.0f;
5146 const float hHueShift = AHueShift * pt2PI;
5147
5148 // we enhance positive saturation values; negative values should give B&W with
5149 // -1.0, at least when not adaptive.
5150 if (ASaturation > 0) ASaturation *= 2.0f;
5151
5152 #pragma omp parallel for default(shared)
5153 for (uint32_t i=0; i< (uint32_t)m_Height*m_Width; i++) {
5154 if (hMask[i] > 0.0f) {
5155 m_ImageL[i] = m_ImageL[i] * (1.0f - hMask[i]) +
5156 ACurve->Curve[m_ImageL[i]] * hMask[i];
5157
5158 if (hSatAdjust) {
5159 if (ASatAdaptive)
5160 m_ImageC[i] = m_ImageC[i] * (1.0f + hMask[i] * ASaturation * (2.0f - m_ImageC[i]/(float)0x3FFF));
5161 else
5162 m_ImageC[i] = m_ImageC[i] * (1.0f + hMask[i] * ASaturation);
5163 }
5164
5165 if (hHueAdjust) {
5166 m_ImageH[i] = m_ImageH[i] + hHueShift * hMask[i];
5167 }
5168 }
5169 }
5170
5171 FREE(hMask);
5172 return this;
5173 }
5174
5175 ////////////////////////////////////////////////////////////////////////////////
5176 //
5177 // GetVignetteMask
5178 //
5179 ////////////////////////////////////////////////////////////////////////////////
5180
GetVignetteMask(const bool Inverted,const TVignetteShape Shape,const double InnerRadius,const double OuterRadius,const double Roundness,const double CenterX,const double CenterY,const double Softness)5181 pt::c_unique_ptr<float> ptImage::GetVignetteMask(
5182 const bool Inverted,
5183 const TVignetteShape Shape,
5184 const double InnerRadius,
5185 const double OuterRadius,
5186 const double Roundness,
5187 const double CenterX,
5188 const double CenterY,
5189 const double Softness)
5190 {
5191 float Radius = MIN(m_Width, m_Height)/2;
5192 float Exponent = static_cast<int>(Shape);
5193 bool Switch = OuterRadius < InnerRadius;
5194
5195 float OR = Radius*(Switch?InnerRadius:OuterRadius);
5196 float IR = Radius*(Switch?OuterRadius:InnerRadius);
5197 float Black = Inverted?1.0:0.0;
5198 float White = Inverted?0.0:1.0;
5199
5200 if (Switch) {
5201 Black = 1.0 - Black;
5202 White = 1.0 - White;
5203 }
5204 float ColorDiff = White - Black;
5205
5206 float CX = (1+CenterX)*m_Width/2;
5207 float CY = (1-CenterY)*m_Height/2;
5208
5209 float InversExponent = 1.0/ Exponent;
5210 float coordinate = 0;
5211 float Value = 0;
5212
5213 pt::c_unique_ptr<float> VignetteMask(static_cast<float*>(CALLOC(m_Width*m_Height,sizeof(float))));
5214 ptMemoryError(VignetteMask.get(),__FILE__,__LINE__);
5215
5216 float dist = 0;
5217 float Denom = 1/MAX((OR-IR),0.0001f);
5218 float Factor1 = 1/powf(2,Roundness);
5219 float Factor2 = 1/powf(2,-Roundness);
5220
5221 #pragma omp parallel for schedule(static) default(shared) firstprivate(dist, Value, coordinate)
5222 for (uint16_t Row=0; Row<m_Height; Row++) {
5223 int32_t Idx = Row*m_Width;
5224
5225 for (uint16_t Col=0; Col<m_Width; Col++) {
5226 dist = powf(powf(fabsf((float)Col-CX)*Factor1,Exponent)
5227 + powf(fabsf((float)Row-CY)*Factor2,Exponent),InversExponent);
5228
5229 if (dist <= IR) {
5230 VignetteMask.get()[Idx+Col] = Black;
5231 } else if (dist >= OR) {
5232 VignetteMask.get()[Idx+Col] = White;
5233 } else {
5234 coordinate = 1.0f-(OR-dist)*Denom;
5235 Value = (1.0f-powf(cosf(coordinate*ptPI/2.0f),50.0f*Softness))
5236 * powf(coordinate,0.07f*Softness);
5237 VignetteMask.get()[Idx+Col] = LIM(Value*ColorDiff+Black,0.0f,1.0f);
5238 }
5239 }
5240 }
5241
5242 return VignetteMask;
5243 }
5244
5245 ////////////////////////////////////////////////////////////////////////////////
5246 //
5247 // Gradual Blur
5248 //
5249 ////////////////////////////////////////////////////////////////////////////////
5250
GradualBlur(const TGradualBlurMode Mode,const double MaxRadius,const double LowerLevel,const double UpperLevel,const double Softness,const double Angle,const TVignetteShape Vignette,const double Roundness,const double CenterX,const double CenterY)5251 ptImage* ptImage::GradualBlur(const TGradualBlurMode Mode,
5252 const double MaxRadius,
5253 const double LowerLevel,
5254 const double UpperLevel,
5255 const double Softness,
5256 const double Angle,
5257 const TVignetteShape Vignette,
5258 const double Roundness,
5259 const double CenterX,
5260 const double CenterY)
5261 {
5262
5263 pt::c_unique_ptr<float> Mask;
5264
5265 if (Mode == TGradualBlurMode::Linear || Mode == TGradualBlurMode::LinearMask) {
5266 Mask = GetGradualMask(Angle, LowerLevel, UpperLevel, Softness);
5267 } else if (Mode == TGradualBlurMode::Vignette || Mode == TGradualBlurMode::VignetteMask) {
5268 Mask = GetVignetteMask(false, Vignette, LowerLevel, UpperLevel, Roundness, CenterX, CenterY, Softness);
5269 }
5270
5271 if ((Mode == TGradualBlurMode::LinearMask) || (Mode == TGradualBlurMode::VignetteMask)) {
5272 Overlay(m_Image, 1, Mask.get(), TOverlayMode::ShowMask);
5273 } else {
5274 Box((uint16_t)(ceil(MaxRadius)), Mask.get());
5275 }
5276
5277 return this;
5278 }
5279
5280 ////////////////////////////////////////////////////////////////////////////////
5281 //
5282 // Box
5283 //
5284 ////////////////////////////////////////////////////////////////////////////////
5285
Box(const uint16_t MaxRadius,float * Mask)5286 ptImage* ptImage::Box(const uint16_t MaxRadius, float* Mask) {
5287 ptImage* Source = new ptImage();
5288 Source->Set(this);
5289
5290 const uint16_t Height1 = m_Height - 1;
5291 const uint16_t Width1 = m_Width - 1;
5292 const uint16_t Height1_2 = 2*Height1;
5293 const uint16_t Width1_2 = 2*Width1;
5294
5295 // Construct the distance matrix
5296 const uint16_t DistSize = MaxRadius+1;
5297 float** Dist = (float**) CALLOC(DistSize,sizeof(float*));
5298 for(int16_t i = 0; i < DistSize; i++) Dist[i] = (float*) CALLOC(DistSize,sizeof(float));
5299
5300 // fill distance matrix
5301 for(int16_t i = 0; i <= MaxRadius; i++) {
5302 for(int16_t j = 0; j <= MaxRadius; j++) {
5303 Dist[i][j] = powf((float) i*i + (float) j*j, 0.5f);
5304 }
5305 }
5306
5307 float Sum0, Sum1, Sum2, Radius;
5308 int16_t i, j;
5309 int32_t Temp, NewRow, NewCol, Index;
5310 uint16_t Row, Col, Count, IntRadius, *PtrSource, *PtrTarget;
5311
5312 #pragma omp parallel for private(Sum0, Sum1, Sum2, i, j, Temp, Row, Col, Count, Radius, IntRadius, NewRow, NewCol, PtrSource, PtrTarget, Index)
5313 for (Row=0; Row<m_Height; Row++) {
5314 Temp = Row*m_Width;
5315 for (Col=0; Col<m_Width; Col++) {
5316 Sum0 = 0.0f;
5317 Sum1 = 0.0f;
5318 Sum2 = 0.0f;
5319 Count = 0;
5320
5321 Index = Temp + Col;
5322 Radius = MaxRadius * Mask[Index];
5323 IntRadius = ceil(Radius);
5324 if (IntRadius == 0) continue;
5325
5326 for(i = -IntRadius; i <= IntRadius; i++) {
5327 NewRow = Row+i;
5328 NewRow = NewRow < 0? -NewRow : NewRow > Height1? Height1_2-NewRow : NewRow ;
5329 NewRow *= m_Width;
5330 for(j = -IntRadius; j <= IntRadius; j++) {
5331 if (Dist[abs(i)][abs(j)] < Radius) {
5332 NewCol = Col+j;
5333 NewCol = NewCol < 0? -NewCol : NewCol > Width1? Width1_2-NewCol : NewCol ;
5334
5335 Count++;
5336 PtrSource = Source->m_Image[NewRow + NewCol];
5337 Sum0 += PtrSource[0];
5338 Sum1 += PtrSource[1];
5339 Sum2 += PtrSource[2];
5340 }
5341 }
5342 }
5343
5344 PtrTarget = m_Image[Index];
5345 PtrTarget[0] = Sum0 / Count;
5346 PtrTarget[1] = Sum1 / Count;
5347 PtrTarget[2] = Sum2 / Count;
5348 }
5349 }
5350
5351 for(uint16_t i = 0; i < DistSize; i++) FREE(Dist[i]);
5352 FREE(Dist);
5353
5354 delete Source;
5355
5356 return this;
5357 }
5358
5359 ////////////////////////////////////////////////////////////////////////////////
5360 //
5361 // ViewLAB
5362 //
5363 ////////////////////////////////////////////////////////////////////////////////
5364
ViewLab(const TViewLabChannel Channel)5365 ptImage* ptImage::ViewLab(const TViewLabChannel Channel) {
5366
5367 assert (m_ColorSpace == ptSpace_Lab);
5368
5369 const double WPH = 0x7fff; // WPH=WP/2
5370
5371 ptImage *ContrastLayer = new ptImage;
5372 ptCurve* MyContrastCurve = new ptCurve();
5373
5374 switch(Channel) {
5375 case TViewLabChannel::L:
5376 # pragma omp parallel for schedule(static)
5377 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5378 m_Image[i][1]=0x8080;
5379 m_Image[i][2]=0x8080;
5380 }
5381 break;
5382
5383 case TViewLabChannel::a:
5384 # pragma omp parallel for schedule(static)
5385 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5386 m_Image[i][0]=m_Image[i][1];
5387 m_Image[i][1]=0x8080;
5388 m_Image[i][2]=0x8080;
5389 }
5390 break;
5391
5392 case TViewLabChannel::b:
5393 # pragma omp parallel for schedule(static)
5394 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5395 m_Image[i][0]=m_Image[i][2];
5396 m_Image[i][1]=0x8080;
5397 m_Image[i][2]=0x8080;
5398 }
5399 break;
5400
5401 case TViewLabChannel::LStructure:
5402 ContrastLayer->Set(this);
5403 ContrastLayer->fastBilateralChannel(4, 0.2, 2, ChMask_L);
5404
5405 # pragma omp parallel for default(shared) schedule(static)
5406 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5407 m_Image[i][0] = CLIP((int32_t) ((WPH-(int32_t)ContrastLayer->m_Image[i][0])+m_Image[i][0]));
5408 }
5409 MyContrastCurve->setFromFunc(ptCurve::Sigmoidal,0.5,30);
5410 ApplyCurve(MyContrastCurve,1);
5411
5412 # pragma omp parallel for schedule(static)
5413 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5414 m_Image[i][1]=0x8080;
5415 m_Image[i][2]=0x8080;
5416 }
5417 break;
5418
5419 case TViewLabChannel::C:
5420 ContrastLayer->Set(this);
5421 ContrastLayer->LabToLch();
5422 # pragma omp parallel for schedule(static)
5423 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5424 m_Image[i][0]=CLIP((int32_t)(2.0f*ContrastLayer->m_ImageC[i]));
5425 m_Image[i][1]=0x8080;
5426 m_Image[i][2]=0x8080;
5427 }
5428 break;
5429
5430 case TViewLabChannel::H:
5431 ContrastLayer->Set(this);
5432 ContrastLayer->LabToLch();
5433 # pragma omp parallel for schedule(static)
5434 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5435 m_Image[i][0]=CLIP((int32_t)(ContrastLayer->m_ImageH[i]/pt2PI*(float)0xffff));
5436 m_Image[i][1]=0x8080;
5437 m_Image[i][2]=0x8080;
5438 }
5439 break;
5440
5441 default:
5442 assert("unexpected ViewLab channel mode");
5443 }
5444
5445 delete MyContrastCurve;
5446 delete ContrastLayer;
5447
5448 return this;
5449 }
5450
5451 ////////////////////////////////////////////////////////////////////////////////
5452 //
5453 // Special Preview
5454 //
5455 ////////////////////////////////////////////////////////////////////////////////
5456
SpecialPreview(const short Mode,const int Intent)5457 ptImage* ptImage::SpecialPreview(const short Mode, const int Intent) {
5458
5459 if (Mode == ptSpecialPreview_L ||
5460 Mode == ptSpecialPreview_A ||
5461 Mode == ptSpecialPreview_B) {
5462 int InColorSpace = m_ColorSpace;
5463 cmsHPROFILE InternalProfile = 0;
5464 if (InColorSpace != ptSpace_Profiled) {
5465 // linear case
5466 cmsToneCurve* Gamma = cmsBuildGamma(NULL, 1.0);
5467 cmsToneCurve* Gamma3[3];
5468 Gamma3[0] = Gamma3[1] = Gamma3[2] = Gamma;
5469
5470 cmsCIExyY DFromReference;
5471
5472 switch (m_ColorSpace) {
5473 case ptSpace_sRGB_D65 :
5474 case ptSpace_AdobeRGB_D65 :
5475 DFromReference = D65;
5476 break;
5477 case ptSpace_WideGamutRGB_D50 :
5478 case ptSpace_ProPhotoRGB_D50 :
5479 DFromReference = D50;
5480 break;
5481 default:
5482 assert(0);
5483 }
5484
5485 InternalProfile = cmsCreateRGBProfile(&DFromReference,
5486 (cmsCIExyYTRIPLE*)&RGBPrimaries[m_ColorSpace],
5487 Gamma3);
5488
5489 if (!InternalProfile) {
5490 ptLogError(ptError_Profile,"Could not open InternalProfile profile.");
5491 return this;
5492 }
5493
5494 cmsFreeToneCurve(Gamma);
5495 } else {
5496 // profiled case
5497 InternalProfile = PreviewColorProfile;
5498 }
5499
5500 cmsHPROFILE LabProfile = 0;
5501 LabProfile = cmsCreateLab4Profile(NULL);
5502 // to Lab
5503 cmsHTRANSFORM Transform;
5504 Transform = cmsCreateTransform(InternalProfile,
5505 TYPE_RGB_16,
5506 LabProfile,
5507 TYPE_Lab_16,
5508 Intent,
5509 cmsFLAGS_NOOPTIMIZE | cmsFLAGS_BLACKPOINTCOMPENSATION);
5510
5511 int32_t Size = m_Width*m_Height;
5512 int32_t Step = 100000;
5513 #pragma omp parallel for schedule(static)
5514 for (int32_t i = 0; i < Size; i+=Step) {
5515 int32_t Length = (i+Step)<Size ? Step : Size - i;
5516 uint16_t* Tile = &(m_Image[i][0]);
5517 cmsDoTransform(Transform,Tile,Tile,Length);
5518 }
5519 m_ColorSpace = ptSpace_Lab;
5520 // ViewLAB
5521 if (Mode == ptSpecialPreview_L) ViewLab(TViewLabChannel::L);
5522 if (Mode == ptSpecialPreview_A) ViewLab(TViewLabChannel::a);
5523 if (Mode == ptSpecialPreview_B) ViewLab(TViewLabChannel::b);
5524
5525 // to RGB
5526 Transform = cmsCreateTransform(LabProfile,
5527 TYPE_Lab_16,
5528 InternalProfile,
5529 TYPE_RGB_16,
5530 Intent,
5531 cmsFLAGS_NOOPTIMIZE | cmsFLAGS_BLACKPOINTCOMPENSATION);
5532
5533 #pragma omp parallel for schedule(static)
5534 for (int32_t i = 0; i < Size; i+=Step) {
5535 int32_t Length = (i+Step)<Size ? Step : Size - i;
5536 uint16_t* Tile = &(m_Image[i][0]);
5537 cmsDoTransform(Transform,Tile,Tile,Length);
5538 }
5539
5540 cmsDeleteTransform(Transform);
5541 cmsCloseProfile(LabProfile);
5542 if (InColorSpace != ptSpace_Profiled)
5543 cmsCloseProfile(InternalProfile);
5544 m_ColorSpace = InColorSpace;
5545 } else if (Mode==ptSpecialPreview_Structure) {
5546 #pragma omp parallel for default(shared) schedule(static)
5547 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5548 m_Image[i][0] = RGB_2_L(m_Image[i]);
5549 }
5550 const double WPH = 0x7fff; // WPH=WP/2
5551 ptImage *ContrastLayer = new ptImage;
5552 ContrastLayer->Set(this);
5553 ContrastLayer->fastBilateralChannel(4, 0.2, 2, ChMask_L);
5554 #pragma omp parallel for default(shared) schedule(static)
5555 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5556 m_Image[i][0] = CLIP((int32_t) ((WPH-(int32_t)ContrastLayer->m_Image[i][0]) + m_Image[i][0]));
5557 }
5558 ptCurve* MyContrastCurve = new ptCurve();
5559 MyContrastCurve->setFromFunc(ptCurve::Sigmoidal,0.5,30);
5560 ApplyCurve(MyContrastCurve,1);
5561 #pragma omp parallel for default(shared) schedule(static)
5562 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5563 m_Image[i][1] = m_Image[i][0];
5564 m_Image[i][2] = m_Image[i][0];
5565 }
5566 delete MyContrastCurve;
5567 delete ContrastLayer;
5568 } else if (Mode==ptSpecialPreview_Gradient) {
5569 #pragma omp parallel for default(shared) schedule(static)
5570 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5571 m_Image[i][0] = RGB_2_L(m_Image[i]);
5572 }
5573 ptCimgEdgeDetection(this,1);
5574 #pragma omp parallel for default(shared) schedule(static)
5575 for (uint32_t i=0; i<(uint32_t) m_Height*m_Width; i++) {
5576 m_Image[i][1] = m_Image[i][0];
5577 m_Image[i][2] = m_Image[i][0];
5578 }
5579 }
5580 return this;
5581 }
5582
5583 // -----------------------------------------------------------------------------
5584
mixChannels(const TChannelMatrix mixFactors)5585 ptImage* ptImage::mixChannels(const TChannelMatrix mixFactors) {
5586 assert (m_ColorSpace != ptSpace_Lab);
5587 double value[3];
5588
5589 # pragma omp parallel for default(shared) private(value)
5590 for (uint32_t px = 0; px < static_cast<uint32_t>(m_Height*m_Width); ++px) {
5591 for (int to = 0; to < 3; ++to) {
5592 value[to] = 0;
5593 for (int from = 0; from < 3; ++from) {
5594 value[to] += mixFactors[to][from] * m_Image[px][from];
5595 }
5596 }
5597 for (int to = 0; to < 3; ++to) {
5598 m_Image[px][to] = CLIP(static_cast<int32_t>(value[to]));
5599 }
5600 }
5601 return this;
5602 }
5603
5604 //==============================================================================
5605
Highlights(const float AHlRed,const float AHlGreen,const float AHlBlue)5606 ptImage *ptImage::Highlights(const float AHlRed, const float AHlGreen, const float AHlBlue) {
5607 TAnchorList hAnchors;
5608 for (int i=0; i<12; ++i) { // 12th anchor is a dummy
5609 hAnchors.push_back(TAnchor((double) i/33.0, (double) i/33.0));
5610 }
5611
5612 ptCurve HighlightsCurve(hAnchors, false); // no immediate curve calc
5613
5614 if (AHlRed < 0) {
5615 HighlightsCurve.setAnchor(11, TAnchor(1.0, 1.0+0.5*AHlRed));
5616 HighlightsCurve.calcCurve();
5617 ApplyCurve(&HighlightsCurve,1);
5618 } else if (AHlRed > 0) {
5619 HighlightsCurve.setAnchor(11, TAnchor(1.0-0.5*AHlRed, 1.0));
5620 HighlightsCurve.calcCurve();
5621 ApplyCurve(&HighlightsCurve,1);
5622 }
5623
5624 if (AHlGreen < 0) {
5625 HighlightsCurve.setAnchor(11, TAnchor(1.0, 1.0+0.5*AHlGreen));
5626 HighlightsCurve.calcCurve();
5627 ApplyCurve(&HighlightsCurve,2);
5628 } else if (AHlGreen > 0) {
5629 HighlightsCurve.setAnchor(11, TAnchor(1.0-0.5*AHlGreen, 1.0));
5630 HighlightsCurve.calcCurve();
5631 ApplyCurve(&HighlightsCurve,2);
5632 }
5633
5634 if (AHlBlue < 0) {
5635 HighlightsCurve.setAnchor(11, TAnchor(1.0, 1.0+0.5*AHlBlue));
5636 HighlightsCurve.calcCurve();
5637 ApplyCurve(&HighlightsCurve,4);
5638 } else if (AHlBlue > 0) {
5639 HighlightsCurve.setAnchor(11, TAnchor(1.0-0.5*AHlBlue, 1.0));
5640 HighlightsCurve.calcCurve();
5641 ApplyCurve(&HighlightsCurve,4);
5642 }
5643
5644 return this;
5645 }
5646
5647 //==============================================================================
5648
GammaTool(const float AStrength,const float ALinearity)5649 ptImage *ptImage::GammaTool(const float AStrength, const float ALinearity)
5650 {
5651 ptCurve RGBGammaCurve;
5652 RGBGammaCurve.setFromFunc(ptCurve::GammaTool, AStrength, ALinearity);
5653 ApplyCurve(&RGBGammaCurve, ChMask_RGB);
5654
5655 return this;
5656 }
5657
5658 ////////////////////////////////////////////////////////////////////////////////
5659 //
5660 // WriteAsPpm
5661 //
5662 ////////////////////////////////////////////////////////////////////////////////
5663
WriteAsPpm(const char * FileName,const short BitsPerColor)5664 short ptImage::WriteAsPpm(const char* FileName,
5665 const short BitsPerColor) {
5666
5667 assert ((8==BitsPerColor)||(16==BitsPerColor));
5668 assert (ptSpace_Lab != m_ColorSpace);
5669 assert (ptSpace_XYZ != m_ColorSpace);
5670
5671 FILE *OutputFile = fopen(FileName,"wb");
5672 if (!OutputFile) {
5673 ptLogError(ptError_FileOpen,FileName);
5674 return ptError_FileOpen;
5675 }
5676
5677 fprintf(OutputFile, "P%d\n%d %d\n%d\n", m_Colors/2+5, m_Width, m_Height,
5678 (1 << BitsPerColor)-1);
5679
5680 uint8_t* PpmRow = (uint8_t *) CALLOC(m_Width,m_Colors*BitsPerColor/8);
5681 ptMemoryError(PpmRow,__FILE__,__LINE__);
5682
5683 // Same buffer, interpreted as short though (16bit)
5684 uint16_t* PpmRow2= (uint16_t*) PpmRow;
5685
5686 for (uint16_t Row=0; Row<m_Height; Row++) {
5687 for (uint16_t Col=0; Col<m_Width; Col++) {
5688 if (8 == BitsPerColor) {
5689 for (short c=0;c<3;c++) {
5690 PpmRow [Col*m_Colors+c] = m_Image[Row*m_Width+Col][c] >>8;
5691 }
5692 } else {
5693 for (short c=0;c<3;c++) {
5694 PpmRow2[Col*m_Colors+c] = m_Image[Row*m_Width+Col][c];
5695 }
5696 }
5697 }
5698 if (16 == BitsPerColor && htons(0x55aa) != 0x55aa) {
5699 swab((char *)PpmRow,(char *)PpmRow,m_Width*m_Colors*2);
5700 }
5701 uint16_t w = fwrite(PpmRow,m_Colors*BitsPerColor/8,m_Width,OutputFile);
5702 assert(m_Width == w);
5703 }
5704
5705 FREE(PpmRow);
5706 FCLOSE(OutputFile);
5707 return 0;
5708 }
5709
5710 ////////////////////////////////////////////////////////////////////////////////
5711 //
5712 // WriteAsJpeg
5713 //
5714 ////////////////////////////////////////////////////////////////////////////////
5715 #ifndef NO_JPEG
WriteAsJpeg(const char * FileName,const short Quality,const uint8_t * ExifBuffer,const unsigned ExifBufferLen)5716 short ptImage::WriteAsJpeg(const char* FileName,
5717 const short Quality,
5718 const uint8_t* ExifBuffer,
5719 const unsigned ExifBufferLen ) {
5720
5721 assert (ptSpace_Lab != m_ColorSpace);
5722 assert (ptSpace_XYZ != m_ColorSpace);
5723
5724 FILE *OutputFile = fopen(FileName,"wb");
5725 if (!OutputFile) {
5726 ptLogError(ptError_FileOpen,FileName);
5727 return ptError_FileOpen;
5728 }
5729
5730 struct jpeg_compress_struct cinfo;
5731 struct jpeg_error_mgr jerr;
5732 cinfo.err = jpeg_std_error(&jerr);
5733 jpeg_create_compress(&cinfo);
5734 jpeg_stdio_dest(&cinfo,OutputFile);
5735 cinfo.image_width = m_Width;
5736 cinfo.image_height = m_Height;
5737 cinfo.input_components = m_Colors;
5738 cinfo.in_color_space = JCS_RGB; // m_colors = 1 grey ???
5739 jpeg_set_defaults(&cinfo);
5740 jpeg_set_quality(&cinfo,Quality,TRUE);
5741 jpeg_start_compress(&cinfo,TRUE);
5742 // TODO Embed exif and profile stuff here in between.
5743
5744 if (ExifBuffer) {
5745 if (ExifBufferLen > 65533) {
5746 ptLogWarning(ptWarning_Argument,
5747 "Exif buffer length %d is too long. Ignored",
5748 ExifBufferLen);
5749 } else {
5750 jpeg_write_marker(&cinfo,JPEG_APP0+1,ExifBuffer,ExifBufferLen);
5751 }
5752 }
5753
5754 uint8_t* PpmRow = (uint8_t *) CALLOC(m_Width,m_Colors);
5755 ptMemoryError(PpmRow,__FILE__,__LINE__);
5756
5757 for (uint16_t Row=0; Row<m_Height; Row++) {
5758 for (uint16_t Col=0; Col<m_Width; Col++) {
5759 for (short c=0;c<3;c++) {
5760 PpmRow [Col*m_Colors+c] = m_Image[Row*m_Width+Col][c] >>8;
5761 }
5762 }
5763 jpeg_write_scanlines(&cinfo,&PpmRow,1);
5764 }
5765
5766 jpeg_finish_compress(&cinfo);
5767 jpeg_destroy_compress(&cinfo);
5768
5769 FREE(PpmRow);
5770 FCLOSE(OutputFile);
5771 return 0;
5772 }
5773 #endif
5774
5775 ////////////////////////////////////////////////////////////////////////////////
5776 //
5777 // WaveletDenoise()
5778 // Basically from from a gimp plugin from dcraw to do a 'wavelet' denoise.
5779 //
5780 // Threshold : 0..0xffff
5781 // New implementation based on the wavelet denoise plugin for gimp 0.3.1
5782 // which is still based on dcraw...
5783 //
5784 ////////////////////////////////////////////////////////////////////////////////
5785
5786 // forward declaration of helper.
5787 static void hat_transform(float *temp, float *base,int st,int size,int sc);
5788
WaveletDenoise(const uint8_t ChannelMask,const double Threshold,const double low,const short WithMask,const double Sharpness,const double Anisotropy,const double Alpha,const double Sigma)5789 ptImage* ptImage::WaveletDenoise(const uint8_t ChannelMask,
5790 const double Threshold,
5791 const double low,
5792 const short WithMask,
5793 const double Sharpness,
5794 const double Anisotropy,
5795 const double Alpha,
5796 const double Sigma) {
5797
5798 assert(m_Colors==3);
5799 if (WithMask) assert(m_ColorSpace == ptSpace_Lab);
5800 uint32_t Size = m_Width*m_Height;
5801
5802 const float WP = 0xffff;
5803 float stdev[5];
5804 uint32_t samples[5];
5805
5806 // 3 : Image, lpass 0, lpass 1 + some tail for temporary processing.
5807 // tail is shifted for multithreading
5808 float *fImage = (float *) MALLOC((Size*3)*sizeof(*fImage));
5809 ptMemoryError(fImage,__FILE__,__LINE__);
5810
5811 for (short Channel=0; Channel<m_Colors; Channel++) {
5812
5813 // Channel supposed to handle ?
5814 if (! (ChannelMask & (1<<Channel))) continue;
5815
5816 // algorithm works between 0..1
5817 #pragma omp parallel for default(shared)
5818 for (uint32_t i=0; i<Size; i++) {
5819 fImage[i] = ToFloatTable[m_Image[i][Channel]];
5820 }
5821
5822 uint32_t lpass;
5823 uint32_t hpass = 0;
5824 for (uint16_t lev = 0; lev < 5; lev++) {
5825 lpass = Size*((lev & 1) + 1);
5826
5827 #pragma omp parallel default(shared)
5828 {
5829 #ifdef _OPENMP
5830 // We need a thread-private copy.
5831 float *TempTemp = (float *) MALLOC((m_Height+m_Width)*sizeof(*fImage));
5832 #else
5833 // Working in the trail (foreseen at CALLOC) done here for multithreading
5834 float *Temp = (float *) MALLOC((m_Height+m_Width)*sizeof(*fImage));
5835 #endif
5836 #pragma omp for
5837 for (uint16_t Row=0; Row<m_Height; Row++) {
5838 #ifdef _OPENMP
5839 hat_transform(TempTemp,fImage+hpass+Row*m_Width,1,m_Width,1<<lev);
5840 for (uint16_t Col = 0; Col<m_Width; Col++) {
5841 fImage[lpass+Row*m_Width+Col] = TempTemp[Col] * 0.25;
5842 }
5843 #else
5844 hat_transform(Temp,fImage+hpass+Row*m_Width,1,m_Width,1<<lev);
5845 for (uint16_t Col = 0; Col<m_Width; Col++) {
5846 fImage[lpass+Row*m_Width+Col] = Temp[Col] * 0.25;
5847 }
5848 #endif
5849 }
5850 #pragma omp for
5851 for (uint16_t Col=0; Col<m_Width; Col++) {
5852 #ifdef _OPENMP
5853 hat_transform(TempTemp,fImage+lpass+Col,m_Width,m_Height,1 << lev);
5854 for (uint16_t Row = 0; Row < m_Height; Row++) {
5855 fImage[lpass+Row*m_Width+Col] = TempTemp[Row] * 0.25;
5856 }
5857 #else
5858 hat_transform(Temp,fImage+lpass+Col,m_Width,m_Height,1 << lev);
5859 for (uint16_t Row = 0; Row < m_Height; Row++) {
5860 fImage[lpass+Row*m_Width+Col] = Temp[Row] * 0.25;
5861 }
5862 #endif
5863 }
5864 #ifdef _OPENMP
5865 // Free thread-private copy.
5866 free (TempTemp);
5867 #else
5868 // free
5869 free (Temp);
5870 #endif
5871
5872 } // End omp parallel zone.
5873
5874 float THold = 5.0 / (1 << 6) * exp (-2.6 * sqrt (lev + 1)) * 0.8002 / exp (-2.6);
5875
5876 /* initialize stdev values for all intensities */
5877 stdev[0] = stdev[1] = stdev[2] = stdev[3] = stdev[4] = 0.0;
5878 samples[0] = samples[1] = samples[2] = samples[3] = samples[4] = 0;
5879
5880 #pragma omp parallel default(shared)
5881 {
5882 /* calculate stdevs for all intensities */
5883 #ifdef _OPENMP
5884 // We need a thread-private copy.
5885 float Tempstdev[5] = {0.0,0.0,0.0,0.0,0.0};
5886 uint32_t Tempsamples[5] = {0,0,0,0,0};
5887 #endif
5888 #pragma omp for
5889 for (uint32_t i = 0; i < Size; i++) {
5890 fImage[hpass+i] -= fImage[lpass+i];
5891 if (fImage[hpass+i] < THold && fImage[hpass+i] > -THold) {
5892 #ifdef _OPENMP
5893 if (fImage[lpass+i] > 0.8) {
5894 Tempstdev[4] += fImage[hpass+i] * fImage[hpass+i];
5895 Tempsamples[4]++;
5896 } else if (fImage[lpass+i] > 0.6) {
5897 Tempstdev[3] += fImage[hpass+i] * fImage[hpass+i];
5898 Tempsamples[3]++;
5899 } else if (fImage[lpass+i] > 0.4) {
5900 Tempstdev[2] += fImage[hpass+i] * fImage[hpass+i];
5901 Tempsamples[2]++;
5902 } else if (fImage[lpass+i] > 0.2) {
5903 Tempstdev[1] += fImage[hpass+i] * fImage[hpass+i];
5904 Tempsamples[1]++;
5905 } else {
5906 Tempstdev[0] += fImage[hpass+i] * fImage[hpass+i];
5907 Tempsamples[0]++;
5908 }
5909 #else
5910 if (fImage[lpass+i] > 0.8) {
5911 stdev[4] += fImage[hpass+i] * fImage[hpass+i];
5912 samples[4]++;
5913 } else if (fImage[lpass+i] > 0.6) {
5914 stdev[3] += fImage[hpass+i] * fImage[hpass+i];
5915 samples[3]++;
5916 } else if (fImage[lpass+i] > 0.4) {
5917 stdev[2] += fImage[hpass+i] * fImage[hpass+i];
5918 samples[2]++;
5919 } else if (fImage[lpass+i] > 0.2) {
5920 stdev[1] += fImage[hpass+i] * fImage[hpass+i];
5921 samples[1]++;
5922 } else {
5923 stdev[0] += fImage[hpass+i] * fImage[hpass+i];
5924 samples[0]++;
5925 }
5926 #endif
5927 }
5928 }
5929 #ifdef _OPENMP
5930 #pragma omp critical
5931 for (int i = 0; i < 5; i++) {
5932 stdev[i] += Tempstdev[i];
5933 samples[i] += Tempsamples[i];
5934 }
5935 #endif
5936 } // End omp parallel zone.
5937 stdev[0] = sqrt (stdev[0] / (samples[0] + 1));
5938 stdev[1] = sqrt (stdev[1] / (samples[1] + 1));
5939 stdev[2] = sqrt (stdev[2] / (samples[2] + 1));
5940 stdev[3] = sqrt (stdev[3] / (samples[3] + 1));
5941 stdev[4] = sqrt (stdev[4] / (samples[4] + 1));
5942
5943 /* do thresholding */
5944 #pragma omp parallel for default(shared) private(THold)
5945 for (uint32_t i = 0; i < Size; i++) {
5946 if (fImage[lpass+i] > 0.8) {
5947 THold = Threshold * stdev[4];
5948 } else if (fImage[lpass+i] > 0.6) {
5949 THold = Threshold * stdev[3];
5950 } else if (fImage[lpass+i] > 0.4) {
5951 THold = Threshold * stdev[2];
5952 } else if (fImage[lpass+i] > 0.2) {
5953 THold = Threshold * stdev[1];
5954 } else {
5955 THold = Threshold * stdev[0];
5956 }
5957
5958 if (fImage[hpass+i] < -THold)
5959 fImage[hpass+i] += THold - THold * low;
5960 else if (fImage[hpass+i] > THold)
5961 fImage[hpass+i] -= THold - THold * low;
5962 else
5963 fImage[hpass+i] *= low;
5964
5965 if (hpass)
5966 fImage[0+i] += fImage[hpass+i];
5967 }
5968
5969 hpass = lpass;
5970 }
5971
5972 if (Channel==0 && WithMask && (Sharpness || (Anisotropy > 0.5))) {
5973 ptImage *MaskLayer = new ptImage;
5974 MaskLayer->Set(this);
5975 ptCimgEdgeTensors(MaskLayer,Sharpness,Anisotropy,Alpha,Sigma);
5976 #pragma omp parallel for default(shared)
5977 for (uint32_t i=0; i<Size; i++) {
5978 m_Image[i][Channel] = CLIP((int32_t)((fImage[i] + fImage[lpass+i])*MaskLayer->m_Image[i][0]+
5979 m_Image[i][Channel]*(1-(float)MaskLayer->m_Image[i][0]/(float)0xffff)));
5980 }
5981 delete MaskLayer;
5982 } else {
5983 #pragma omp parallel for default(shared)
5984 for (uint32_t i=0; i<Size; i++) {
5985 m_Image[i][Channel] = CLIP((int32_t)((fImage[i] + fImage[lpass+i])*WP));
5986 }
5987 }
5988 }
5989
5990 FREE(fImage);
5991
5992 return this;
5993 }
5994
hat_transform(float * temp,float * base,int st,int size,int sc)5995 static void hat_transform (float *temp, float *base, int st, int size, int sc) {
5996 int i;
5997 for (i = 0; i < sc; i++)
5998 temp[i] = 2 * base[st * i] + base[st * (sc - i)] + base[st * (i + sc)];
5999 for (; i + sc < size; i++)
6000 temp[i] = 2 * base[st * i] + base[st * (i - sc)] + base[st * (i + sc)];
6001 for (; i < size; i++)
6002 temp[i] = 2 * base[st * i] + base[st * (i - sc)]
6003 + base[st * (2 * size - 2 - (i + sc))];
6004 }
6005
createAmpAnchors(const double Amount,const double HaloControl)6006 TAnchorList ptImage::createAmpAnchors(const double Amount,
6007 const double HaloControl)
6008 {
6009 // TODO: mike: Anpassen der Verstärkungskurve an den Sigmoidalen Kontrast.
6010 // Der Berech mit Halocontrol soll linear sein, der andere sigmoidal.
6011 // Wichtig ist, dass bei kleinen Werten der Halocontrol im schwächeren Bereich
6012 // keine Verstärkung auftritt, da evtl die Ableitung der sigmoidalen Kurve
6013 // schwächer ist, als bei der linearen Kurve.
6014 // Anpassen bei allen Filtern, die dieses Verfahren benutzen.
6015
6016 const double t = (1.0 - Amount)/2;
6017 const double mHC = Amount*(1.0-fabs(HaloControl)); // m with HaloControl
6018 const double tHC = (1.0 - mHC)/2; // t with HaloControl
6019 const int Steps = 20;
6020
6021 TAnchorList hAnchors;
6022 for (int i = 0; i<= Steps; i++) {
6023 auto x = (float)i/(float)Steps;
6024 float YAnchor = 0;
6025 if (x < 0.5) {
6026 if (HaloControl > 0) YAnchor=mHC*x+tHC;
6027 else YAnchor=Amount*x+t;
6028 } else if (x > 0.5) {
6029 if (HaloControl < 0) YAnchor=mHC*x+tHC;
6030 else YAnchor=Amount*x+t;
6031 } else {
6032 YAnchor=0.5;
6033 }
6034 hAnchors.push_back(TAnchor(x, YAnchor));
6035 }
6036 return hAnchors;
6037 }
6038
setCurrentRGB(const short ARGB)6039 void ptImage::setCurrentRGB(const short ARGB)
6040 {
6041 CurrentRGBMode = ARGB;
6042 }
6043
getCurrentRGB()6044 short ptImage::getCurrentRGB()
6045 {
6046 return CurrentRGBMode;
6047 }
6048
6049 // -----------------------------------------------------------------------------
6050
6051 typedef Array_2D<float> image_type;
6052 extern float ToFloatTable[0x10000];
6053
6054 // From the theoretical part, the bilateral filter should blur more when values
6055 // are closer together. Since we use it with linear data, an additional gamma
6056 // correction could give better results.
6057
fastBilateralChannel(const float Sigma_s,const float Sigma_r,const int Iterations,const TChannelMask ChannelMask)6058 ptImage* ptImage::fastBilateralChannel(
6059 const float Sigma_s,
6060 const float Sigma_r,
6061 const int Iterations,
6062 const TChannelMask ChannelMask)
6063 {
6064 uint16_t Width = m_Width;
6065 uint16_t Height = m_Height;
6066 int32_t hIdx = 0;
6067
6068 image_type InImage(Width,Height);
6069 image_type FilteredImage(Width,Height);
6070
6071 for (short Channel = 0; Channel<3; Channel++) {
6072 // Is it a channel we are supposed to handle ?
6073 if (! (ChannelMask & (1<<Channel))) continue;
6074 #pragma omp parallel for default(shared) schedule(static) private(hIdx)
6075 for (uint16_t Row=0; Row<m_Height; Row++) {
6076 hIdx = Row*Width;
6077 for (uint16_t Col=0; Col<m_Width; Col++) {
6078 InImage(Col,Row) = ToFloatTable[m_Image[hIdx+Col][Channel]];
6079 }
6080 }
6081
6082 for (int i=0;i<Iterations;i++)
6083 Image_filter::fast_LBF(InImage,InImage,
6084 Sigma_s,Sigma_r,
6085 1,
6086 &FilteredImage,&FilteredImage);
6087
6088 #pragma omp parallel for default(shared) schedule(static) private(hIdx)
6089 for (uint16_t Row=0; Row<m_Height; Row++) {
6090 hIdx = Row*Width;
6091 for (uint16_t Col=0; Col<m_Width; Col++) {
6092 m_Image[hIdx+Col][Channel] = CLIP((int32_t)(FilteredImage(Col,Row)*0xffff));
6093 }
6094 }
6095 }
6096
6097 return this;
6098 }
6099