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