Home | History | Annotate | Download | only in src
      1 //---------------------------------------------------------------------------------
      2 //
      3 //  Little Color Management System
      4 //  Copyright (c) 1998-2014 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 
     80         y = -3.000*(x*x) + 2.870*x - 0.275;
     81 
     82         // wave factors (not used, but here for futures extensions)
     83 
     84         // M1 = (-1.3515 - 1.7703*x + 5.9114 *y)/(0.0241 + 0.2562*x - 0.7341*y);
     85         // M2 = (0.0300 - 31.4424*x + 30.0717*y)/(0.0241 + 0.2562*x - 0.7341*y);
     86 
     87         WhitePoint -> x = x;
     88         WhitePoint -> y = y;
     89         WhitePoint -> Y = 1.0;
     90 
     91         return TRUE;
     92 }
     93 
     94 
     95 
     96 typedef struct {
     97 
     98     cmsFloat64Number mirek;  // temp (in microreciprocal kelvin)
     99     cmsFloat64Number ut;     // u coord of intersection w/ blackbody locus
    100     cmsFloat64Number vt;     // v coord of intersection w/ blackbody locus
    101     cmsFloat64Number tt;     // slope of ISOTEMPERATURE. line
    102 
    103     } ISOTEMPERATURE;
    104 
    105 static ISOTEMPERATURE isotempdata[] = {
    106 //  {Mirek, Ut,       Vt,      Tt      }
    107     {0,     0.18006,  0.26352,  -0.24341},
    108     {10,    0.18066,  0.26589,  -0.25479},
    109     {20,    0.18133,  0.26846,  -0.26876},
    110     {30,    0.18208,  0.27119,  -0.28539},
    111     {40,    0.18293,  0.27407,  -0.30470},
    112     {50,    0.18388,  0.27709,  -0.32675},
    113     {60,    0.18494,  0.28021,  -0.35156},
    114     {70,    0.18611,  0.28342,  -0.37915},
    115     {80,    0.18740,  0.28668,  -0.40955},
    116     {90,    0.18880,  0.28997,  -0.44278},
    117     {100,   0.19032,  0.29326,  -0.47888},
    118     {125,   0.19462,  0.30141,  -0.58204},
    119     {150,   0.19962,  0.30921,  -0.70471},
    120     {175,   0.20525,  0.31647,  -0.84901},
    121     {200,   0.21142,  0.32312,  -1.0182 },
    122     {225,   0.21807,  0.32909,  -1.2168 },
    123     {250,   0.22511,  0.33439,  -1.4512 },
    124     {275,   0.23247,  0.33904,  -1.7298 },
    125     {300,   0.24010,  0.34308,  -2.0637 },
    126     {325,   0.24702,  0.34655,  -2.4681 },
    127     {350,   0.25591,  0.34951,  -2.9641 },
    128     {375,   0.26400,  0.35200,  -3.5814 },
    129     {400,   0.27218,  0.35407,  -4.3633 },
    130     {425,   0.28039,  0.35577,  -5.3762 },
    131     {450,   0.28863,  0.35714,  -6.7262 },
    132     {475,   0.29685,  0.35823,  -8.5955 },
    133     {500,   0.30505,  0.35907,  -11.324 },
    134     {525,   0.31320,  0.35968,  -15.628 },
    135     {550,   0.32129,  0.36011,  -23.325 },
    136     {575,   0.32931,  0.36038,  -40.770 },
    137     {600,   0.33724,  0.36051,  -116.45  }
    138 };
    139 
    140 #define NISO sizeof(isotempdata)/sizeof(ISOTEMPERATURE)
    141 
    142 
    143 // Robertson's method
    144 cmsBool  CMSEXPORT cmsTempFromWhitePoint(cmsFloat64Number* TempK, const cmsCIExyY* WhitePoint)
    145 {
    146     cmsUInt32Number j;
    147     cmsFloat64Number us,vs;
    148     cmsFloat64Number uj,vj,tj,di,dj,mi,mj;
    149     cmsFloat64Number xs, ys;
    150 
    151     _cmsAssert(WhitePoint != NULL);
    152     _cmsAssert(TempK != NULL);
    153 
    154     di = mi = 0;
    155     xs = WhitePoint -> x;
    156     ys = WhitePoint -> y;
    157 
    158     // convert (x,y) to CIE 1960 (u,WhitePoint)
    159 
    160     us = (2*xs) / (-xs + 6*ys + 1.5);
    161     vs = (3*ys) / (-xs + 6*ys + 1.5);
    162 
    163 
    164     for (j=0; j < NISO; j++) {
    165 
    166         uj = isotempdata[j].ut;
    167         vj = isotempdata[j].vt;
    168         tj = isotempdata[j].tt;
    169         mj = isotempdata[j].mirek;
    170 
    171         dj = ((vs - vj) - tj * (us - uj)) / sqrt(1.0 + tj * tj);
    172 
    173         if ((j != 0) && (di/dj < 0.0)) {
    174 
    175             // Found a match
    176             *TempK = 1000000.0 / (mi + (di / (di - dj)) * (mj - mi));
    177             return TRUE;
    178         }
    179 
    180         di = dj;
    181         mi = mj;
    182     }
    183 
    184     // Not found
    185     return FALSE;
    186 }
    187 
    188 
    189 // Compute chromatic adaptation matrix using Chad as cone matrix
    190 
    191 static
    192 cmsBool ComputeChromaticAdaptation(cmsMAT3* Conversion,
    193                                 const cmsCIEXYZ* SourceWhitePoint,
    194                                 const cmsCIEXYZ* DestWhitePoint,
    195                                 const cmsMAT3* Chad)
    196 
    197 {
    198 
    199     cmsMAT3 Chad_Inv;
    200     cmsVEC3 ConeSourceXYZ, ConeSourceRGB;
    201     cmsVEC3 ConeDestXYZ, ConeDestRGB;
    202     cmsMAT3 Cone, Tmp;
    203 
    204 
    205     Tmp = *Chad;
    206     if (!_cmsMAT3inverse(&Tmp, &Chad_Inv)) return FALSE;
    207 
    208     _cmsVEC3init(&ConeSourceXYZ, SourceWhitePoint -> X,
    209                              SourceWhitePoint -> Y,
    210                              SourceWhitePoint -> Z);
    211 
    212     _cmsVEC3init(&ConeDestXYZ,   DestWhitePoint -> X,
    213                              DestWhitePoint -> Y,
    214                              DestWhitePoint -> Z);
    215 
    216     _cmsMAT3eval(&ConeSourceRGB, Chad, &ConeSourceXYZ);
    217     _cmsMAT3eval(&ConeDestRGB,   Chad, &ConeDestXYZ);
    218 
    219     // Build matrix
    220     _cmsVEC3init(&Cone.v[0], ConeDestRGB.n[0]/ConeSourceRGB.n[0],    0.0,  0.0);
    221     _cmsVEC3init(&Cone.v[1], 0.0,   ConeDestRGB.n[1]/ConeSourceRGB.n[1],   0.0);
    222     _cmsVEC3init(&Cone.v[2], 0.0,   0.0,   ConeDestRGB.n[2]/ConeSourceRGB.n[2]);
    223 
    224 
    225     // Normalize
    226     _cmsMAT3per(&Tmp, &Cone, Chad);
    227     _cmsMAT3per(Conversion, &Chad_Inv, &Tmp);
    228 
    229     return TRUE;
    230 }
    231 
    232 // Returns the final chrmatic adaptation from illuminant FromIll to Illuminant ToIll
    233 // The cone matrix can be specified in ConeMatrix. If NULL, Bradford is assumed
    234 cmsBool  _cmsAdaptationMatrix(cmsMAT3* r, const cmsMAT3* ConeMatrix, const cmsCIEXYZ* FromIll, const cmsCIEXYZ* ToIll)
    235 {
    236     cmsMAT3 LamRigg   = {{ // Bradford matrix
    237         {{  0.8951,  0.2664, -0.1614 }},
    238         {{ -0.7502,  1.7135,  0.0367 }},
    239         {{  0.0389, -0.0685,  1.0296 }}
    240     }};
    241 
    242     if (ConeMatrix == NULL)
    243         ConeMatrix = &LamRigg;
    244 
    245     return ComputeChromaticAdaptation(r, FromIll, ToIll, ConeMatrix);
    246 }
    247 
    248 // Same as anterior, but assuming D50 destination. White point is given in xyY
    249 static
    250 cmsBool _cmsAdaptMatrixToD50(cmsMAT3* r, const cmsCIExyY* SourceWhitePt)
    251 {
    252     cmsCIEXYZ Dn;
    253     cmsMAT3 Bradford;
    254     cmsMAT3 Tmp;
    255 
    256     cmsxyY2XYZ(&Dn, SourceWhitePt);
    257 
    258     if (!_cmsAdaptationMatrix(&Bradford, NULL, &Dn, cmsD50_XYZ())) return FALSE;
    259 
    260     Tmp = *r;
    261     _cmsMAT3per(r, &Bradford, &Tmp);
    262 
    263     return TRUE;
    264 }
    265 
    266 // Build a White point, primary chromas transfer matrix from RGB to CIE XYZ
    267 // This is just an approximation, I am not handling all the non-linear
    268 // aspects of the RGB to XYZ process, and assumming that the gamma correction
    269 // has transitive property in the tranformation chain.
    270 //
    271 // the alghoritm:
    272 //
    273 //            - First I build the absolute conversion matrix using
    274 //              primaries in XYZ. This matrix is next inverted
    275 //            - Then I eval the source white point across this matrix
    276 //              obtaining the coeficients of the transformation
    277 //            - Then, I apply these coeficients to the original matrix
    278 //
    279 cmsBool _cmsBuildRGB2XYZtransferMatrix(cmsMAT3* r, const cmsCIExyY* WhitePt, const cmsCIExyYTRIPLE* Primrs)
    280 {
    281     cmsVEC3 WhitePoint, Coef;
    282     cmsMAT3 Result, Primaries;
    283     cmsFloat64Number xn, yn;
    284     cmsFloat64Number xr, yr;
    285     cmsFloat64Number xg, yg;
    286     cmsFloat64Number xb, yb;
    287 
    288     xn = WhitePt -> x;
    289     yn = WhitePt -> y;
    290     xr = Primrs -> Red.x;
    291     yr = Primrs -> Red.y;
    292     xg = Primrs -> Green.x;
    293     yg = Primrs -> Green.y;
    294     xb = Primrs -> Blue.x;
    295     yb = Primrs -> Blue.y;
    296 
    297     // Build Primaries matrix
    298     _cmsVEC3init(&Primaries.v[0], xr,        xg,         xb);
    299     _cmsVEC3init(&Primaries.v[1], yr,        yg,         yb);
    300     _cmsVEC3init(&Primaries.v[2], (1-xr-yr), (1-xg-yg),  (1-xb-yb));
    301 
    302 
    303     // Result = Primaries ^ (-1) inverse matrix
    304     if (!_cmsMAT3inverse(&Primaries, &Result))
    305         return FALSE;
    306 
    307 
    308     _cmsVEC3init(&WhitePoint, xn/yn, 1.0, (1.0-xn-yn)/yn);
    309 
    310     // Across inverse primaries ...
    311     _cmsMAT3eval(&Coef, &Result, &WhitePoint);
    312 
    313     // Give us the Coefs, then I build transformation matrix
    314     _cmsVEC3init(&r -> v[0], Coef.n[VX]*xr,          Coef.n[VY]*xg,          Coef.n[VZ]*xb);
    315     _cmsVEC3init(&r -> v[1], Coef.n[VX]*yr,          Coef.n[VY]*yg,          Coef.n[VZ]*yb);
    316     _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));
    317 
    318 
    319     return _cmsAdaptMatrixToD50(r, WhitePt);
    320 
    321 }
    322 
    323 
    324 // Adapts a color to a given illuminant. Original color is expected to have
    325 // a SourceWhitePt white point.
    326 cmsBool CMSEXPORT cmsAdaptToIlluminant(cmsCIEXYZ* Result,
    327                                        const cmsCIEXYZ* SourceWhitePt,
    328                                        const cmsCIEXYZ* Illuminant,
    329                                        const cmsCIEXYZ* Value)
    330 {
    331     cmsMAT3 Bradford;
    332     cmsVEC3 In, Out;
    333 
    334     _cmsAssert(Result != NULL);
    335     _cmsAssert(SourceWhitePt != NULL);
    336     _cmsAssert(Illuminant != NULL);
    337     _cmsAssert(Value != NULL);
    338 
    339     if (!_cmsAdaptationMatrix(&Bradford, NULL, SourceWhitePt, Illuminant)) return FALSE;
    340 
    341     _cmsVEC3init(&In, Value -> X, Value -> Y, Value -> Z);
    342     _cmsMAT3eval(&Out, &Bradford, &In);
    343 
    344     Result -> X = Out.n[0];
    345     Result -> Y = Out.n[1];
    346     Result -> Z = Out.n[2];
    347 
    348     return TRUE;
    349 }
    350