jdk/src/share/native/sun/java2d/cmm/lcms/cmscam97.c
changeset 6483 8f9f87a564f6
parent 6390 a6442d6bc38a
parent 6482 0f6a4442b29e
child 6488 404882083757
equal deleted inserted replaced
6390:a6442d6bc38a 6483:8f9f87a564f6
     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 //  Little cms
       
    32 //  Copyright (C) 1998-2007 Marti Maria
       
    33 //
       
    34 // Permission is hereby granted, free of charge, to any person obtaining
       
    35 // a copy of this software and associated documentation files (the "Software"),
       
    36 // to deal in the Software without restriction, including without limitation
       
    37 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
       
    38 // and/or sell copies of the Software, and to permit persons to whom the Software
       
    39 // is furnished to do so, subject to the following conditions:
       
    40 //
       
    41 // The above copyright notice and this permission notice shall be included in
       
    42 // all copies or substantial portions of the Software.
       
    43 //
       
    44 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
       
    45 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
       
    46 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
       
    47 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
       
    48 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
       
    49 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
       
    50 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
       
    51 
       
    52 
       
    53 #include "lcms.h"
       
    54 
       
    55 
       
    56 /*
       
    57 typedef struct {
       
    58                double J;
       
    59                double C;
       
    60                double h;
       
    61 
       
    62                } cmsJCh, FAR* LPcmsJCh;
       
    63 
       
    64 
       
    65 #define AVG_SURROUND_4     0
       
    66 #define AVG_SURROUND       1
       
    67 #define DIM_SURROUND       2
       
    68 #define DARK_SURROUND      3
       
    69 #define CUTSHEET_SURROUND  4
       
    70 
       
    71 
       
    72 typedef struct {
       
    73 
       
    74               cmsCIEXYZ whitePoint;
       
    75               double    Yb;
       
    76               double    La;
       
    77               int       surround;
       
    78               double    D_value;
       
    79 
       
    80               } cmsViewingConditions, FAR* LPcmsViewingConditions;
       
    81 
       
    82 
       
    83 
       
    84 LCMSAPI LCMSHANDLE LCMSEXPORT cmsCIECAM97sInit(LPcmsViewingConditions pVC);
       
    85 LCMSAPI void   LCMSEXPORT cmsCIECAM97sDone(LCMSHANDLE hModel);
       
    86 LCMSAPI void   LCMSEXPORT cmsCIECAM97sForward(LCMSHANDLE hModel, LPcmsCIEXYZ pIn, LPcmsJCh pOut);
       
    87 LCMSAPI void   LCMSEXPORT cmsCIECAM97sReverse(LCMSHANDLE hModel, LPcmsJCh pIn,    LPcmsCIEXYZ pOut);
       
    88 
       
    89 */
       
    90 
       
    91 // ---------- Implementation --------------------------------------------
       
    92 
       
    93 // #define USE_CIECAM97s2  1
       
    94 
       
    95 #ifdef USE_CIECAM97s2
       
    96 
       
    97 #       define NOISE_CONSTANT   3.05
       
    98 #else
       
    99 #       define NOISE_CONSTANT   2.05
       
   100 #endif
       
   101 
       
   102 
       
   103 /*
       
   104   The model input data are the adapting field luminance in cd/m2
       
   105   (normally taken to be 20% of the luminance of white in the adapting field),
       
   106   LA , the relative tristimulus values of the stimulus, XYZ, the relative
       
   107   tristimulus values of white in the same viewing conditions, Xw Yw Zw ,
       
   108   and the relative luminance of the background, Yb . Relative tristimulus
       
   109   values should be expressed on a scale from Y = 0 for a perfect black
       
   110   to Y = 100 for a perfect reflecting diffuser. Additionally, the
       
   111   parameters c, for the impact of surround, Nc , a chromatic induction factor,
       
   112   and F, a factor for degree of adaptation, must be selected according to the
       
   113   guidelines in table
       
   114 
       
   115   All CIE tristimulus values are obtained using the CIE 1931
       
   116   Standard Colorimetric Observer (2°).
       
   117 
       
   118 */
       
   119 
       
   120 typedef struct {
       
   121 
       
   122     cmsCIEXYZ WP;
       
   123     int surround;
       
   124     int calculate_D;
       
   125 
       
   126     double  Yb;         // rel. luminance of background
       
   127 
       
   128     cmsCIEXYZ RefWhite;
       
   129 
       
   130     double La;    // The adapting field luminance in cd/m2
       
   131 
       
   132     double c;     // Impact of surround
       
   133     double Nc;    // Chromatic induction factor
       
   134     double Fll;   // Lightness contrast factor (Removed on rev 2)
       
   135     double F;     // Degree of adaptation
       
   136 
       
   137 
       
   138     double k;
       
   139     double Fl;
       
   140 
       
   141     double Nbb;  // The background and chromatic brightness induction factors.
       
   142     double Ncb;
       
   143     double z;    // base exponential nonlinearity
       
   144     double n;    // background induction factor
       
   145     double D;
       
   146 
       
   147     MAT3 MlamRigg;
       
   148     MAT3 MlamRigg_1;
       
   149 
       
   150     MAT3 Mhunt;
       
   151     MAT3 Mhunt_1;
       
   152 
       
   153     MAT3 Mhunt_x_MlamRigg_1;
       
   154     MAT3 MlamRigg_x_Mhunt_1;
       
   155 
       
   156 
       
   157     VEC3 RGB_subw;
       
   158     VEC3 RGB_subw_prime;
       
   159 
       
   160     double p;
       
   161 
       
   162     VEC3 RGB_subwc;
       
   163 
       
   164     VEC3 RGB_subaw_prime;
       
   165     double A_subw;
       
   166     double Q_subw;
       
   167 
       
   168     } cmsCIECAM97s,FAR *LPcmsCIECAM97s;
       
   169 
       
   170 
       
   171 
       
   172 // Free model structure
       
   173 
       
   174 LCMSAPI void LCMSEXPORT cmsCIECAM97sDone(LCMSHANDLE hModel)
       
   175 {
       
   176     LPcmsCIECAM97s lpMod = (LPcmsCIECAM97s) (LPSTR) hModel;
       
   177     if (lpMod) _cmsFree(lpMod);
       
   178 }
       
   179 
       
   180 // Partial discounting for adaptation degree computation
       
   181 
       
   182 static
       
   183 double discount(double d, double chan)
       
   184 {
       
   185     return (d * chan + 1 - d);
       
   186 }
       
   187 
       
   188 
       
   189 // This routine does model exponential nonlinearity on the short wavelenght
       
   190 // sensitive channel. On CIECAM97s rev 2 this has been reverted to linear.
       
   191 
       
   192 static
       
   193 void FwAdaptationDegree(LPcmsCIECAM97s lpMod, LPVEC3 RGBc, LPVEC3 RGB)
       
   194 {
       
   195 
       
   196 
       
   197 #ifdef USE_CIECAM97s2
       
   198     RGBc->n[0] = RGB->n[0]* discount(lpMod->D, 100.0/lpMod->RGB_subw.n[0]);
       
   199     RGBc->n[1] = RGB->n[1]* discount(lpMod->D, 100.0/lpMod->RGB_subw.n[1]);
       
   200     RGBc->n[2] = RGB->n[2]* discount(lpMod->D, 100.0/lpMod->RGB_subw.n[2]);
       
   201 #else
       
   202 
       
   203     RGBc->n[0] = RGB->n[0]* discount(lpMod->D, 1.0/lpMod->RGB_subw.n[0]);
       
   204     RGBc->n[1] = RGB->n[1]* discount(lpMod->D, 1.0/lpMod->RGB_subw.n[1]);
       
   205 
       
   206     RGBc->n[2] = pow(fabs(RGB->n[2]), lpMod ->p) * discount(lpMod->D, (1.0/pow(lpMod->RGB_subw.n[2], lpMod->p)));
       
   207 
       
   208     // If B happens to be negative, Then Bc is also set to be negative
       
   209 
       
   210     if (RGB->n[2] < 0)
       
   211            RGBc->n[2] = -RGBc->n[2];
       
   212 #endif
       
   213 }
       
   214 
       
   215 
       
   216 static
       
   217 void RvAdaptationDegree(LPcmsCIECAM97s lpMod, LPVEC3 RGBc, LPVEC3 RGB)
       
   218 {
       
   219 
       
   220 
       
   221 #ifdef USE_CIECAM97s2
       
   222     RGBc->n[0] = RGB->n[0]/discount(lpMod->D, 100.0/lpMod->RGB_subw.n[0]);
       
   223     RGBc->n[1] = RGB->n[1]/discount(lpMod->D, 100.0/lpMod->RGB_subw.n[1]);
       
   224     RGBc->n[2] = RGB->n[2]/discount(lpMod->D, 100.0/lpMod->RGB_subw.n[2]);
       
   225 #else
       
   226 
       
   227     RGBc->n[0] = RGB->n[0]/discount(lpMod->D, 1.0/lpMod->RGB_subw.n[0]);
       
   228     RGBc->n[1] = RGB->n[1]/discount(lpMod->D, 1.0/lpMod->RGB_subw.n[1]);
       
   229     RGBc->n[2] = pow(fabs(RGB->n[2]), 1.0/lpMod->p)/pow(discount(lpMod->D, 1.0/pow(lpMod->RGB_subw.n[2], lpMod->p)), 1.0/lpMod->p);
       
   230     if (RGB->n[2] < 0)
       
   231            RGBc->n[2] = -RGBc->n[2];
       
   232 #endif
       
   233 }
       
   234 
       
   235 
       
   236 
       
   237 static
       
   238 void PostAdaptationConeResponses(LPcmsCIECAM97s lpMod, LPVEC3 RGBa_prime, LPVEC3 RGBprime)
       
   239 {
       
   240      if (RGBprime->n[0]>=0.0) {
       
   241 
       
   242             RGBa_prime->n[0]=((40.0*pow(lpMod -> Fl * RGBprime->n[0]/100.0, 0.73))/(pow(lpMod -> Fl * RGBprime->n[0]/100.0, 0.73)+2))+1;
       
   243      }
       
   244      else
       
   245      {
       
   246             RGBa_prime->n[0]=((-40.0*pow((-lpMod -> Fl * RGBprime->n[0])/100.0, 0.73))/(pow((-lpMod -> Fl * RGBprime->n[0])/100.0, 0.73)+2))+1;
       
   247      }
       
   248 
       
   249      if (RGBprime->n[1]>=0.0)
       
   250      {
       
   251             RGBa_prime->n[1]=((40.0*pow(lpMod -> Fl * RGBprime->n[1]/100.0, 0.73))/(pow(lpMod -> Fl * RGBprime->n[1]/100.0, 0.73)+2))+1;
       
   252      }
       
   253      else
       
   254      {
       
   255             RGBa_prime->n[1]=((-40.0*pow((-lpMod -> Fl * RGBprime->n[1])/100.0, 0.73))/(pow((-lpMod -> Fl * RGBprime->n[1])/100.0, 0.73)+2))+1;
       
   256      }
       
   257 
       
   258      if (RGBprime->n[2]>=0.0)
       
   259      {
       
   260             RGBa_prime->n[2]=((40.0*pow(lpMod -> Fl * RGBprime->n[2]/100.0, 0.73))/(pow(lpMod -> Fl * RGBprime->n[2]/100.0, 0.73)+2))+1;
       
   261      }
       
   262      else
       
   263      {
       
   264             RGBa_prime->n[2]=((-40.0*pow((-lpMod -> Fl * RGBprime->n[2])/100.0, 0.73))/(pow((-lpMod -> Fl * RGBprime->n[2])/100.0, 0.73)+2))+1;
       
   265      }
       
   266 }
       
   267 
       
   268 
       
   269 // Compute hue quadrature, eccentricity factor, e
       
   270 
       
   271 static
       
   272 void ComputeHueQuadrature(double h, double* H, double* e)
       
   273 {
       
   274 
       
   275 
       
   276 #define IRED    0
       
   277 #define IYELLOW 1
       
   278 #define IGREEN  2
       
   279 #define IBLUE   3
       
   280 
       
   281       double e_tab[] = {0.8, 0.7, 1.0, 1.2};
       
   282       double H_tab[] = {  0, 100, 200, 300};
       
   283       int p1, p2;
       
   284       double e1, e2, h1, h2;
       
   285 
       
   286 
       
   287        if (h >= 20.14 && h < 90.0) { // Red
       
   288 
       
   289                         p1 = IRED;
       
   290                         p2 = IYELLOW;
       
   291        }
       
   292        else
       
   293        if (h >= 90.0 && h < 164.25) { // Yellow
       
   294 
       
   295                         p1 = IYELLOW;
       
   296                         p2 = IGREEN;
       
   297        }
       
   298        else
       
   299        if (h >= 164.25 && h < 237.53) { // Green
       
   300 
       
   301                         p1 = IGREEN;
       
   302                         p2 = IBLUE;       }
       
   303        else {                         // Blue
       
   304 
       
   305                         p1 = IBLUE;
       
   306                         p2 = IRED;
       
   307        }
       
   308 
       
   309        e1 = e_tab[p1]; e2 = e_tab[p2];
       
   310        h1 = H_tab[p1]; h2 = H_tab[p2];
       
   311 
       
   312 
       
   313 
       
   314        *e = e1 + ((e2-e1)*(h-h1)/(h2 - h1));
       
   315        *H = h1 + (100. * (h - h1) / e1) / ((h - h1)/e1 + (h2 - h) / e2);
       
   316 
       
   317 #undef IRED
       
   318 #undef IYELLOW
       
   319 #undef IGREEN
       
   320 #undef IBLUE
       
   321 
       
   322 }
       
   323 
       
   324 
       
   325 
       
   326 
       
   327 
       
   328 
       
   329 LCMSAPI LCMSHANDLE LCMSEXPORT cmsCIECAM97sInit(LPcmsViewingConditions pVC)
       
   330 {
       
   331     LPcmsCIECAM97s lpMod;
       
   332     VEC3 tmp;
       
   333 
       
   334     if((lpMod = (LPcmsCIECAM97s) _cmsMalloc(sizeof(cmsCIECAM97s))) == NULL) {
       
   335         return (LCMSHANDLE) NULL;
       
   336     }
       
   337 
       
   338 
       
   339     lpMod->WP.X = pVC->whitePoint.X;
       
   340     lpMod->WP.Y = pVC->whitePoint.Y;
       
   341     lpMod->WP.Z = pVC->whitePoint.Z;
       
   342 
       
   343     lpMod->Yb   = pVC->Yb;
       
   344     lpMod->La   = pVC->La;
       
   345 
       
   346     lpMod->surround = pVC->surround;
       
   347 
       
   348     lpMod->RefWhite.X = 100.0;
       
   349     lpMod->RefWhite.Y = 100.0;
       
   350     lpMod->RefWhite.Z = 100.0;
       
   351 
       
   352 #ifdef USE_CIECAM97s2
       
   353 
       
   354     VEC3init(&lpMod->MlamRigg.v[0],  0.8562, 0.3372, -0.1934);
       
   355     VEC3init(&lpMod->MlamRigg.v[1], -0.8360, 1.8327,  0.0033);
       
   356     VEC3init(&lpMod->MlamRigg.v[2],  0.0357,-0.0469,  1.0112);
       
   357 
       
   358     VEC3init(&lpMod->MlamRigg_1.v[0], 0.9874, -0.1768, 0.1894);
       
   359     VEC3init(&lpMod->MlamRigg_1.v[1], 0.4504,  0.4649, 0.0846);
       
   360     VEC3init(&lpMod->MlamRigg_1.v[2],-0.0139,  0.0278, 0.9861);
       
   361 
       
   362 #else
       
   363     // Bradford transform: Lam-Rigg cone responses
       
   364     VEC3init(&lpMod->MlamRigg.v[0],  0.8951,  0.2664, -0.1614);
       
   365     VEC3init(&lpMod->MlamRigg.v[1], -0.7502,  1.7135,  0.0367);
       
   366     VEC3init(&lpMod->MlamRigg.v[2],  0.0389, -0.0685,  1.0296);
       
   367 
       
   368 
       
   369     // Inverse of Lam-Rigg
       
   370     VEC3init(&lpMod->MlamRigg_1.v[0],  0.98699, -0.14705,  0.15996);
       
   371     VEC3init(&lpMod->MlamRigg_1.v[1],  0.43231,  0.51836,  0.04929);
       
   372     VEC3init(&lpMod->MlamRigg_1.v[2], -0.00853,  0.04004,  0.96849);
       
   373 
       
   374 #endif
       
   375 
       
   376     // Hunt-Pointer-Estevez cone responses
       
   377     VEC3init(&lpMod->Mhunt.v[0],   0.38971,  0.68898, -0.07868);
       
   378     VEC3init(&lpMod->Mhunt.v[1],  -0.22981,  1.18340,  0.04641);
       
   379     VEC3init(&lpMod->Mhunt.v[2],   0.0,      0.0,      1.0);
       
   380 
       
   381     // Inverse of Hunt-Pointer-Estevez
       
   382     VEC3init(&lpMod->Mhunt_1.v[0],     1.91019, -1.11214, 0.20195);
       
   383     VEC3init(&lpMod->Mhunt_1.v[1],     0.37095,  0.62905, 0.0);
       
   384     VEC3init(&lpMod->Mhunt_1.v[2],     0.0,      0.0,     1.0);
       
   385 
       
   386 
       
   387     if (pVC->D_value == -1.0)
       
   388           lpMod->calculate_D = 1;
       
   389     else
       
   390     if (pVC->D_value == -2.0)
       
   391            lpMod->calculate_D = 2;
       
   392     else {
       
   393         lpMod->calculate_D = 0;
       
   394         lpMod->D = pVC->D_value;
       
   395     }
       
   396 
       
   397    // Table I (revised)
       
   398 
       
   399    switch (lpMod->surround) {
       
   400 
       
   401     case AVG_SURROUND_4:
       
   402        lpMod->F = 1.0;
       
   403        lpMod->c = 0.69;
       
   404        lpMod->Fll = 0.0;    // Not included on Rev 2
       
   405        lpMod->Nc = 1.0;
       
   406        break;
       
   407     case AVG_SURROUND:
       
   408        lpMod->F = 1.0;
       
   409        lpMod->c = 0.69;
       
   410        lpMod->Fll = 1.0;
       
   411        lpMod->Nc = 1.0;
       
   412        break;
       
   413     case DIM_SURROUND:
       
   414        lpMod->F = 0.99;
       
   415        lpMod->c = 0.59;
       
   416        lpMod->Fll = 1.0;
       
   417        lpMod->Nc = 0.95;
       
   418        break;
       
   419     case DARK_SURROUND:
       
   420        lpMod->F = 0.9;
       
   421        lpMod->c = 0.525;
       
   422        lpMod->Fll = 1.0;
       
   423        lpMod->Nc = 0.8;
       
   424        break;
       
   425     case CUTSHEET_SURROUND:
       
   426        lpMod->F = 0.9;
       
   427        lpMod->c = 0.41;
       
   428        lpMod->Fll = 1.0;
       
   429        lpMod->Nc = 0.8;
       
   430        break;
       
   431     default:
       
   432        lpMod->F = 1.0;
       
   433        lpMod->c = 0.69;
       
   434        lpMod->Fll = 1.0;
       
   435        lpMod->Nc = 1.0;
       
   436        break;
       
   437     }
       
   438 
       
   439     lpMod->k = 1 / (5 * lpMod->La  + 1);
       
   440     lpMod->Fl = lpMod->La * pow(lpMod->k, 4) + 0.1*pow(1 - pow(lpMod->k, 4), 2.0) * pow(5*lpMod->La, 1.0/3.0);
       
   441 
       
   442     if (lpMod->calculate_D > 0) {
       
   443 
       
   444        lpMod->D = lpMod->F * (1 - 1 / (1 + 2*pow(lpMod->La, 0.25) + pow(lpMod->La, 2)/300.0));
       
   445        if (lpMod->calculate_D > 1)
       
   446            lpMod->D = (lpMod->D + 1.0) / 2;
       
   447     }
       
   448 
       
   449 
       
   450     // RGB_subw = [MlamRigg][WP/YWp]
       
   451 #ifdef USE_CIECAM97s2
       
   452     MAT3eval(&lpMod -> RGB_subw, &lpMod -> MlamRigg, &lpMod -> WP);
       
   453 #else
       
   454     VEC3divK(&tmp, (LPVEC3) &lpMod -> WP, lpMod->WP.Y);
       
   455     MAT3eval(&lpMod -> RGB_subw, &lpMod -> MlamRigg, &tmp);
       
   456 #endif
       
   457 
       
   458 
       
   459 
       
   460     MAT3per(&lpMod -> Mhunt_x_MlamRigg_1,   &lpMod -> Mhunt,   &lpMod->MlamRigg_1  );
       
   461     MAT3per(&lpMod -> MlamRigg_x_Mhunt_1,   &lpMod -> MlamRigg, &lpMod -> Mhunt_1  );
       
   462 
       
   463     // p is used on forward model
       
   464     lpMod->p = pow(lpMod->RGB_subw.n[2], 0.0834);
       
   465 
       
   466     FwAdaptationDegree(lpMod, &lpMod->RGB_subwc, &lpMod->RGB_subw);
       
   467 
       
   468 #if USE_CIECAM97s2
       
   469     MAT3eval(&lpMod->RGB_subw_prime, &lpMod->Mhunt_x_MlamRigg_1, &lpMod -> RGB_subwc);
       
   470 #else
       
   471     VEC3perK(&tmp, &lpMod -> RGB_subwc, lpMod->WP.Y);
       
   472     MAT3eval(&lpMod->RGB_subw_prime, &lpMod->Mhunt_x_MlamRigg_1, &tmp);
       
   473 #endif
       
   474 
       
   475     lpMod->n = lpMod-> Yb / lpMod-> WP.Y;
       
   476 
       
   477     lpMod->z = 1 + lpMod->Fll * sqrt(lpMod->n);
       
   478     lpMod->Nbb = lpMod->Ncb = 0.725 / pow(lpMod->n, 0.2);
       
   479 
       
   480     PostAdaptationConeResponses(lpMod, &lpMod->RGB_subaw_prime, &lpMod->RGB_subw_prime);
       
   481 
       
   482     lpMod->A_subw=lpMod->Nbb*(2.0*lpMod->RGB_subaw_prime.n[0]+lpMod->RGB_subaw_prime.n[1]+lpMod->RGB_subaw_prime.n[2]/20.0-NOISE_CONSTANT);
       
   483 
       
   484     return (LCMSHANDLE) lpMod;
       
   485 }
       
   486 
       
   487 
       
   488 
       
   489 
       
   490 //
       
   491 // The forward model: XYZ -> JCh
       
   492 //
       
   493 
       
   494 LCMSAPI void LCMSEXPORT cmsCIECAM97sForward(LCMSHANDLE hModel, LPcmsCIEXYZ inPtr, LPcmsJCh outPtr)
       
   495 {
       
   496 
       
   497         LPcmsCIECAM97s lpMod = (LPcmsCIECAM97s) (LPSTR) hModel;
       
   498         double a, b, h, s, H1val, es, A;
       
   499         VEC3 In, RGB, RGBc, RGBprime, RGBa_prime;
       
   500 
       
   501         if (inPtr -> Y <= 0.0) {
       
   502 
       
   503       outPtr -> J = outPtr -> C = outPtr -> h = 0.0;
       
   504           return;
       
   505         }
       
   506 
       
   507        // An initial chromatic adaptation transform is used to go from the source
       
   508        // viewing conditions to corresponding colours under the equal-energy-illuminant
       
   509        // reference viewing conditions. This is handled differently on rev 2
       
   510 
       
   511        VEC3init(&In, inPtr -> X, inPtr -> Y, inPtr -> Z);    // 2.1
       
   512 
       
   513 #ifdef USE_CIECAM97s2
       
   514        // Since the chromatic adaptation transform has been linearized, it
       
   515        // is no longer required to divide the stimulus tristimulus values
       
   516        // by their own Y tristimulus value prior to the chromatic adaptation.
       
   517 #else
       
   518        VEC3divK(&In, &In, inPtr -> Y);
       
   519 #endif
       
   520 
       
   521        MAT3eval(&RGB, &lpMod -> MlamRigg, &In);              // 2.2
       
   522 
       
   523        FwAdaptationDegree(lpMod, &RGBc, &RGB);
       
   524 
       
   525        // The post-adaptation signals for both the sample and the white are then
       
   526        // transformed from the sharpened cone responses to the Hunt-Pointer-Estevez
       
   527        // cone responses.
       
   528 #ifdef USE_CIECAM97s2
       
   529 #else
       
   530        VEC3perK(&RGBc, &RGBc, inPtr->Y);
       
   531 #endif
       
   532 
       
   533        MAT3eval(&RGBprime, &lpMod->Mhunt_x_MlamRigg_1, &RGBc);
       
   534 
       
   535        // The post-adaptation cone responses (for both the stimulus and the white)
       
   536        // are then calculated.
       
   537 
       
   538        PostAdaptationConeResponses(lpMod, &RGBa_prime, &RGBprime);
       
   539 
       
   540        // Preliminary red-green and yellow-blue opponent dimensions are calculated
       
   541 
       
   542        a = RGBa_prime.n[0] - (12.0 * RGBa_prime.n[1] / 11.0) + RGBa_prime.n[2]/11.0;
       
   543        b = (RGBa_prime.n[0] + RGBa_prime.n[1] - 2.0 * RGBa_prime.n[2]) / 9.0;
       
   544 
       
   545 
       
   546        // The CIECAM97s hue angle, h, is then calculated
       
   547        h = (180.0/M_PI)*(atan2(b, a));
       
   548 
       
   549 
       
   550        while (h < 0)
       
   551               h += 360.0;
       
   552 
       
   553        outPtr->h = h;
       
   554 
       
   555        // hue quadrature and eccentricity factors, e, are calculated
       
   556 
       
   557        ComputeHueQuadrature(h, &H1val, &es);
       
   558 
       
   559        // ComputeHueQuadrature(h, &H1val, &h1, &e1, &h2, &e2, &es);
       
   560 
       
   561 
       
   562       // The achromatic response A
       
   563       A = lpMod->Nbb * (2.0 * RGBa_prime.n[0] + RGBa_prime.n[1] + RGBa_prime.n[2]/20.0 - NOISE_CONSTANT);
       
   564 
       
   565       // CIECAM97s Lightness J
       
   566       outPtr -> J = 100.0 * pow(A / lpMod->A_subw, lpMod->c * lpMod->z);
       
   567 
       
   568       // CIECAM97s saturation s
       
   569       s =  (50 * hypot (a, b) * 100 * es * (10.0/13.0) * lpMod-> Nc * lpMod->Ncb) / (RGBa_prime.n[0] + RGBa_prime.n[1] + 1.05 * RGBa_prime.n[2]);
       
   570 
       
   571       // CIECAM97s Chroma C
       
   572 
       
   573 #ifdef USE_CIECAM97s2
       
   574       // Eq. 26 has been modified to allow accurate prediction of the Munsell chroma scales.
       
   575       outPtr->C = 0.7487 * pow(s, 0.973) * pow(outPtr->J/100.0, 0.945 * lpMod->n) * (1.64 - pow(0.29, lpMod->n));
       
   576 
       
   577 #else
       
   578       outPtr->C = 2.44 * pow(s, 0.69) * pow(outPtr->J/100.0, 0.67 * lpMod->n) * (1.64 - pow(0.29, lpMod->n));
       
   579 #endif
       
   580 }
       
   581 
       
   582 
       
   583 //
       
   584 // The reverse model JCh -> XYZ
       
   585 //
       
   586 
       
   587 
       
   588 LCMSAPI void LCMSEXPORT cmsCIECAM97sReverse(LCMSHANDLE hModel, LPcmsJCh inPtr, LPcmsCIEXYZ outPtr)
       
   589 {
       
   590     LPcmsCIECAM97s lpMod = (LPcmsCIECAM97s) (LPSTR) hModel;
       
   591     double J, C, h, A, H1val, es, s, a, b;
       
   592     double tan_h, sec_h;
       
   593     double R_suba_prime, G_suba_prime, B_suba_prime;
       
   594     double R_prime, G_prime, B_prime;
       
   595     double Y_subc, Y_prime, B_term;
       
   596     VEC3 tmp;
       
   597     VEC3 RGB_prime, RGB_subc_Y;
       
   598     VEC3 Y_over_Y_subc_RGB;
       
   599     VEC3 XYZ_primeprime_over_Y_subc;
       
   600 #ifdef USE_CIECAM92s2
       
   601     VEC3 RGBY;
       
   602     VEC3 Out;
       
   603 #endif
       
   604 
       
   605     J = inPtr->J;
       
   606     h = inPtr->h;
       
   607     C = inPtr->C;
       
   608 
       
   609     if (J <= 0) {
       
   610 
       
   611         outPtr->X =  0.0;
       
   612         outPtr->Y =  0.0;
       
   613         outPtr->Z =  0.0;
       
   614         return;
       
   615     }
       
   616 
       
   617 
       
   618 
       
   619     // (2) From J Obtain A
       
   620 
       
   621     A =  pow(J/100.0, 1/(lpMod->c * lpMod->z)) * lpMod->A_subw;
       
   622 
       
   623 
       
   624     // (3), (4), (5) Using H Determine h1, h2, e1, e2
       
   625     // e1 and h1 are the values  of e and h for the unique hue having the
       
   626     // nearest lower valur of h and e2 and h2 are the values of e and h for
       
   627     // the unique hue having the nearest higher value of h.
       
   628 
       
   629 
       
   630     ComputeHueQuadrature(h, &H1val, &es);
       
   631 
       
   632     // (7) Calculate s
       
   633 
       
   634     s = pow(C / (2.44 * pow(J/100.0, 0.67*lpMod->n) * (1.64 - pow(0.29, lpMod->n))) , (1./0.69));
       
   635 
       
   636 
       
   637     // (8) Calculate a and b.
       
   638     // NOTE: sqrt(1 + tan^2) == sec(h)
       
   639 
       
   640     tan_h = tan ((M_PI/180.)*(h));
       
   641     sec_h = sqrt(1 + tan_h * tan_h);
       
   642 
       
   643     if ((h > 90) && (h < 270))
       
   644             sec_h = -sec_h;
       
   645 
       
   646     a = s * ( A/lpMod->Nbb + NOISE_CONSTANT) / ( sec_h * 50000.0 * es * lpMod->Nc * lpMod->Ncb/ 13.0 +
       
   647            s * (11.0 / 23.0 + (108.0/23.0) * tan_h));
       
   648 
       
   649     b = a * tan_h;
       
   650 
       
   651     //(9) Calculate R'a G'a and B'a
       
   652 
       
   653     R_suba_prime = (20.0/61.0) * (A/lpMod->Nbb + NOISE_CONSTANT) + (41.0/61.0) * (11.0/23.0) * a + (288.0/61.0) / 23.0 * b;
       
   654     G_suba_prime = (20.0/61.0) * (A/lpMod->Nbb + NOISE_CONSTANT) - (81.0/61.0) * (11.0/23.0) * a - (261.0/61.0) / 23.0 * b;
       
   655     B_suba_prime = (20.0/61.0) * (A/lpMod->Nbb + NOISE_CONSTANT) - (20.0/61.0) * (11.0/23.0) * a - (20.0/61.0) * (315.0/23.0) * b;
       
   656 
       
   657     // (10) Calculate R', G' and B'
       
   658 
       
   659     if ((R_suba_prime - 1) < 0) {
       
   660 
       
   661          R_prime = -100.0 * pow((2.0 - 2.0 * R_suba_prime) /
       
   662                             (39.0 + R_suba_prime), 1.0/0.73);
       
   663     }
       
   664     else
       
   665     {
       
   666          R_prime = 100.0 * pow((2.0 * R_suba_prime - 2.0) /
       
   667                             (41.0 - R_suba_prime), 1.0/0.73);
       
   668     }
       
   669 
       
   670     if ((G_suba_prime - 1) < 0)
       
   671     {
       
   672          G_prime = -100.0 * pow((2.0 - 2.0 * G_suba_prime) /
       
   673                             (39.0 + G_suba_prime), 1.0/0.73);
       
   674     }
       
   675     else
       
   676     {
       
   677          G_prime = 100.0 * pow((2.0 * G_suba_prime - 2.0) /
       
   678                             (41.0 - G_suba_prime), 1.0/0.73);
       
   679     }
       
   680 
       
   681     if ((B_suba_prime - 1) < 0)
       
   682     {
       
   683          B_prime = -100.0 * pow((2.0 - 2.0 * B_suba_prime) /
       
   684                             (39.0 + B_suba_prime), 1.0/0.73);
       
   685     }
       
   686     else
       
   687     {
       
   688          B_prime = 100.0 * pow((2.0 * B_suba_prime - 2.0) /
       
   689                             (41.0 - B_suba_prime), 1.0/0.73);
       
   690     }
       
   691 
       
   692 
       
   693     // (11) Calculate RcY, GcY and BcY
       
   694 
       
   695     VEC3init(&RGB_prime, R_prime, G_prime, B_prime);
       
   696     VEC3divK(&tmp, &RGB_prime, lpMod -> Fl);
       
   697 
       
   698     MAT3eval(&RGB_subc_Y, &lpMod->MlamRigg_x_Mhunt_1, &tmp);
       
   699 
       
   700 
       
   701 
       
   702 
       
   703 #ifdef USE_CIECAM97s2
       
   704 
       
   705        // (12)
       
   706 
       
   707 
       
   708            RvAdaptationDegree(lpMod, &RGBY, &RGB_subc_Y);
       
   709            MAT3eval(&Out, &lpMod->MlamRigg_1, &RGBY);
       
   710 
       
   711            outPtr -> X = Out.n[0];
       
   712            outPtr -> Y = Out.n[1];
       
   713            outPtr -> Z = Out.n[2];
       
   714 
       
   715 #else
       
   716 
       
   717            // (12) Calculate Yc
       
   718 
       
   719        Y_subc = 0.43231*RGB_subc_Y.n[0]+0.51836*RGB_subc_Y.n[1]+0.04929*RGB_subc_Y.n[2];
       
   720 
       
   721            // (13) Calculate (Y/Yc)R, (Y/Yc)G and (Y/Yc)B
       
   722 
       
   723            VEC3divK(&RGB_subc_Y, &RGB_subc_Y, Y_subc);
       
   724            RvAdaptationDegree(lpMod, &Y_over_Y_subc_RGB, &RGB_subc_Y);
       
   725 
       
   726            // (14) Calculate Y'
       
   727        Y_prime = 0.43231*(Y_over_Y_subc_RGB.n[0]*Y_subc) + 0.51836*(Y_over_Y_subc_RGB.n[1]*Y_subc) + 0.04929 * (Y_over_Y_subc_RGB.n[2]*Y_subc);
       
   728 
       
   729            if (Y_prime < 0 || Y_subc < 0)
       
   730            {
       
   731                 // Discard to near black point
       
   732 
       
   733                 outPtr -> X = 0;
       
   734                 outPtr -> Y = 0;
       
   735                 outPtr -> Z = 0;
       
   736                 return;
       
   737            }
       
   738 
       
   739        B_term = pow(Y_prime / Y_subc, (1.0 / lpMod->p) - 1);
       
   740 
       
   741           // (15) Calculate X'', Y'' and Z''
       
   742            Y_over_Y_subc_RGB.n[2] /= B_term;
       
   743            MAT3eval(&XYZ_primeprime_over_Y_subc, &lpMod->MlamRigg_1, &Y_over_Y_subc_RGB);
       
   744 
       
   745            outPtr->X =  XYZ_primeprime_over_Y_subc.n[0] * Y_subc;
       
   746            outPtr->Y =  XYZ_primeprime_over_Y_subc.n[1] * Y_subc;
       
   747            outPtr->Z =  XYZ_primeprime_over_Y_subc.n[2] * Y_subc;
       
   748 #endif
       
   749 
       
   750 }