1 /*
2  */
3 
4 /*
5 
6     Copyright (C) 2014 Ferrero Andrea
7 
8     This program is free software: you can redistribute it and/or modify
9     it under the terms of the GNU General Public License as published by
10     the Free Software Foundation, either version 3 of the License, or
11     (at your option) any later version.
12 
13     This program is distributed in the hope that it will be useful,
14     but WITHOUT ANY WARRANTY; without even the implied warranty of
15     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16     GNU General Public License for more details.
17 
18     You should have received a copy of the GNU General Public License
19     along with this program. If not, see <http://www.gnu.org/licenses/>.
20 
21 
22  */
23 
24 /*
25 
26     These files are distributed with PhotoFlow - http://aferrero2707.github.io/PhotoFlow/
27 
28  */
29 
30 #include <string>
31 #include "photoflow.hh"
32 #include "format_info.hh"
33 #include "iccstore.hh"
34 
35 
36 cmsToneCurve* PF::Lstar_trc;
37 cmsToneCurve* PF::iLstar_trc;
38 
39 static cmsCIExyYTRIPLE rec2020_primaries = {
40     {0.7079, 0.2920, 1.0},
41     {0.1702, 0.7965, 1.0},
42     {0.1314, 0.0459, 1.0}
43 };
44 
45 static cmsCIExyYTRIPLE rec2020_primaries_prequantized = {
46     {0.708012540607, 0.291993664388, 1.0},
47     {0.169991652439, 0.797007778423, 1.0},
48     {0.130997824007, 0.045996550894, 1.0}
49 };
50 
51 static cmsCIExyYTRIPLE aces_cg_primaries =
52 {
53     {0.713, 0.293,  1.0},
54     {0.165, 0.830,  1.0},
55     {0.128, 0.044,  1.0}
56 };
57 
58 
59 /* ************************** WHITE POINTS ************************** */
60 
61 /* D50 WHITE POINTS */
62 
63 static cmsCIExyY d50_romm_spec= {0.3457, 0.3585, 1.0};
64 /* http://photo-lovers.org/pdf/color/romm.pdf */
65 
66 static cmsCIExyY d50_illuminant_specs = {0.345702915, 0.358538597, 1.0};
67 /* calculated from D50 illuminant XYZ values in ICC specs */
68 
69 
70 /* D65 WHITE POINTS */
71 
72 static cmsCIExyY  d65_srgb_adobe_specs = {0.3127, 0.3290, 1.0};
73 /* White point from the sRGB.icm and AdobeRGB1998 profile specs:
74  * http://www.adobe.com/digitalimag/pdfs/AdobeRGB1998.pdf
75  * 4.2.1 Reference Display White Point
76  * The chromaticity coordinates of white displayed on
77  * the reference color monitor shall be x=0.3127, y=0.3290.
78  * . . . [which] correspond to CIE Standard Illuminant D65.
79  *
80  * Wikipedia gives this same white point for SMPTE-C.
81  * This white point is also given in the sRGB color space specs.
82  * It's probably correct for most or all of the standard D65 profiles.
83  *
84  * The D65 white point values used in the LCMS virtual sRGB profile
85  * is slightly different than the D65 white point values given in the
86  * sRGB color space specs, so the LCMS virtual sRGB profile
87  * doesn't match sRGB profiles made using the values given in the
88  * sRGB color space specs.
89  *
90  * */
91 
92 
93 /* Various C and E WHITE POINTS */
94 
95 static cmsCIExyY c_astm  = {0.310060511, 0.316149551, 1.0};
96 /* see http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html */
97 static cmsCIExyY e_astm  = {0.333333333, 0.333333333, 1.0};
98 /* see http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html */
99 
100 static cmsCIExyY c_cie= {0.310, 0.316};
101 /* https://en.wikipedia.org/wiki/NTSC#Colorimetry */
102 static cmsCIExyY e_cie= {0.333, 0.333};
103 
104 static cmsCIExyY c_6774_robertson= {0.308548930, 0.324928102, 1.0};
105 /* see http://en.wikipedia.org/wiki/Standard_illuminant#White_points_of_standard_illuminants
106  * also see  http://www.brucelindbloom.com/index.html?Eqn_T_to_xy.html for the equations */
107 static cmsCIExyY e_5454_robertson= {0.333608970, 0.348572909, 1.0};
108 /* see http://en.wikipedia.org/wiki/Standard_illuminant#White_points_of_standard_illuminants
109  * also see http://www.brucelindbloom.com/index.html?Eqn_T_to_xy.html for the equations */
110 
111 
112 /* ACES white point, taken from
113  * Specification S-2014-004
114  * ACEScg – A Working Space for CGI Render and Compositing
115  */
116 static cmsCIExyY d60_aces= {0.32168, 0.33767, 1.0};
117 
118 /* Darktable code starts here
119  */
120 
121 #define generate_mat3inv_body(c_type, A, B)                                                                  \
122     static int mat3inv_##c_type(c_type dst[3][3], const c_type src[3][3])                                           \
123     {                                                                                                          \
124   \
125   const c_type det = A(1, 1) * (A(3, 3) * A(2, 2) - A(3, 2) * A(2, 3))                                     \
126   - A(2, 1) * (A(3, 3) * A(1, 2) - A(3, 2) * A(1, 3))                                   \
127   + A(3, 1) * (A(2, 3) * A(1, 2) - A(2, 2) * A(1, 3));                                  \
128   \
129   const c_type epsilon = 1e-7f;                                                                            \
130   if(fabs(det) < epsilon) return 1;                                                                        \
131   \
132   const c_type invDet = 1.0 / det;                                                                         \
133   \
134   B(1, 1) = invDet * (A(3, 3) * A(2, 2) - A(3, 2) * A(2, 3));                                              \
135   B(1, 2) = -invDet * (A(3, 3) * A(1, 2) - A(3, 2) * A(1, 3));                                             \
136   B(1, 3) = invDet * (A(2, 3) * A(1, 2) - A(2, 2) * A(1, 3));                                              \
137   \
138   B(2, 1) = -invDet * (A(3, 3) * A(2, 1) - A(3, 1) * A(2, 3));                                             \
139   B(2, 2) = invDet * (A(3, 3) * A(1, 1) - A(3, 1) * A(1, 3));                                              \
140   B(2, 3) = -invDet * (A(2, 3) * A(1, 1) - A(2, 1) * A(1, 3));                                             \
141   \
142   B(3, 1) = invDet * (A(3, 2) * A(2, 1) - A(3, 1) * A(2, 2));                                              \
143   B(3, 2) = -invDet * (A(3, 2) * A(1, 1) - A(3, 1) * A(1, 2));                                             \
144   B(3, 3) = invDet * (A(2, 2) * A(1, 1) - A(2, 1) * A(1, 2));                                              \
145   return 0;                                                                                                \
146     }
147 
148 #define A(y, x) src[y-1][x-1]
149 #define B(y, x) dst[y-1][x-1]
150 /** inverts the given 3x3 matrix */
generate_mat3inv_body(float,A,B)151 generate_mat3inv_body(float, A, B)
152 
153 static int mat3inv(float dst[3][3], const float src[3][3])
154 {
155   return mat3inv_float(dst, src);
156 }
157 /* Darktable code ends here
158  */
159 
160 
161 
162 
ICCProfile()163 PF::ICCProfile::ICCProfile()
164 {
165   has_colorants = false;
166   profile_data = NULL;
167   profile_size = 0;
168   profile_type = PROF_TYPE_CUSTOM;
169   trc_type = PF_TRC_UNKNOWN;
170   perceptual_trc = NULL;
171   perceptual_trc_inv = NULL;
172   parametric_trc = false;
173   gamut_boundary = NULL;
174   gamut_boundary_out = NULL;
175   to_lab = NULL;
176   from_lab = NULL;
177 }
178 
179 
~ICCProfile()180 PF::ICCProfile::~ICCProfile()
181 {
182 
183 }
184 
185 
set_profile(cmsHPROFILE p)186 void PF::ICCProfile::set_profile( cmsHPROFILE p )
187 {
188   profile = p;
189   cmsSaveProfileToMem( profile, NULL, &profile_size);
190   profile_data = malloc( profile_size );
191   cmsSaveProfileToMem( profile, profile_data, &profile_size);
192 
193   char tstr[1024];
194   cmsGetProfileInfoASCII(profile, cmsInfoDescription, "en", "US", tstr, 1024);
195 //#ifndef NDEBUG
196   //std::cout<<"ICCProfile::set_profile(): data="<<profile_data<<" size="<<profile_size<<"  name="<<tstr<<std::endl;
197 //#endif
198 
199   init_colorants();
200   init_trc();
201 }
202 
203 
init_colorants()204 void PF::ICCProfile::init_colorants()
205 {
206   has_colorants = false;
207   if( !profile ) return;
208   /* get the profile colorant information and fill in colorants */
209   cmsCIEXYZ *red            = (cmsCIEXYZ*)cmsReadTag(profile, cmsSigRedColorantTag);
210   if( !red ) return;
211   cmsCIEXYZ  red_colorant   = *red;
212   cmsCIEXYZ *green          = (cmsCIEXYZ*)cmsReadTag(profile, cmsSigGreenColorantTag);
213   if( !green ) return;
214   cmsCIEXYZ  green_colorant = *green;
215   cmsCIEXYZ *blue           = (cmsCIEXYZ*)cmsReadTag(profile, cmsSigBlueColorantTag);
216   if( !blue ) return;
217   cmsCIEXYZ  blue_colorant  = *blue;
218 
219   /* Get the Red channel XYZ values */
220   colorants[0]=red_colorant.X;
221   colorants[1]=red_colorant.Y;
222   colorants[2]=red_colorant.Z;
223 
224   /* Get the Green channel XYZ values */
225   colorants[3]=green_colorant.X;
226   colorants[4]=green_colorant.Y;
227   colorants[5]=green_colorant.Z;
228 
229   /* Get the Blue channel XYZ values */
230   colorants[6]=blue_colorant.X;
231   colorants[7]=blue_colorant.Y;
232   colorants[8]=blue_colorant.Z;
233 
234 #ifndef NDEBUG
235   //for( int i = 0; i < 9; i++ ) std::cout<<"colorants["<<i<<"]="<<colorants[i]<<std::endl;
236   std::cout<<"colorants:"<<std::endl;
237   for( int i = 0; i < 3; i++ ) std::cout<<colorants[i]<<" ";
238   std::cout<<std::endl;
239   for( int i = 3; i < 6; i++ ) std::cout<<colorants[i]<<" ";
240   std::cout<<std::endl;
241   for( int i = 6; i < 9; i++ ) std::cout<<colorants[i]<<" ";
242   std::cout<<std::endl;
243   //getchar();
244 
245   std::cout<<"RGB -> XYZ:"<<std::endl;
246   for( int i = 0; i < 3; i++ ) std::cout<<colorants[i*3]<<" ";
247   std::cout<<std::endl;
248   for( int i = 0; i < 3; i++ ) std::cout<<colorants[i*3+i]<<" ";
249   std::cout<<std::endl;
250   for( int i = 0; i < 3; i++ ) std::cout<<colorants[i*3+2]<<" ";
251   std::cout<<std::endl;
252   //getchar();
253 #endif
254 
255   Y_R = colorants[1];
256   Y_G = colorants[4];
257   Y_B = colorants[7];
258 
259   has_colorants = true;
260 
261   //float rgb2xyz[3][3];
262   double* pc = colorants;
263   for(int i=0;i<3;i++){
264     for(int j=0;j<3;j++){
265       rgb2xyz[j][i] = static_cast<float>( *pc );
266       pc += 1;
267     }
268   }
269 #ifndef NDEBUG
270   std::cout<<"RGB -> XYZ:"<<std::endl;
271   for( int i = 0; i < 3; i++ ) std::cout<<rgb2xyz[0][i]<<" ";
272   std::cout<<std::endl;
273   for( int i = 0; i < 3; i++ ) std::cout<<rgb2xyz[1][i]<<" ";
274   std::cout<<std::endl;
275   for( int i = 0; i < 3; i++ ) std::cout<<rgb2xyz[2][i]<<" ";
276   std::cout<<std::endl;
277 #endif
278 
279 
280   float D50_to_D65[3][3] = {
281       {0.955576656, -0.023039343, 0.063163668},
282       {-0.028289547, 1.009941621, 0.021007661},
283       {0.012298179, -0.020483004, 1.329909891}
284   };
285   // multiply rgb2xyz * D50_to_D65 to obtain the rgb2xyz100_D65 matrix
286   for(int i=0;i<3;i++){
287     for(int j=0;j<3;j++){
288       rgb2xyz100_D65[i][j]=0;
289       for(int k=0;k<3;k++){
290         rgb2xyz100_D65[i][j]=rgb2xyz100_D65[i][j]+(D50_to_D65[i][k] * rgb2xyz[k][j]);
291       }
292     }
293   }
294   mat3inv( xyz1002rgb_D65, rgb2xyz100_D65 );
295 #ifndef NDEBUG
296   std::cout<<"RGB -> XYZ_D65:"<<std::endl;
297   for( int i = 0; i < 3; i++ ) std::cout<<rgb2xyz100_D65[0][i]<<" ";
298   std::cout<<std::endl;
299   for( int i = 0; i < 3; i++ ) std::cout<<rgb2xyz100_D65[1][i]<<" ";
300   std::cout<<std::endl;
301   for( int i = 0; i < 3; i++ ) std::cout<<rgb2xyz100_D65[2][i]<<" ";
302   std::cout<<std::endl;
303 #endif
304 
305 
306 
307   cmsCIEXYZ *chad            = (cmsCIEXYZ*)cmsReadTag(profile, cmsSigChromaticAdaptationTag);
308 #ifndef NDEBUG
309   std::cout<<"chad tag: "<<chad<<std::endl;
310   if( chad ) {
311     for( int i = 0; i < 3; i++ ) {
312       std::cout<<chad[i].X<<" "<<chad[i].Y<<" "<<chad[i].Z<<std::endl;
313     }
314   }
315 #endif
316 
317 /*
318   std::string xyzprofname = PF::PhotoFlow::Instance().get_data_dir() + "/icc/XYZ-D50-Identity-elle-V4.icc";
319   cmsHPROFILE xyzprof = cmsOpenProfileFromFile( xyzprofname.c_str(), "r" );
320   cmsHTRANSFORM transform = cmsCreateTransform( profile,
321       TYPE_RGB_FLT,
322       xyzprof,
323       TYPE_RGB_FLT,
324       INTENT_RELATIVE_COLORIMETRIC,
325       cmsFLAGS_NOCACHE );
326 
327   float RGB[3], XYZ[3];
328   RGB[0] = 1; RGB[1] = 0; RGB[2] = 0;
329   cmsDoTransform(transform, RGB, XYZ, 1);
330   Y_R = XYZ[1];
331   RGB[0] = 0; RGB[1] = 1; RGB[2] = 0;
332   cmsDoTransform(transform, RGB, XYZ, 1);
333   Y_G = XYZ[1];
334   RGB[0] = 0; RGB[1] = 0; RGB[2] = 1;
335   cmsDoTransform(transform, RGB, XYZ, 1);
336   Y_B = XYZ[1];
337 */
338 }
339 
340 
init_trc()341 void PF::ICCProfile::init_trc()
342 {
343   if( !profile ) return;
344   /* get the profile colorant information and fill in colorants */
345   cmsToneCurve *red_trc   = (cmsToneCurve*)cmsReadTag(profile, cmsSigRedTRCTag);
346   //cmsToneCurve *green_trc = (cmsToneCurve*)cmsReadTag(profile, cmsSigGreenTRCTag);
347   //cmsToneCurve *blue_trc  = (cmsToneCurve*)cmsReadTag(profile, cmsSigBlueTRCTag);
348   //std::cout<<"ICCProfile::init_trc(): red_trc="<<red_trc<<std::endl;
349 
350   //if( !red_trc ) return;
351 
352   cmsBool is_linear = red_trc ? cmsIsToneCurveLinear(red_trc) : true;
353   cmsInt32Number tcpt = red_trc ? cmsGetToneCurveParametricType(red_trc) : 1;
354   parametric_trc = (tcpt>0) ? true : false;
355 
356   //std::cout<<"ICCProfile::init_trc(): is_linear="<<is_linear<<"  is_parametric="<<is_parametric()<<std::endl;
357 
358   if( is_linear ) {
359     /* LAB "L" (perceptually uniform) TRC */
360     cmsFloat64Number labl_parameters[5] =
361     { 3.0, 0.862076,  0.137924, 0.110703, 0.080002 };
362     cmsToneCurve *labl_parametic_curve =
363         cmsBuildParametricToneCurve(NULL, 4, labl_parameters);
364     cmsToneCurve *labl_parametic_curve_inv = cmsReverseToneCurve( labl_parametic_curve );
365     init_trc( labl_parametic_curve, labl_parametic_curve_inv );
366     if(get_trc_type() == PF::PF_TRC_UNKNOWN) set_trc_type( PF::PF_TRC_LINEAR );
367   } else {
368     cmsToneCurve* red_trc_inv = cmsReverseToneCurve( red_trc );
369     init_trc( red_trc, red_trc_inv );
370     if(get_trc_type() == PF::PF_TRC_UNKNOWN) set_trc_type( PF::PF_TRC_PERCEPTUAL );
371   }
372 }
373 
374 
init_trc(cmsToneCurve * trc,cmsToneCurve * trc_inv)375 void PF::ICCProfile::init_trc( cmsToneCurve* trc, cmsToneCurve* trc_inv )
376 {
377   perceptual_trc = trc;
378   perceptual_trc_inv = trc_inv;
379 
380   p2l_lut(256,0);
381   l2p_lut(256,0);
382 
383   for( int i = 0; i < 256; i++ ) {
384     cmsFloat32Number in=i, out;
385     in /= 255;
386     out = cmsEvalToneCurveFloat( perceptual_trc, in );
387     if( out > 1 ) out = 1;
388     if( out < 0 ) out = 0;
389     p2l_lut[i] = out;
390     //std::cout<<"ICCProfile::init_trc(): p2l_lut["<<i<<"] = "<<out<<std::endl;
391     out = cmsEvalToneCurveFloat( perceptual_trc_inv, in );
392     if( out > 1 ) out = 1;
393     if( out < 0 ) out = 0;
394     l2p_lut[i] = out;
395   }
396 
397   /*
398   for( int i = 0; i < 65536; i++ ) {
399     cmsFloat32Number in = i, out;
400     in /= 65535;
401     out = cmsEvalToneCurveFloat( perceptual_trc, in );
402     if( out > 1 ) out = 1;
403     if( out < 0 ) out = 0;
404     perceptual_trc_vec[i] = out;
405     out = cmsEvalToneCurveFloat( perceptual_trc_inv, in );
406     if( out > 1 ) out = 1;
407     if( out < 0 ) out = 0;
408     perceptual_trc_inv_vec[i] = out;
409   }
410   */
411 }
412 
413 
init_Lab_conversions(PF::ICCProfile * plab)414 void PF::ICCProfile::init_Lab_conversions( PF::ICCProfile* plab )
415 {
416   to_lab = new PF::ICCTransform;
417   from_lab = new PF::ICCTransform;
418   to_lab->init(this, plab, VIPS_FORMAT_FLOAT, INTENT_RELATIVE_COLORIMETRIC, false, 0);
419   from_lab->init(plab, this, VIPS_FORMAT_FLOAT, INTENT_RELATIVE_COLORIMETRIC, false, 0);
420 }
421 
422 
is_matrix()423 bool PF::ICCProfile::is_matrix()
424 {
425   if( !has_colorants ) return false;
426   if( !cmsIsMatrixShaper(get_profile()) ) return false;
427   if( cmsIsCLUT(get_profile(), INTENT_PERCEPTUAL, LCMS_USED_AS_INPUT) ) return false;
428   if( cmsIsCLUT(get_profile(), INTENT_PERCEPTUAL, LCMS_USED_AS_OUTPUT) ) return false;
429   return true;
430 }
431 
432 
is_rgb()433 bool PF::ICCProfile::is_rgb()
434 {
435   cmsColorSpaceSignature cs = cmsGetColorSpace( get_profile() );
436   return( cs == cmsSigRgbData );
437 }
438 
is_grayscale()439 bool PF::ICCProfile::is_grayscale()
440 {
441   cmsColorSpaceSignature cs = cmsGetColorSpace( get_profile() );
442   return( cs == cmsSigGrayData );
443 }
444 
is_lab()445 bool PF::ICCProfile::is_lab()
446 {
447   cmsColorSpaceSignature cs = cmsGetColorSpace( get_profile() );
448   return( cs == cmsSigLabData );
449 }
450 
is_cmyk()451 bool PF::ICCProfile::is_cmyk()
452 {
453   cmsColorSpaceSignature cs = cmsGetColorSpace( get_profile() );
454   return( cs == cmsSigCmykData );
455 }
456 
457 
get_profile()458 cmsHPROFILE PF::ICCProfile::get_profile()
459 {
460   return profile;
461 
462   cmsHPROFILE result;
463   cmsUInt32Number out_length;
464   cmsSaveProfileToMem( profile, NULL, &out_length);
465   void* buf = malloc( out_length );
466   cmsSaveProfileToMem( profile, buf, &out_length);
467 
468   result = cmsOpenProfileFromMem( buf, out_length );
469   //std::cout<<"ICCProfile::get_profile(): buf="<<buf<<std::endl;
470   //free( buf );
471   return result;
472 }
473 
474 
linear2perceptual(cmsFloat32Number val)475 cmsFloat32Number PF::ICCProfile::linear2perceptual( cmsFloat32Number val )
476 {
477   //return cmsEvalToneCurveFloat( perceptual_trc_inv, val );
478   //if( l2p_lut.getSize()==0 ) {
479     //std::cout<<"ICCProfile::linear2perceptual(): WARNING: empty LUT"<<std::endl;
480     //return 0; //cmsEvalToneCurveFloat( perceptual_trc_inv, val );
481   //}
482   cmsFloat32Number r = (val>1) ? cmsEvalToneCurveFloat( perceptual_trc_inv, val ) :
483       ( (val<0) ? cmsEvalToneCurveFloat( perceptual_trc_inv, val ) : l2p_lut[val*255] );
484   return r;
485 }
486 
perceptual2linear(cmsFloat32Number val)487 cmsFloat32Number PF::ICCProfile::perceptual2linear( cmsFloat32Number val )
488 {
489   //return cmsEvalToneCurveFloat( perceptual_trc, val );
490   //if( p2l_lut.getSize()==0 ) {
491     //std::cout<<"ICCProfile::perceptual2linear(): WARNING: empty LUT"<<std::endl;
492     //return 0; //cmsEvalToneCurveFloat( perceptual_trc, val );
493   //}
494   cmsFloat32Number r = (val>1) ? cmsEvalToneCurveFloat( perceptual_trc, val ) :
495       ( (val<0) ? cmsEvalToneCurveFloat( perceptual_trc, val ) : p2l_lut[val*255] );
496   return r;
497 }
498 
499 
500 
501 
get_lightness(const float & R,const float & G,const float & B)502 float PF::ICCProfile::get_lightness( const float& R, const float& G, const float& B )
503 {
504   if( !has_colorants ) return 0;
505   if( is_linear() ) {
506     return( Y_R*R + Y_G*G + Y_B*B );
507   } else {
508     float lR = perceptual2linear(R);
509     float lG = perceptual2linear(G);
510     float lB = perceptual2linear(B);
511     float L = Y_R*lR + Y_G*lG + Y_B*lB;
512     return linear2perceptual(L);
513   }
514 }
515 
get_lightness(float * RGBv,float * Lv,size_t size)516 void PF::ICCProfile::get_lightness( float* RGBv, float* Lv, size_t size )
517 {
518   if( !has_colorants ) {
519     for( size_t i = 0; i < size; i++ ) {
520       Lv[i] = 0;
521     }
522   }
523   if( is_linear() ) {
524     for( size_t i = 0, pos = 0; i < size; i++, pos += 3 ) {
525       Lv[i] = Y_R * RGBv[pos] + Y_G * RGBv[pos+1] + Y_B * RGBv[pos+2];
526     }
527   } else {
528     float lR, lG, lB, L;
529     for( size_t i = 0, pos = 0; i < size; i++, pos += 3 ) {
530       lR = perceptual2linear(RGBv[pos]);
531       lG = perceptual2linear(RGBv[pos+1]);
532       lB = perceptual2linear(RGBv[pos+2]);
533       L = Y_R*lR + Y_G*lG + Y_B*lB;
534       Lv[i] = linear2perceptual(L);
535     }
536   }
537 }
538 
539 
540 
to_Jzazbz(const float & R,const float & G,const float & B,float & Jz,float & az,float & bz)541 void PF::ICCProfile::to_Jzazbz( const float& R, const float& G, const float& B, float& Jz, float& az, float& bz )
542 {
543   float LAB[3];
544   float RGB[3] = {R, G, B};
545   to_lab->apply(RGB, LAB, 1);
546   Jz = LAB[0]/100; az = LAB[1]; bz = LAB[2];
547   return;
548 
549   float r = R; //(R >= 0) ? R : 0;
550   float g = G; //(G >= 0) ? G : 0;
551   float b = B; //(B >= 0) ? B : 0;
552   jabz::xyz100 xyz = {0, 0, 0};
553   xyz.x = (rgb2xyz100_D65[0][0]*r + rgb2xyz100_D65[0][1]*g + rgb2xyz100_D65[0][2]*b) * 100;
554   xyz.y = (rgb2xyz100_D65[1][0]*r + rgb2xyz100_D65[1][1]*g + rgb2xyz100_D65[1][2]*b) * 100;
555   xyz.z = (rgb2xyz100_D65[2][0]*r + rgb2xyz100_D65[2][1]*g + rgb2xyz100_D65[2][2]*b) * 100;
556   //std::cout<<"to_Jzazbz: xyz="<<xyz.x<<" "<<xyz.y<<" "<<xyz.z<<" "<<std::endl;
557 
558   jabz::jzazbz jab = jabz::jab_rgb::forth(xyz);
559   Jz = jab.jz;
560   az = jab.az;
561   bz = jab.bz;
562 }
563 
564 
565 
from_Jzazbz(const float & Jz,const float & az,const float & bz,float & R,float & G,float & B)566 void PF::ICCProfile::from_Jzazbz( const float& Jz, const float& az, const float& bz, float& R, float& G, float& B )
567 {
568   float LAB[3] = {Jz*100, az, bz};
569   float RGB[3];
570   from_lab->apply(LAB, RGB, 1);
571   R = RGB[0]; G = RGB[1]; B = RGB[2];
572   return;
573 
574   jabz::jzazbz jab = {Jz, az, bz};
575   jabz::xyz100 xyz = jabz::jab_rgb::back(jab);
576   //std::cout<<"from_Jzazbz: Jz="<<Jz<<" az="<<az<<"  bz="<<bz<<"    X="<<xyz.x<<" Y="<<xyz.y<<" Z="<<xyz.z<<std::endl;
577 
578   R = (xyz1002rgb_D65[0][0]*xyz.x + xyz1002rgb_D65[0][1]*xyz.y + xyz1002rgb_D65[0][2]*xyz.z) / 100;
579   G = (xyz1002rgb_D65[1][0]*xyz.x + xyz1002rgb_D65[1][1]*xyz.y + xyz1002rgb_D65[1][2]*xyz.z) / 100;
580   B = (xyz1002rgb_D65[2][0]*xyz.x + xyz1002rgb_D65[2][1]*xyz.y + xyz1002rgb_D65[2][2]*xyz.z) / 100;
581 }
582 
583 
init_gamut_mapping()584 void PF::ICCProfile::init_gamut_mapping()
585 {
586   if( gamut_boundary != NULL ) return;
587 
588   gamut_boundary = new float*[PF_GAMUT_MAP_NJZ+1];
589   for(int j = 0; j < PF_GAMUT_MAP_NJZ+1; j++) {
590     gamut_boundary[j] = new float[360];
591   }
592   gamut_Lid_Cmax = new float[360];
593 
594   float Jab[3], JCH[3], RGB[3] = {0.5, 0.2, 0.1};
595 
596   jabz::xyz100 xyz = {0.20654008*100, 0.12197225*0+100, 0.05136952*100};
597   jabz::jzazbz jab = jabz::jab_rgb::forth(xyz);
598   std::cout<<"init_gamut_mapping: xyz="<<xyz.x<<" "<<xyz.y<<" "<<xyz.z<<" "<<std::endl;
599   std::cout<<"init_gamut_mapping: Jab="<<jab.jz<<" "<<jab.az<<" "<<jab.bz<<" "<<std::endl;
600   xyz = jabz::jab_rgb::back(jab);
601   std::cout<<"init_gamut_mapping: xyz="<<xyz.x<<" "<<xyz.y<<" "<<xyz.z<<" "<<std::endl;
602   std::cout<<"init_gamut_mapping: Jab="<<jab.jz<<" "<<jab.az<<" "<<jab.bz<<" "<<std::endl;
603   //getchar(); return;
604 
605 
606   std::cout<<"init_gamut_mapping: RGB="<<RGB[0]<<" "<<RGB[1]<<" "<<RGB[2]<<" "<<std::endl;
607   to_Jzazbz( RGB[0], RGB[1], RGB[2], Jab[0], Jab[1], Jab[2] );
608   std::cout<<"init_gamut_mapping: Jab="<<Jab[0]<<" "<<Jab[1]<<" "<<Jab[2]<<" "<<std::endl;
609   std::cout<<"init_gamut_mapping: JCh="<<JCH[0]<<" "<<JCH[1]<<" "<<JCH[2]<<" "<<std::endl;
610   //getchar(); return;
611 
612   float vCmax[360] = { 0 };
613   float vLidmax[360] = { 0 };
614   // re-compute the gamut boundaries
615   int NJsteps = PF_GAMUT_MAP_NJZ;
616   float Jdelta = 1.0f / (NJsteps);
617   float J = Jdelta;
618   for(int h = 0; h < 360; h++) {
619     gamut_boundary[0][h] = gamut_boundary[NJsteps][h] = 0;
620   }
621   float delta = 1.0e-5;
622   for(int j = 1; j < NJsteps; j++) {
623     J = (Jdelta*j);
624     Jab[0] = J; Jab[1] = Jab[2] = 0;
625     from_Jzazbz( Jab[0], Jab[1], Jab[2], RGB[0], RGB[1], RGB[2] );
626     if( RGB[0] >= (1.0f-delta) && RGB[1] >= (1.0f-delta) && RGB[2] >= (1.0f-delta) ) {
627       for(int h = 0; h < 360; h++) {
628         gamut_boundary[j][h] = 0;
629       }
630       continue;
631     }
632 
633     for(int h = 0; h < 360; h++) {
634       float C = 0.1, Cmin = 0, Cmax = 0;
635       bool found = false;
636       gamut_boundary[j][h] = 0;
637       int iter = 0;
638       while( !found ) {
639         JCH[0] = J; JCH[1] = C; JCH[2] = h*M_PI/180.f;
640         PF::LCH2Lab(JCH, Jab, 1);
641         from_Jzazbz( Jab[0], Jab[1], Jab[2], RGB[0], RGB[1], RGB[2] );
642         //std::cout<<"Gamut: JCh="<<JCH[0]<<" "<<JCH[1]<<" "<<JCH[2]<<" "<<std::endl;
643         //std::cout<<"Gamut: Jab="<<Jab[0]<<" "<<Jab[1]<<" "<<Jab[2]<<" "<<std::endl;
644         //std::cout<<"Gamut: RGB="<<RGB[0]<<" "<<RGB[1]<<" "<<RGB[2]<<" "<<std::endl;
645         if( RGB[0] < (1.0f-delta) && RGB[1] < (1.0f-delta) && RGB[2] < (1.0f-delta) &&
646             RGB[0] > delta && RGB[1] > delta && RGB[2] > delta ) {
647           // we are still within gamut, increase C
648           //std::cout<<"Gamut: incresing C: "<<C<<" -> ";
649           if( C > Cmax ) C *= 2;
650           else {Cmin = C; C = (Cmax+Cmin)/2;}
651           //std::cout<<C<<std::endl;
652         } else if( RGB[0] > (1.0f+delta) || RGB[1] > (1.0f+delta) || RGB[2] > (1.0f+delta) ||
653             RGB[0] < -delta || RGB[1] < -delta || RGB[2] < -delta ) {
654           // at least one of the components is too much out-of-gamut, decrease C
655           //std::cout<<"Gamut: decreasing C: "<<C<<" -> ";
656           Cmax = C; C = (Cmax+Cmin)/2;
657           //std::cout<<C<<std::endl;
658         } else
659           found = true;
660 
661         if(found) {
662           gamut_boundary[j][h] = C;
663           if( C > vCmax[h] ) {
664             vCmax[h] = C;
665             vLidmax[h] = j;
666           }
667         }
668         iter++;
669         if(iter == 100) {break;}
670         //break;
671       }
672       std::cout<<"iter="<<iter<<"  J="<<J<<"  h="<<h<<"  C="<<C<<"  (min="<<Cmin<<" max="<<Cmax<<")  RGB="<<RGB[0]<<" "<<RGB[1]<<" "<<RGB[2]<<std::endl;
673       //break;
674     }
675     //std::cout<<"Y="<<Y<<" done."<<std::endl;
676     //break;
677   }
678 
679   for(int h = 0; h < 360; h++) {
680     gamut_Lid_Cmax[h] = vLidmax[h];
681   }
682   //getchar();
683 }
684 
685 
set_destination_gamut(ICCProfile * pout)686 void PF::ICCProfile::set_destination_gamut(ICCProfile* pout)
687 {
688   gamut_boundary_out = pout->get_gamut_boundary();
689 }
690 
691 
692 
gamut_mapping(float & R,float & G,float & B,float ** gamut_boundary_out,float * gamut_Lid_Cmax_out,float saturation)693 void PF::ICCProfile::gamut_mapping( float& R, float& G, float& B, float** gamut_boundary_out, float* gamut_Lid_Cmax_out, float saturation )
694 {
695   float Jab[3], JCH[3], RGB[3];
696   float Jz, az, bz, C, H;
697   // convert to Jzazbz in polar coordinates
698   //std::cout<<std::endl<<"gamut_mapping: RGB="<<R<<" "<<G<<" "<<B<<std::endl;
699   to_Jzazbz(R, G, B, Jab[0], Jab[1], Jab[2]);
700   PF::Lab2LCH(Jab, JCH, 1);
701   //std::cout<<"  J="<<Jab[0]<<"  a="<<Jab[1]<<"  b="<<Jab[2]
702   //    <<"  J="<<JCH[0]<<"  C="<<JCH[1]<<"  h="<<JCH[2]<<std::endl;
703 
704   if( !chroma_compression(JCH[0], JCH[1], JCH[2], gamut_boundary_out, gamut_Lid_Cmax_out, saturation) ) return;
705 
706   // re-calculate the RGB values
707   PF::LCH2Lab(JCH, Jab, 1);
708   from_Jzazbz( Jab[0], Jab[1], Jab[2], R, G, B );
709 }
710 
711 
chroma_compression(float & J,float & C,float & H,float ** gamut_boundary_out,float * gamut_Lid_Cmax_out,float saturation)712 bool PF::ICCProfile::chroma_compression( float& J, float& C, float& H, float** gamut_boundary_out, float* gamut_Lid_Cmax_out, float saturation )
713 {
714   //saturation = 0.5;
715   //if( std::isnan(J) ) {
716   //  J = 0;
717   //  C = 0;
718   //  return true;
719   //  getchar();
720   //}
721 
722   //if( J >= 1 ) {
723   //  C = 0;
724   //  return true;
725   //}
726 
727   if( J <= 0 ) {
728     C = 0;
729     return true;
730   }
731 
732   // get the index in the gamut mapping LUT
733   int j = static_cast<int>(J*PF_GAMUT_MAP_NJZ); if( j >= PF_GAMUT_MAP_NJZ ) j = PF_GAMUT_MAP_NJZ-1;
734   int h = static_cast<int>(H*180.0f/M_PI);
735   //std::cout<<"chroma_compression: J="<<J<<"  j="<<j<<"  h="<<h<<std::endl;
736 
737 
738   float Cmax2 = gamut_boundary_out[j][h];
739   int j2 = j;
740   if( saturation > 0 && j >= gamut_Lid_Cmax_out[h] ) {
741     // preserve saturation via luminance scaling
742     for(/*j2 = PF_GAMUT_MAP_NJZ-1*/; j2 >= 0; j2--) {
743       if( j2 < gamut_Lid_Cmax_out[h] ) break;
744       if( gamut_boundary_out[j2][h] > Cmax2 ) {
745         Cmax2 = gamut_boundary_out[j2][h];
746       }
747       if( Cmax2 > C ) break;
748     }
749     float J2 = ((float)j2)/PF_GAMUT_MAP_NJZ;
750     //std::cout<<"sat="<<saturation<<"  Jin="<<J<<"  J2="<<J2<<"  J="<<J-saturation * (J - J2)
751     //    <<"  j2="<<static_cast<int>(J*PF_GAMUT_MAP_NJZ)<<std::endl;
752     J -= saturation * (J - J2);
753     j2 = static_cast<int>(J*PF_GAMUT_MAP_NJZ);
754   }
755 
756   if( J >= 1 ) {
757     C = 0;
758     return true;
759   }
760 
761 
762   // compress all chroma values above 90% of Cmax
763   float Cmax = gamut_boundary_out[j2][h];
764   float C0 = Cmax * 0.9f;
765   float C1 = Cmax - C0;
766   //if( C>gamut_boundary[j][h]) std::cout<<"chroma_compression: j="<<j<<"  h="<<h<<"  C0="<<C0<<std::endl;
767   if( C < C0 ) return false;
768 
769   float Cout = ( C < C0 ) ? C : C0 + C1*( 1.0f - exp((C0-C)/C1) );
770   C = Cout;
771   return true;
772 }
773 
774 
775 
776 
equals_to(PF::ICCProfile * prof)777 bool PF::ICCProfile::equals_to( PF::ICCProfile* prof)
778 {
779   if( !prof ) return false;
780   if( prof->get_profile_size() == get_profile_size() &&
781       memcmp(prof->get_profile_data(), get_profile_data(), get_profile_size()) == 0 ) {
782     return true;
783   }
784   return false;
785 }
786 
787 
788 /*
789 PF::ICCProfileData* PF::get_icc_profile_data( VipsImage* img )
790 {
791   if( !img ) return NULL;
792   PF::ICCProfileData* data;
793   size_t data_length;
794   if( vips_image_get_blob( img, "pf-icc-profile-data",
795           (void**)(&data), &data_length ) ) {
796     std::cout<<"get_icc_profile_data(): cannot find ICC profile data"<<std::endl;
797     data = NULL;
798   }
799   if( data_length != sizeof(PF::ICCProfileData) ) {
800     std::cout<<"get_icc_profile_data(): wrong size of ICC profile data"<<std::endl;
801     data = NULL;
802   }
803   return data;
804 }
805 
806 
807 void PF::free_icc_profile_data( ICCProfileData* data )
808 {
809   cmsFreeToneCurve( data->perceptual_trc );
810   cmsFreeToneCurve( data->perceptual_trc_inv );
811 
812   delete data;
813 }
814 
815 
816 cmsFloat32Number PF::linear2perceptual( ICCProfileData* data, cmsFloat32Number val )
817 {
818   return cmsEvalToneCurveFloat( data->perceptual_trc_inv, val );
819 }
820 
821 
822 cmsFloat32Number PF::perceptual2linear( ICCProfileData* data, cmsFloat32Number val )
823 {
824   return cmsEvalToneCurveFloat( data->perceptual_trc, val );
825 }
826 */
827 
828 
829 //#define PF_TRC_LUT_SIZE 255
830 #define PF_TRC_LUT_SIZE 65535
831 
832 
init(ICCProfile * pin,ICCProfile * pout,VipsBandFormat band_fmt,cmsUInt32Number _intent,bool _bpc,float _adaptation_state)833 void PF::ICCTransform::init(ICCProfile* pin, ICCProfile* pout, VipsBandFormat band_fmt,
834     cmsUInt32Number _intent, bool _bpc, float _adaptation_state)
835 {
836   if( transform ) {
837     cmsDeleteTransform( transform );
838     transform = NULL;
839   }
840   trc_lut.reset();
841   itrc_lut.reset();
842 
843   in_profile = pin;
844   out_profile = pout;
845   if( !pin || !pout) return;
846 
847   bpc = _bpc;
848   intent = _intent;
849   adaptation_state = _adaptation_state;
850 
851   input_cs_type = cmsGetColorSpace( in_profile->get_profile() );
852   output_cs_type = cmsGetColorSpace( out_profile->get_profile() );
853 
854   transform = NULL;
855   is_rgb2rgb = false;
856 
857   //std::cout<<"ICCTransform::init() called"<<std::endl;
858 
859   bool is_matrix = in_profile->is_matrix() && out_profile->is_matrix();
860   bool is_parametric = (in_profile->is_linear() || in_profile->is_parametric()) &&
861       (out_profile->is_linear() || out_profile->is_parametric());
862   bool is_linear = in_profile->is_linear() && out_profile->is_linear();
863 
864   if( true && in_profile->is_rgb() && out_profile->is_rgb() &&
865       is_matrix && is_parametric && is_linear &&
866       !bpc && intent==INTENT_RELATIVE_COLORIMETRIC ) {
867     // fast path for linear RGB -> RGB matrix conversions
868     // get input profile colorants
869     double* in_colorants = in_profile->get_colorants();
870     float rgb2xyz[3][3];
871     for(int i=0;i<3;i++){
872       for(int j=0;j<3;j++){
873         rgb2xyz[j][i] = static_cast<float>( *in_colorants );
874         in_colorants += 1;
875       }
876     }
877     //printf("rgb2xyz:\n");
878     //printf("  %0.5f %0.5f %0.5f \n", rgb2xyz[0][0], rgb2xyz[0][1], rgb2xyz[0][2]);
879     //printf("  %0.5f %0.5f %0.5f \n", rgb2xyz[1][0], rgb2xyz[1][1], rgb2xyz[1][2]);
880     //printf("  %0.5f %0.5f %0.5f \n", rgb2xyz[2][0], rgb2xyz[2][1], rgb2xyz[2][2]);
881     // get output profile colorants, and invert them
882     double* out_colorants = out_profile->get_colorants();
883     float rgb2xyz2[3][3];
884     for(int i=0;i<3;i++){
885       for(int j=0;j<3;j++){
886         rgb2xyz2[j][i] = static_cast<float>( *out_colorants );
887         out_colorants += 1;
888       }
889     }
890     float xyz2rgb[3][3];
891     mat3inv( xyz2rgb, rgb2xyz2 );
892     //printf("rgb2xyz2:\n");
893     //printf("  %0.5f %0.5f %0.5f \n", rgb2xyz2[0][0], rgb2xyz2[0][1], rgb2xyz2[0][2]);
894     //printf("  %0.5f %0.5f %0.5f \n", rgb2xyz2[1][0], rgb2xyz2[1][1], rgb2xyz2[1][2]);
895     //printf("  %0.5f %0.5f %0.5f \n", rgb2xyz2[2][0], rgb2xyz2[2][1], rgb2xyz2[2][2]);
896     //printf("xyz2rgb:\n");
897     //printf("  %0.5f %0.5f %0.5f \n", xyz2rgb[0][0], xyz2rgb[0][1], xyz2rgb[0][2]);
898     //printf("  %0.5f %0.5f %0.5f \n", xyz2rgb[1][0], xyz2rgb[1][1], xyz2rgb[1][2]);
899     //printf("  %0.5f %0.5f %0.5f \n", xyz2rgb[2][0], xyz2rgb[2][1], xyz2rgb[2][2]);
900 
901     // multiply rgb2xyz * xyz2rgb to obtain the rgb2rgb matrix
902     for(int i=0;i<3;i++){
903       for(int j=0;j<3;j++){
904         rgb2rgb[i][j]=0;
905         for(int k=0;k<3;k++){
906           rgb2rgb[i][j]=rgb2rgb[i][j]+(xyz2rgb[i][k] * rgb2xyz[k][j]);
907         }
908       }
909     }
910     //printf("rgb2rgb:\n");
911     //printf("  %0.5f %0.5f %0.5f \n", rgb2rgb[0][0], rgb2rgb[0][1], rgb2rgb[0][2]);
912     //printf("  %0.5f %0.5f %0.5f \n", rgb2rgb[1][0], rgb2rgb[1][1], rgb2rgb[1][2]);
913     //printf("  %0.5f %0.5f %0.5f \n", rgb2rgb[2][0], rgb2rgb[2][1], rgb2rgb[2][2]);
914     is_rgb2rgb = true;
915 
916     if( !in_profile->is_linear() ) {
917       itrc_lut(PF_TRC_LUT_SIZE+1,0);
918       itrc_lut.setClip(0);
919       for(int i = 0; i < PF_TRC_LUT_SIZE+1; i++) {
920         cmsFloat32Number in = i, out;
921         in /= PF_TRC_LUT_SIZE;
922         out = cmsEvalToneCurveFloat( in_profile->get_p2l_trc(), in );
923         itrc_lut[i] = out;
924       }
925     }
926 
927     if( !out_profile->is_linear() ) {
928       trc_lut(PF_TRC_LUT_SIZE+1,0);
929       trc_lut.setClip(0);
930       for(int i = 0; i < PF_TRC_LUT_SIZE+1; i++) {
931         cmsFloat32Number in = i, out;
932         in /= PF_TRC_LUT_SIZE;
933         out = cmsEvalToneCurveFloat( out_profile->get_l2p_trc(), in );
934         trc_lut[i] = out;
935       }
936     }
937   } //else {
938     if( in_profile && out_profile && out_profile->get_profile() ) {
939       //std::cout<<"ICCTransform::init(): getting input profile format"<<std::endl;
940       cmsUInt32Number infmt = vips2lcms_pixel_format( band_fmt, in_profile->get_profile() );
941       //std::cout<<"ICCTransform::init(): getting output profile format"<<std::endl;
942       cmsUInt32Number outfmt = vips2lcms_pixel_format( band_fmt, out_profile->get_profile() );
943 
944       cmsUInt32Number flags = cmsFLAGS_NOOPTIMIZE | cmsFLAGS_NOCACHE;
945       //std::cout<<"ICCTransform::init(): bpc="<<bpc<<std::endl;
946       if( bpc ) flags |= cmsFLAGS_BLACKPOINTCOMPENSATION;
947       cmsFloat64Number old_state = cmsSetAdaptationState( adaptation_state );
948       transform = cmsCreateTransform( in_profile->get_profile(), infmt,
949           out_profile->get_profile(), outfmt, intent, flags );
950       cmsSetAdaptationState( old_state );
951       //std::cout<<"ICCTransform::init(): transform: "<<transform<<std::endl;
952       //std::cout<<"ICCTransform::init(): in_profile: "<<in_profile<<std::endl;
953       //std::cout<<"ICCTransform::init(): infmt: "<<infmt<<std::endl;
954       //std::cout<<"ICCTransform::init(): outfmt: "<<outfmt<<std::endl;
955     }
956   //}
957 }
958 
959 
apply(float * in,float * out,int n)960 void PF::ICCTransform::apply(float* in, float* out, int n)
961 {
962   if( is_rgb2rgb ) {
963     /* std::cout<<"ICCTransform::apply(): in="<<(void*)in<<"  out="<<(void*)out<<std::endl;
964     size_t addr = (size_t)in;
965     float faddr = (float)addr;
966     printf("    %f / 16 = %f\n",faddr,faddr/16);*/
967     float* in2 = in; float* out2 = out;
968     for(int i = 0; i < n; i++) {
969       float r = in2[0], g = in2[1], b = in2[2];
970       if(std::isnan(r)) r = 0;
971       if(std::isnan(g)) g = 0;
972       if(std::isnan(b)) b = 0;
973       float outr, outg, outb;
974       if(itrc_lut) {
975         r = (r>1) ? cmsEvalToneCurveFloat( in_profile->get_p2l_trc(), r ) :
976             ( (r<0) ? cmsEvalToneCurveFloat( in_profile->get_p2l_trc(), r ) : itrc_lut[r*PF_TRC_LUT_SIZE] );
977         g = (g>1) ? cmsEvalToneCurveFloat( in_profile->get_p2l_trc(), g ) :
978             ( (g<0) ? cmsEvalToneCurveFloat( in_profile->get_p2l_trc(), g ) : itrc_lut[g*PF_TRC_LUT_SIZE] );
979         b = (b>1) ? cmsEvalToneCurveFloat( in_profile->get_p2l_trc(), b ) :
980             ( (b<0) ? cmsEvalToneCurveFloat( in_profile->get_p2l_trc(), b ) : itrc_lut[b*PF_TRC_LUT_SIZE] );
981       }
982       outr = rgb2rgb[0][0]*r + rgb2rgb[0][1]*g + rgb2rgb[0][2]*b;
983       outg = rgb2rgb[1][0]*r + rgb2rgb[1][1]*g + rgb2rgb[1][2]*b;
984       outb = rgb2rgb[2][0]*r + rgb2rgb[2][1]*g + rgb2rgb[2][2]*b;
985       if(std::isnan(outr)) outr = 0;
986       if(std::isnan(outg)) outg = 0;
987       if(std::isnan(outb)) outb = 0;
988       if(trc_lut) {
989         outr = (outr>1) ? cmsEvalToneCurveFloat( out_profile->get_l2p_trc(), outr ) :
990             ( (outr<0) ? cmsEvalToneCurveFloat( out_profile->get_l2p_trc(), outr ) : trc_lut[outr*PF_TRC_LUT_SIZE] );
991         outg = (outg>1) ? cmsEvalToneCurveFloat( out_profile->get_l2p_trc(), outg ) :
992             ( (outg<0) ? cmsEvalToneCurveFloat( out_profile->get_l2p_trc(), outg ) : trc_lut[outg*PF_TRC_LUT_SIZE] );
993         outb = (outb>1) ? cmsEvalToneCurveFloat( out_profile->get_l2p_trc(), outb ) :
994             ( (outb<0) ? cmsEvalToneCurveFloat( out_profile->get_l2p_trc(), outb ) : trc_lut[outb*PF_TRC_LUT_SIZE] );
995       }
996       out2[0] = outr;
997       out2[1] = outg;
998       out2[2] = outb;
999       in2 += 3; out2 += 3;
1000     }
1001     return;
1002     //std::cout<<"out(1): "<<out[0]<<","<<out[1]<<","<<out[2]<<std::endl;
1003   }
1004 
1005   if( !transform ) return;
1006   cmsDoTransform( transform, in, out, n );
1007   //std::cout<<"out(2): "<<out[0]<<","<<out[1]<<","<<out[2]<<std::endl;
1008 }
1009 
1010 
1011 struct ICCProfileContainer
1012 {
1013   PF::ICCProfile* profile;
1014 };
1015 
1016 
set_icc_profile(VipsImage * img,PF::ICCProfile * prof)1017 void PF::set_icc_profile( VipsImage* img, PF::ICCProfile* prof )
1018 {
1019   if( !prof ) {
1020     if( vips_image_remove( img, "pf-icc-profile" ) == FALSE ) {
1021       std::cout<<"set_icc_profile: failed to remove \"pf-icc-profile\" metadata"<<std::endl;
1022     }
1023     return;
1024   }
1025 
1026   //std::cout<<"set_icc_profile: prof="<<prof<<std::endl;
1027   if( prof->get_profile() ) {
1028     char tstr[1024];
1029     cmsGetProfileInfoASCII(prof->get_profile(), cmsInfoDescription, "en", "US", tstr, 1024);
1030     //std::cout<<"set_icc_profile: profile name : "<<tstr<<std::endl;
1031   } else {
1032     std::cout<<"set_icc_profile: prof->get_profile() == NULL"<<std::endl;
1033   }
1034 
1035   ICCProfileContainer* pc = (ICCProfileContainer*)g_malloc( sizeof(ICCProfileContainer) );
1036   if( pc ) {
1037     pc->profile = prof;
1038     vips_image_set_blob( img, "pf-icc-profile",
1039         (VipsCallbackFn)g_free, pc, sizeof(ICCProfileContainer) );
1040 
1041     void* buf = g_malloc( prof->get_profile_size() );
1042     if( buf ) {
1043       //std::cout<<"PF::set_icc_profile(): VIPS_META_ICC blob set (size="<<prof->get_profile_size()<<")"<<std::endl;
1044       memcpy( buf, prof->get_profile_data(), prof->get_profile_size() );
1045       vips_image_set_blob( img, VIPS_META_ICC_NAME,
1046           (VipsCallbackFn) g_free, buf, prof->get_profile_size() );
1047     } else {
1048       std::cerr<<"PF::set_icc_profile(): cannot allocate "<<prof->get_profile_size()<<" bytes"<<std::endl;
1049     }
1050   } else {
1051     std::cerr<<"PF::set_icc_profile(): cannot allocate "<<prof->get_profile_size()<<" bytes"<<std::endl;
1052   }
1053 }
1054 
1055 
1056 
get_icc_profile(VipsImage * img)1057 PF::ICCProfile* PF::get_icc_profile( VipsImage* img )
1058 {
1059   if( !img ) return NULL;
1060   ICCProfileContainer* pc;
1061   size_t data_size;
1062   if( PF_VIPS_IMAGE_GET_BLOB( img, "pf-icc-profile", &pc, &data_size ) ) {
1063     //std::cout<<"get_icc_profile_data(): cannot find ICC profile data"<<std::endl;
1064     return( NULL );
1065   }
1066   if( data_size != sizeof(ICCProfileContainer) ) {
1067     //std::cout<<"get_icc_profile(): wrong size of ICC profile data"<<std::endl;
1068     return( NULL );
1069   }
1070   if( !pc ) return NULL;
1071   return(pc->profile);
1072 }
1073 
1074 
iccprofile_unref(void * prof)1075 void PF::iccprofile_unref( void* prof )
1076 {
1077   PF::ICCProfile* iccprof = (PF::ICCProfile*)( prof );
1078   if( !iccprof ) return;
1079   if( iccprof->get_ref_count() < 1 ) {
1080     std::cerr<<"WARNING!!! PF::iccprofile_unref("<<prof<<"): wrong refcount: "<<iccprof->get_ref_count()<<std::endl;
1081     return;
1082   }
1083 
1084   iccprof->unref();
1085   if( iccprof->get_ref_count() == 0 ) {
1086     std::cout<<"PF::iccprofile_unref("<<prof<<"): deleting profile"<<std::endl;
1087     delete iccprof;
1088   }
1089 }
1090 
1091 
1092 
1093 
ICCStore()1094 PF::ICCStore::ICCStore()
1095 {
1096   Lab_profile = new LabProfile( PF::PF_TRC_PERCEPTUAL );
1097   Lab_profile->ref(); profiles.push_back( Lab_profile );
1098 
1099   //std::cout<<"Initializing sRGB profiles"<<std::endl;
1100   srgb_profiles[0] = new sRGBProfile( PF::PF_TRC_STANDARD );
1101   srgb_profiles[1] = new sRGBProfile( PF::PF_TRC_PERCEPTUAL );
1102   srgb_profiles[2] = new sRGBProfile( PF::PF_TRC_LINEAR );
1103   srgb_profiles[3] = new sRGBProfile( PF::PF_TRC_sRGB );
1104   srgb_profiles[4] = new sRGBProfile( PF::PF_TRC_GAMMA_22 );
1105   srgb_profiles[5] = new sRGBProfile( PF::PF_TRC_GAMMA_18 );
1106   srgb_profiles[0]->ref(); profiles.push_back( srgb_profiles[0] );
1107   srgb_profiles[1]->ref(); profiles.push_back( srgb_profiles[1] );
1108   srgb_profiles[2]->ref(); profiles.push_back( srgb_profiles[2] );
1109   srgb_profiles[3]->ref(); profiles.push_back( srgb_profiles[3] );
1110   srgb_profiles[4]->ref(); profiles.push_back( srgb_profiles[4] );
1111   srgb_profiles[5]->ref(); profiles.push_back( srgb_profiles[5] );
1112   //std::cout<<"sRGB profiles initialization finished"<<std::endl
1113   //    <<"====================================="<<std::endl;
1114 
1115   srgb_d50_profiles[0] = new sRGBProfileD50( PF::PF_TRC_STANDARD );
1116   srgb_d50_profiles[1] = new sRGBProfileD50( PF::PF_TRC_PERCEPTUAL );
1117   srgb_d50_profiles[2] = new sRGBProfileD50( PF::PF_TRC_LINEAR );
1118   srgb_d50_profiles[0]->ref(); profiles.push_back( srgb_d50_profiles[0] );
1119   srgb_d50_profiles[1]->ref(); profiles.push_back( srgb_d50_profiles[1] );
1120   srgb_d50_profiles[2]->ref(); profiles.push_back( srgb_d50_profiles[2] );
1121 
1122   rec2020_profiles[0] = new Rec2020Profile( PF::PF_TRC_STANDARD );
1123   rec2020_profiles[1] = new Rec2020Profile( PF::PF_TRC_PERCEPTUAL );
1124   rec2020_profiles[2] = new Rec2020Profile( PF::PF_TRC_LINEAR );
1125   rec2020_profiles[0]->ref(); profiles.push_back( rec2020_profiles[0] );
1126   rec2020_profiles[1]->ref(); profiles.push_back( rec2020_profiles[1] );
1127   rec2020_profiles[2]->ref(); profiles.push_back( rec2020_profiles[2] );
1128 
1129   aces_profiles[0] = new ACESProfile( PF::PF_TRC_STANDARD );
1130   aces_profiles[1] = new ACESProfile( PF::PF_TRC_PERCEPTUAL );
1131   aces_profiles[2] = new ACESProfile( PF::PF_TRC_LINEAR );
1132   aces_profiles[0]->ref(); profiles.push_back( aces_profiles[0] );
1133   aces_profiles[1]->ref(); profiles.push_back( aces_profiles[1] );
1134   aces_profiles[2]->ref(); profiles.push_back( aces_profiles[2] );
1135 
1136   acescg_profiles[0] = get_profile(PF::PhotoFlow::Instance().get_data_dir() + "/icc/ACEScg-elle-V4-srgbtrc.icc");
1137   acescg_profiles[1] = get_profile(PF::PhotoFlow::Instance().get_data_dir() + "/icc/ACEScg-elle-V4-labl.icc");
1138   acescg_profiles[2] = get_profile(PF::PhotoFlow::Instance().get_data_dir() + "/icc/ACEScg-elle-V4-g10.icc");
1139 
1140   adobe_profiles[0] = get_profile(PF::PhotoFlow::Instance().get_data_dir() + "/icc/ClayRGB-elle-V4-g22.icc");
1141   adobe_profiles[1] = get_profile(PF::PhotoFlow::Instance().get_data_dir() + "/icc/ClayRGB-elle-V4-labl.icc");
1142   adobe_profiles[2] = get_profile(PF::PhotoFlow::Instance().get_data_dir() + "/icc/ClayRGB-elle-V4-g10.icc");
1143 
1144   //prophoto_profiles[0] = get_profile(PF::PhotoFlow::Instance().get_data_dir() + "/icc/LargeRGB-elle-V4-g18.icc");
1145   //prophoto_profiles[1] = get_profile(PF::PhotoFlow::Instance().get_data_dir() + "/icc/LargeRGB-elle-V4-labl.icc");
1146   //prophoto_profiles[2] = get_profile(PF::PhotoFlow::Instance().get_data_dir() + "/icc/LargeRGB-elle-V4-g10.icc");
1147   prophoto_profiles[0] = new ProPhotoProfile( PF::PF_TRC_STANDARD );
1148   prophoto_profiles[1] = new ProPhotoProfile( PF::PF_TRC_PERCEPTUAL );
1149   prophoto_profiles[2] = new ProPhotoProfile( PF::PF_TRC_LINEAR );
1150   prophoto_profiles[0]->ref(); profiles.push_back( prophoto_profiles[0] );
1151   prophoto_profiles[1]->ref(); profiles.push_back( prophoto_profiles[1] );
1152   prophoto_profiles[2]->ref(); profiles.push_back( prophoto_profiles[2] );
1153 
1154   XYZ_profile = new XYZProfile( PF::PF_TRC_LINEAR );
1155   XYZ_profile->ref(); profiles.push_back( XYZ_profile );
1156 
1157   system_monitor_profile = NULL;
1158 
1159   for(unsigned int pi = 0; pi < profiles.size(); pi++) {
1160     if( profiles[pi]->is_lab() ) continue;
1161     profiles[pi]->init_Lab_conversions(Lab_profile);
1162   }
1163 
1164   /*
1165   std::string wprofname = PF::PhotoFlow::Instance().get_data_dir() + "/icc/Rec2020-elle-V4-g10.icc";
1166   profile = cmsOpenProfileFromFile( wprofname.c_str(), "r" );
1167 
1168   d50_wp_primaries = rec2020_primaries_prequantized;
1169    */
1170 
1171   /* LAB "L" (perceptually uniform) TRC */
1172 
1173   cmsFloat64Number labl_parameters[5] =
1174   { 3.0, 0.862076,  0.137924, 0.110703, 0.080002 };
1175   Lstar_trc = cmsBuildParametricToneCurve(NULL, 4, labl_parameters);
1176   iLstar_trc = cmsBuildParametricToneCurve(NULL, 4, labl_parameters);
1177   iLstar_trc = cmsReverseToneCurve( iLstar_trc );
1178 
1179   /*
1180   for( int i = 0; i < 65536; i++ ) {
1181     cmsFloat32Number in = i, out;
1182     in /= 65535;
1183     out = cmsEvalToneCurveFloat( perceptual_trc, in )*65535;
1184     if( out > 65535 ) out = 65535;
1185     if( out < 0 ) out = 0;
1186     perceptual_trc_vec[i] = (int)out;
1187     out = cmsEvalToneCurveFloat( perceptual_trc_inv, in )*65535;
1188     if( out > 65535 ) out = 65535;
1189     if( out < 0 ) out = 0;
1190     perceptual_trc_inv_vec[i] = (int)out;
1191   }
1192    */
1193 
1194   defaultMonitorProfile = "";
1195 
1196 #ifdef WIN32
1197   /*
1198   // Get current main monitor. Could be fine tuned to get the current windows monitor (multi monitor setup),
1199   // but problem is that we live in RTEngine with no GUI window to query around
1200   HDC hDC = GetDC(NULL);
1201 
1202   if (hDC != NULL) {
1203       if (SetICMMode(hDC, ICM_ON)) {
1204           char profileName[MAX_PATH + 1];
1205           DWORD profileLength = MAX_PATH;
1206 
1207           if (GetICMProfileA(hDC, &profileLength, profileName)) {
1208               defaultMonitorProfile = Glib::ustring(profileName);
1209           }
1210 
1211           // might fail if e.g. the monitor has no profile
1212       }
1213 
1214       ReleaseDC(NULL, hDC);
1215   }
1216   */
1217 #else
1218 // TODO: Add other OS specific code here
1219 #endif
1220 }
1221 
1222 
get_profile(PF::profile_type_t ptype,PF::TRC_type trc_type)1223 PF::ICCProfile* PF::ICCStore::get_profile( PF::profile_type_t ptype, PF::TRC_type trc_type)
1224 {
1225   //std::cout<<"ICCStore::get_profile("<<ptype<<", "<<trc_type<<")"<<std::endl;
1226   switch( ptype ) {
1227   case PF::PROF_TYPE_sRGB: return srgb_profiles[trc_type];
1228   case PF::PROF_TYPE_sRGB_D50: return srgb_d50_profiles[trc_type];
1229   case PF::PROF_TYPE_REC2020: return rec2020_profiles[trc_type];
1230   case PF::PROF_TYPE_ACES: return aces_profiles[trc_type];
1231   case PF::PROF_TYPE_ACEScg: return acescg_profiles[trc_type];
1232   case PF::PROF_TYPE_ADOBE: return adobe_profiles[trc_type];
1233   case PF::PROF_TYPE_PROPHOTO: return prophoto_profiles[trc_type];
1234   case PF::PROF_TYPE_LAB: return Lab_profile;
1235   case PF::PROF_TYPE_XYZ:
1236     //std::cout<<"  PROF_TYPE_XYZ: XYZ_profile="<<(void*)XYZ_profile<<std::endl;
1237     return XYZ_profile;
1238   default: return NULL;
1239   }
1240 }
1241 
1242 
get_profile(Glib::ustring pname)1243 PF::ICCProfile* PF::ICCStore::get_profile( Glib::ustring pname )
1244 {
1245   // First we loop over all the opened profiles to see if there is
1246   // already one associated to the same file
1247   for(unsigned int pi = 0; pi < profiles.size(); pi++ ) {
1248     if( profiles[pi]->get_file_name() == pname ) {
1249       return profiles[pi];
1250     }
1251   }
1252 
1253   // Next, we check if an equivalent profile exists, even if not directly opened from disk
1254   cmsHPROFILE temp_profile = cmsOpenProfileFromFile( pname.c_str(), "r" );
1255   if( temp_profile ) {
1256     cmsUInt32Number psize;
1257     cmsSaveProfileToMem( temp_profile, NULL, &psize);
1258     if( psize > 0 ) {
1259       void* pdata = malloc( psize );
1260       if( pdata ) {
1261         cmsSaveProfileToMem( temp_profile, pdata, &psize);
1262 
1263         // We loop over all the opened profiles to see if there is one with matching data
1264         for(unsigned int pi = 0; pi < profiles.size(); pi++ ) {
1265           if( profiles[pi]->get_profile_size() == psize &&
1266               memcmp(profiles[pi]->get_profile_data(), pdata, psize) == 0 ) {
1267             return profiles[pi];
1268           }
1269         }
1270 
1271         free( pdata );
1272       }
1273     }
1274 
1275     //std::cout<<"ICCStore::get_profile(): loading profile from \""<<pname<<"\""<<std::endl;
1276     PF::ICCProfile* new_profile = new PF::ICCProfile();
1277     new_profile->set_profile( temp_profile );
1278     new_profile->set_file_name( pname );
1279 
1280     profiles.push_back( new_profile );
1281     return new_profile;
1282   }
1283 
1284   return NULL;
1285 }
1286 
1287 
get_profile(void * pdata,cmsUInt32Number psize)1288 PF::ICCProfile* PF::ICCStore::get_profile( void* pdata, cmsUInt32Number psize )
1289 {
1290   if( !pdata || psize==0 ) return NULL;
1291   // We loop over all the opened profiles to see if there is one with matching data
1292   for(unsigned int pi = 0; pi < profiles.size(); pi++ ) {
1293     if( profiles[pi]->get_profile_size() == psize &&
1294         memcmp(profiles[pi]->get_profile_data(), pdata, psize) == 0 ) {
1295       return profiles[pi];
1296     }
1297   }
1298 
1299   void* buf = malloc( psize );
1300   if( !buf ) return NULL;
1301   memcpy( buf, pdata, psize );
1302 
1303   cmsHPROFILE temp_profile = cmsOpenProfileFromMem( buf, psize );
1304   free( buf );
1305 
1306   //std::cout<<"ICCStore::get_profile( void* pdata, cmsUInt32Number psize ): temp_profile="<<temp_profile<<std::endl;
1307 
1308   PF::ICCProfile* new_profile = new PF::ICCProfile();
1309   new_profile->set_profile( temp_profile );
1310 
1311   profiles.push_back( new_profile );
1312   return new_profile;
1313 }
1314 
1315 
get_profile(cmsHPROFILE prof)1316 PF::ICCProfile* PF::ICCStore::get_profile( cmsHPROFILE prof )
1317 {
1318   if( !prof ) return NULL;
1319   cmsUInt32Number psize;
1320   cmsSaveProfileToMem( prof, NULL, &psize);
1321   if( psize > 0 ) {
1322     void* pdata = g_malloc( psize );
1323     if( pdata ) {
1324       cmsSaveProfileToMem( prof, pdata, &psize);
1325 
1326       // We loop over all the opened profiles to see if there is one with matching data
1327       for(unsigned int pi = 0; pi < profiles.size(); pi++ ) {
1328         if( profiles[pi]->get_profile_size() == psize &&
1329             memcmp(profiles[pi]->get_profile_data(), pdata, psize) == 0 ) {
1330           g_free( pdata );
1331           cmsCloseProfile( prof );
1332           return profiles[pi];
1333         }
1334       }
1335 
1336       cmsHPROFILE temp_profile = cmsOpenProfileFromMem( pdata, psize );
1337       PF::ICCProfile* new_profile = new PF::ICCProfile();
1338       new_profile->set_profile( prof );
1339 
1340       profiles.push_back( new_profile );
1341       return new_profile;
1342     } else {
1343       return NULL;
1344     }
1345   } else {
1346     return NULL;
1347   }
1348 
1349   return NULL;
1350 }
1351 
1352 
set_system_monitor_profile(cmsHPROFILE icc_profile)1353 void PF::ICCStore::set_system_monitor_profile( cmsHPROFILE icc_profile )
1354 {
1355   PF::ICCProfile* profile = get_profile( icc_profile );
1356   system_monitor_profile = profile;
1357 }
1358 
1359 
1360 PF::ICCStore* PF::ICCStore::instance = NULL;
1361 
Instance()1362 PF::ICCStore& PF::ICCStore::Instance()
1363 {
1364   if(!PF::ICCStore::instance)
1365     PF::ICCStore::instance = new PF::ICCStore();
1366   return( *instance );
1367 };
1368