libs/lcms2/src/cmswtpnt.c

   1 //---------------------------------------------------------------------------------
   2 //
   3 //  Little Color Management System
   4 //  Copyright (c) 1998-2023 Marti Maria Saguer
   5 //
   6 // Permission is hereby granted, free of charge, to any person obtaining
   7 // a copy of this software and associated documentation files (the "Software"),
   8 // to deal in the Software without restriction, including without limitation
   9 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  10 // and/or sell copies of the Software, and to permit persons to whom the Software
  11 // is furnished to do so, subject to the following conditions:
  12 //
  13 // The above copyright notice and this permission notice shall be included in
  14 // all copies or substantial portions of the Software.
  15 //
  16 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
  18 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  20 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  21 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  22 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  23 //
  24 //---------------------------------------------------------------------------------
  25 //
  26
  27 #include "lcms2_internal.h"
  28
  29
  30 // D50 - Widely used
  31 const cmsCIEXYZ* CMSEXPORT cmsD50_XYZ(void)
  32 {
  33     static cmsCIEXYZ D50XYZ = {cmsD50X, cmsD50Y, cmsD50Z};
  34
  35     return &D50XYZ;
  36 }
  37
  38 const cmsCIExyY* CMSEXPORT cmsD50_xyY(void)
  39 {
  40     static cmsCIExyY D50xyY;
  41
  42     cmsXYZ2xyY(&D50xyY, cmsD50_XYZ());
  43
  44     return &D50xyY;
  45 }
  46
  47 // Obtains WhitePoint from Temperature
  48 cmsBool  CMSEXPORT cmsWhitePointFromTemp(cmsCIExyY* WhitePoint, cmsFloat64Number TempK)
  49 {
  50     cmsFloat64Number x, y;
  51     cmsFloat64Number T, T2, T3;
  52     // cmsFloat64Number M1, M2;
  53
  54     _cmsAssert(WhitePoint != NULL);
  55
  56     T = TempK;
  57     T2 = T*T;            // Square
  58     T3 = T2*T;           // Cube
  59
  60     // For correlated color temperature (T) between 4000K and 7000K:
  61
  62     if (T >= 4000. && T <= 7000.)
  63     {
  64         x = -4.6070*(1E9/T3) + 2.9678*(1E6/T2) + 0.09911*(1E3/T) + 0.244063;
  65     }
  66     else
  67         // or for correlated color temperature (T) between 7000K and 25000K:
  68
  69         if (T > 7000.0 && T <= 25000.0)
  70         {
  71             x = -2.0064*(1E9/T3) + 1.9018*(1E6/T2) + 0.24748*(1E3/T) + 0.237040;
  72         }
  73         else {
  74             cmsSignalError(0, cmsERROR_RANGE, "cmsWhitePointFromTemp: invalid temp");
  75             return FALSE;
  76         }
  77
  78     // Obtain y(x)
  79     y = -3.000*(x*x) + 2.870*x - 0.275;
  80
  81     // wave factors (not used, but here for futures extensions)
  82
  83     // M1 = (-1.3515 - 1.7703*x + 5.9114 *y)/(0.0241 + 0.2562*x - 0.7341*y);
  84     // M2 = (0.0300 - 31.4424*x + 30.0717*y)/(0.0241 + 0.2562*x - 0.7341*y);
  85
  86     WhitePoint -> x = x;
  87     WhitePoint -> y = y;
  88     WhitePoint -> Y = 1.0;
  89
  90     return TRUE;
  91 }
  92
  93
  94
  95 typedef struct {
  96
  97     cmsFloat64Number mirek;  // temp (in microreciprocal kelvin)
  98     cmsFloat64Number ut;     // u coord of intersection w/ blackbody locus
  99     cmsFloat64Number vt;     // v coord of intersection w/ blackbody locus
 100     cmsFloat64Number tt;     // slope of ISOTEMPERATURE. line
 101
 102     } ISOTEMPERATURE;
 103
 104 static const ISOTEMPERATURE isotempdata[] = {
 105 //  {Mirek, Ut,       Vt,      Tt      }
 106     {0,     0.18006,  0.26352,  -0.24341},
 107     {10,    0.18066,  0.26589,  -0.25479},
 108     {20,    0.18133,  0.26846,  -0.26876},
 109     {30,    0.18208,  0.27119,  -0.28539},
 110     {40,    0.18293,  0.27407,  -0.30470},
 111     {50,    0.18388,  0.27709,  -0.32675},
 112     {60,    0.18494,  0.28021,  -0.35156},
 113     {70,    0.18611,  0.28342,  -0.37915},
 114     {80,    0.18740,  0.28668,  -0.40955},
 115     {90,    0.18880,  0.28997,  -0.44278},
 116     {100,   0.19032,  0.29326,  -0.47888},
 117     {125,   0.19462,  0.30141,  -0.58204},
 118     {150,   0.19962,  0.30921,  -0.70471},
 119     {175,   0.20525,  0.31647,  -0.84901},
 120     {200,   0.21142,  0.32312,  -1.0182 },
 121     {225,   0.21807,  0.32909,  -1.2168 },
 122     {250,   0.22511,  0.33439,  -1.4512 },
 123     {275,   0.23247,  0.33904,  -1.7298 },
 124     {300,   0.24010,  0.34308,  -2.0637 },
 125     {325,   0.24702,  0.34655,  -2.4681 },
 126     {350,   0.25591,  0.34951,  -2.9641 },
 127     {375,   0.26400,  0.35200,  -3.5814 },
 128     {400,   0.27218,  0.35407,  -4.3633 },
 129     {425,   0.28039,  0.35577,  -5.3762 },
 130     {450,   0.28863,  0.35714,  -6.7262 },
 131     {475,   0.29685,  0.35823,  -8.5955 },
 132     {500,   0.30505,  0.35907,  -11.324 },
 133     {525,   0.31320,  0.35968,  -15.628 },
 134     {550,   0.32129,  0.36011,  -23.325 },
 135     {575,   0.32931,  0.36038,  -40.770 },
 136     {600,   0.33724,  0.36051,  -116.45  }
 137 };
 138
 139 #define NISO sizeof(isotempdata)/sizeof(ISOTEMPERATURE)
 140
 141
 142 // Robertson's method
 143 cmsBool  CMSEXPORT cmsTempFromWhitePoint(cmsFloat64Number* TempK, const cmsCIExyY* WhitePoint)
 144 {
 145     cmsUInt32Number j;
 146     cmsFloat64Number us,vs;
 147     cmsFloat64Number uj,vj,tj,di,dj,mi,mj;
 148     cmsFloat64Number xs, ys;
 149
 150     _cmsAssert(WhitePoint != NULL);
 151     _cmsAssert(TempK != NULL);
 152
 153     di = mi = 0;
 154     xs = WhitePoint -> x;
 155     ys = WhitePoint -> y;
 156
 157     // convert (x,y) to CIE 1960 (u,WhitePoint)
 158
 159     us = (2*xs) / (-xs + 6*ys + 1.5);
 160     vs = (3*ys) / (-xs + 6*ys + 1.5);
 161
 162
 163     for (j=0; j < NISO; j++) {
 164
 165         uj = isotempdata[j].ut;
 166         vj = isotempdata[j].vt;
 167         tj = isotempdata[j].tt;
 168         mj = isotempdata[j].mirek;
 169
 170         dj = ((vs - vj) - tj * (us - uj)) / sqrt(1.0 + tj * tj);
 171
 172         if ((j != 0) && (di/dj < 0.0)) {
 173
 174             // Found a match
 175             *TempK = 1000000.0 / (mi + (di / (di - dj)) * (mj - mi));
 176             return TRUE;
 177         }
 178
 179         di = dj;
 180         mi = mj;
 181     }
 182
 183     // Not found
 184     return FALSE;
 185 }
 186
 187
 188 // Compute chromatic adaptation matrix using Chad as cone matrix
 189
 190 static
 191 cmsBool ComputeChromaticAdaptation(cmsMAT3* Conversion,
 192                                 const cmsCIEXYZ* SourceWhitePoint,
 193                                 const cmsCIEXYZ* DestWhitePoint,
 194                                 const cmsMAT3* Chad)
 195
 196 {
 197
 198     cmsMAT3 Chad_Inv;
 199     cmsVEC3 ConeSourceXYZ, ConeSourceRGB;
 200     cmsVEC3 ConeDestXYZ, ConeDestRGB;
 201     cmsMAT3 Cone, Tmp;
 202
 203
 204     Tmp = *Chad;
 205     if (!_cmsMAT3inverse(&Tmp, &Chad_Inv)) return FALSE;
 206
 207     _cmsVEC3init(&ConeSourceXYZ, SourceWhitePoint -> X,
 208                              SourceWhitePoint -> Y,
 209                              SourceWhitePoint -> Z);
 210
 211     _cmsVEC3init(&ConeDestXYZ,   DestWhitePoint -> X,
 212                              DestWhitePoint -> Y,
 213                              DestWhitePoint -> Z);
 214
 215     _cmsMAT3eval(&ConeSourceRGB, Chad, &ConeSourceXYZ);
 216     _cmsMAT3eval(&ConeDestRGB,   Chad, &ConeDestXYZ);
 217
 218     if ((fabs(ConeSourceRGB.n[0]) < MATRIX_DET_TOLERANCE) ||
 219         (fabs(ConeSourceRGB.n[1]) < MATRIX_DET_TOLERANCE) ||
 220         (fabs(ConeSourceRGB.n[2]) < MATRIX_DET_TOLERANCE)) return FALSE;
 221
 222     // Build matrix
 223     _cmsVEC3init(&Cone.v[0], ConeDestRGB.n[0]/ConeSourceRGB.n[0],    0.0,  0.0);
 224     _cmsVEC3init(&Cone.v[1], 0.0,   ConeDestRGB.n[1]/ConeSourceRGB.n[1],   0.0);
 225     _cmsVEC3init(&Cone.v[2], 0.0,   0.0,   ConeDestRGB.n[2]/ConeSourceRGB.n[2]);
 226
 227     // Normalize
 228     _cmsMAT3per(&Tmp, &Cone, Chad);
 229     _cmsMAT3per(Conversion, &Chad_Inv, &Tmp);
 230
 231     return TRUE;
 232 }
 233
 234 // Returns the final chrmatic adaptation from illuminant FromIll to Illuminant ToIll
 235 // The cone matrix can be specified in ConeMatrix. If NULL, Bradford is assumed
 236 cmsBool  _cmsAdaptationMatrix(cmsMAT3* r, const cmsMAT3* ConeMatrix, const cmsCIEXYZ* FromIll, const cmsCIEXYZ* ToIll)
 237 {
 238     cmsMAT3 LamRigg   = {{ // Bradford matrix
 239         {{  0.8951,  0.2664, -0.1614 }},
 240         {{ -0.7502,  1.7135,  0.0367 }},
 241         {{  0.0389, -0.0685,  1.0296 }}
 242     }};
 243
 244     if (ConeMatrix == NULL)
 245         ConeMatrix = &LamRigg;
 246
 247     return ComputeChromaticAdaptation(r, FromIll, ToIll, ConeMatrix);
 248 }
 249
 250 // Same as anterior, but assuming D50 destination. White point is given in xyY
 251 static
 252 cmsBool _cmsAdaptMatrixToD50(cmsMAT3* r, const cmsCIExyY* SourceWhitePt)
 253 {
 254     cmsCIEXYZ Dn;
 255     cmsMAT3 Bradford;
 256     cmsMAT3 Tmp;
 257
 258     cmsxyY2XYZ(&Dn, SourceWhitePt);
 259
 260     if (!_cmsAdaptationMatrix(&Bradford, NULL, &Dn, cmsD50_XYZ())) return FALSE;
 261
 262     Tmp = *r;
 263     _cmsMAT3per(r, &Bradford, &Tmp);
 264
 265     return TRUE;
 266 }
 267
 268 // Build a White point, primary chromas transfer matrix from RGB to CIE XYZ
 269 // This is just an approximation, I am not handling all the non-linear
 270 // aspects of the RGB to XYZ process, and assuming that the gamma correction
 271 // has transitive property in the transformation chain.
 272 //
 273 // the algorithm:
 274 //
 275 //            - First I build the absolute conversion matrix using
 276 //              primaries in XYZ. This matrix is next inverted
 277 //            - Then I eval the source white point across this matrix
 278 //              obtaining the coefficients of the transformation
 279 //            - Then, I apply these coefficients to the original matrix
 280 //
 281 cmsBool _cmsBuildRGB2XYZtransferMatrix(cmsMAT3* r, const cmsCIExyY* WhitePt, const cmsCIExyYTRIPLE* Primrs)
 282 {
 283     cmsVEC3 WhitePoint, Coef;
 284     cmsMAT3 Result, Primaries;
 285     cmsFloat64Number xn, yn;
 286     cmsFloat64Number xr, yr;
 287     cmsFloat64Number xg, yg;
 288     cmsFloat64Number xb, yb;
 289
 290     xn = WhitePt -> x;
 291     yn = WhitePt -> y;
 292     xr = Primrs -> Red.x;
 293     yr = Primrs -> Red.y;
 294     xg = Primrs -> Green.x;
 295     yg = Primrs -> Green.y;
 296     xb = Primrs -> Blue.x;
 297     yb = Primrs -> Blue.y;
 298
 299     // Build Primaries matrix
 300     _cmsVEC3init(&Primaries.v[0], xr,        xg,         xb);
 301     _cmsVEC3init(&Primaries.v[1], yr,        yg,         yb);
 302     _cmsVEC3init(&Primaries.v[2], (1-xr-yr), (1-xg-yg),  (1-xb-yb));
 303
 304
 305     // Result = Primaries ^ (-1) inverse matrix
 306     if (!_cmsMAT3inverse(&Primaries, &Result))
 307         return FALSE;
 308
 309
 310     _cmsVEC3init(&WhitePoint, xn/yn, 1.0, (1.0-xn-yn)/yn);
 311
 312     // Across inverse primaries ...
 313     _cmsMAT3eval(&Coef, &Result, &WhitePoint);
 314
 315     // Give us the Coefs, then I build transformation matrix
 316     _cmsVEC3init(&r -> v[0], Coef.n[VX]*xr,          Coef.n[VY]*xg,          Coef.n[VZ]*xb);
 317     _cmsVEC3init(&r -> v[1], Coef.n[VX]*yr,          Coef.n[VY]*yg,          Coef.n[VZ]*yb);
 318     _cmsVEC3init(&r -> v[2], Coef.n[VX]*(1.0-xr-yr), Coef.n[VY]*(1.0-xg-yg), Coef.n[VZ]*(1.0-xb-yb));
 319
 320
 321     return _cmsAdaptMatrixToD50(r, WhitePt);
 322
 323 }
 324
 325
 326 // Adapts a color to a given illuminant. Original color is expected to have
 327 // a SourceWhitePt white point.
 328 cmsBool CMSEXPORT cmsAdaptToIlluminant(cmsCIEXYZ* Result,
 329                                        const cmsCIEXYZ* SourceWhitePt,
 330                                        const cmsCIEXYZ* Illuminant,
 331                                        const cmsCIEXYZ* Value)
 332 {
 333     cmsMAT3 Bradford;
 334     cmsVEC3 In, Out;
 335
 336     _cmsAssert(Result != NULL);
 337     _cmsAssert(SourceWhitePt != NULL);
 338     _cmsAssert(Illuminant != NULL);
 339     _cmsAssert(Value != NULL);
 340
 341     if (!_cmsAdaptationMatrix(&Bradford, NULL, SourceWhitePt, Illuminant)) return FALSE;
 342
 343     _cmsVEC3init(&In, Value -> X, Value -> Y, Value -> Z);
 344     _cmsMAT3eval(&Out, &Bradford, &In);
 345
 346     Result -> X = Out.n[0];
 347     Result -> Y = Out.n[1];
 348     Result -> Z = Out.n[2];
 349
 350     return TRUE;
 351 }
 352
 353