294 

   295 

   296 // ---------------------------------------------------------------------------------------------------------

   297 

   298 // Least Squares Fit of a Quadratic Curve to Data

   299 // http://www.personal.psu.edu/jhm/f90/lectures/lsq2.html

   300 

   301 static

   302 cmsFloat64Number RootOfLeastSquaresFitQuadraticCurve(int n, cmsFloat64Number x[], cmsFloat64Number y[])

   303 {

   304     double sum_x = 0, sum_x2 = 0, sum_x3 = 0, sum_x4 = 0;

   305     double sum_y = 0, sum_yx = 0, sum_yx2 = 0;

   306     double disc;

   307     int i;

   308     cmsMAT3 m;

   309     cmsVEC3 v, res;

   310 

   311     if (n < 4) return 0;

   312 

   313     for (i=0; i < n; i++) {

   314 

   315         double xn = x[i];

   316         double yn = y[i];

   317 

   318         sum_x  += xn;

   319         sum_x2 += xn*xn;

   320         sum_x3 += xn*xn*xn;

   321         sum_x4 += xn*xn*xn*xn;

   322 

   323         sum_y += yn;

   324         sum_yx += yn*xn;

   325         sum_yx2 += yn*xn*xn;

   326     }

   327 

   328     _cmsVEC3init(&m.v[0], n,      sum_x,  sum_x2);

   329     _cmsVEC3init(&m.v[1], sum_x,  sum_x2, sum_x3);

   330     _cmsVEC3init(&m.v[2], sum_x2, sum_x3, sum_x4);

   331 

   332     _cmsVEC3init(&v, sum_y, sum_yx, sum_yx2);

   333 

   334     if (!_cmsMAT3solve(&res, &m, &v)) return 0;

   335 

   336     // y = t x2 + u x + c

   337     // x = ( - u + Sqrt( u^2 - 4 t c ) ) / ( 2 t )

   338     disc = res.n[1]*res.n[1] - 4.0 * res.n[0] * res.n[2];

   339     if (disc < 0) return -1;

   340 

   341     return ( -1.0 * res.n[1] + sqrt( disc )) / (2.0 * res.n[0]);

   342 }

   343 

   344 static

   345 cmsBool IsMonotonic(int n, const cmsFloat64Number Table[])

   346 {

   347     int i;

   348     cmsFloat64Number last;

   349 

   350     last = Table[n-1];

   351 

   352     for (i = n-2; i >= 0; --i) {

   353 

   354         if (Table[i] > last)

   355 

   356             return FALSE;

   357         else

   358             last = Table[i];

   359 

   360     }

   361 

   362     return TRUE;

   363 }

   364 

   365 // Calculates the black point of a destination profile.

   366 // This algorithm comes from the Adobe paper disclosing its black point compensation method.

   367 cmsBool CMSEXPORT cmsDetectDestinationBlackPoint(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfile, cmsUInt32Number Intent, cmsUInt32Number dwFlags)

   368 {

   369     cmsColorSpaceSignature ColorSpace;

   370     cmsHTRANSFORM hRoundTrip = NULL;

   371     cmsCIELab InitialLab, destLab, Lab;

   372 

   373     cmsFloat64Number MinL, MaxL;

   374     cmsBool NearlyStraightMidRange = FALSE;

   375     cmsFloat64Number L;

   376     cmsFloat64Number x[101], y[101];

   377     cmsFloat64Number lo, hi, NonMonoMin;

   378     int n, l, i, NonMonoIndx;

   379 

   380 

   381     // Make sure intent is adequate

   382     if (Intent != INTENT_PERCEPTUAL &&

   383         Intent != INTENT_RELATIVE_COLORIMETRIC &&

   384         Intent != INTENT_SATURATION) {

   385         BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0;

   386         return FALSE;

   387     }

   388 

   389 

   390     // v4 + perceptual & saturation intents does have its own black point, and it is

   391     // well specified enough to use it. Black point tag is deprecated in V4.

   392     if ((cmsGetEncodedICCversion(hProfile) >= 0x4000000) &&

   393         (Intent == INTENT_PERCEPTUAL || Intent == INTENT_SATURATION)) {

   394 

   395             // Matrix shaper share MRC & perceptual intents

   396             if (cmsIsMatrixShaper(hProfile))

   397                 return BlackPointAsDarkerColorant(hProfile, INTENT_RELATIVE_COLORIMETRIC, BlackPoint, 0);

   398 

   399             // Get Perceptual black out of v4 profiles. That is fixed for perceptual & saturation intents

   400             BlackPoint -> X = cmsPERCEPTUAL_BLACK_X;

   401             BlackPoint -> Y = cmsPERCEPTUAL_BLACK_Y;

   402             BlackPoint -> Z = cmsPERCEPTUAL_BLACK_Z;

   403             return TRUE;

   404     }

   405 

   406 

   407     // Check if the profile is lut based and gray, rgb or cmyk (7.2 in Adobe's document)

   408     ColorSpace = cmsGetColorSpace(hProfile);

   409     if (!cmsIsCLUT(hProfile, Intent, LCMS_USED_AS_OUTPUT ) ||

   410         (ColorSpace != cmsSigGrayData &&

   411          ColorSpace != cmsSigRgbData  &&

   412          ColorSpace != cmsSigCmykData)) {

   413 

   414         // In this case, handle as input case

   415         return cmsDetectBlackPoint(BlackPoint, hProfile, Intent, dwFlags);

   416     }

   417 

   418     // It is one of the valid cases!, presto chargo hocus pocus, go for the Adobe magic

   419 

   420     // Step 1

   421     // ======

   422 

   423     // Set a first guess, that should work on good profiles.

   424     if (Intent == INTENT_RELATIVE_COLORIMETRIC) {

   425 

   426         cmsCIEXYZ IniXYZ;

   427 

   428         // calculate initial Lab as source black point

   429         if (!cmsDetectBlackPoint(&IniXYZ, hProfile, Intent, dwFlags)) {

   430             return FALSE;

   431         }

   432 

   433         // convert the XYZ to lab

   434         cmsXYZ2Lab(NULL, &InitialLab, &IniXYZ);

   435 

   436     } else {

   437 

   438         // set the initial Lab to zero, that should be the black point for perceptual and saturation

   439         InitialLab.L = 0;

   440         InitialLab.a = 0;

   441         InitialLab.b = 0;

   442     }

   443 

   444 

   445     // Step 2

   446     // ======

   447 

   448     // Create a roundtrip. Define a Transform BT for all x in L*a*b*

   449     hRoundTrip = CreateRoundtripXForm(hProfile, Intent);

   450     if (hRoundTrip == NULL)  return FALSE;

   451 

   452     // Calculate Min L*

   453     Lab = InitialLab;

   454     Lab.L = 0;

   455     cmsDoTransform(hRoundTrip, &Lab, &destLab, 1);

   456     MinL = destLab.L;

   457 

   458     // Calculate Max L*

   459     Lab = InitialLab;

   460     Lab.L = 100;

   461     cmsDoTransform(hRoundTrip, &Lab, &destLab, 1);

   462     MaxL = destLab.L;

   463 

   464     // Step 3

   465     // ======

   466 

   467     // check if quadratic estimation needs to be done.

   468     if (Intent == INTENT_RELATIVE_COLORIMETRIC) {

   469 

   470         // Conceptually, this code tests how close the source l and converted L are to one another in the mid-range

   471         // of the values. If the converted ramp of L values is close enough to a straight line y=x, then InitialLab

   472         // is good enough to be the DestinationBlackPoint,

   473         NearlyStraightMidRange = TRUE;

   474 

   475         for (l=0; l <= 100; l++) {

   476 

   477             Lab.L = l;

   478             Lab.a = InitialLab.a;

   479             Lab.b = InitialLab.b;

   480 

   481             cmsDoTransform(hRoundTrip, &Lab, &destLab, 1);

   482 

   483             L = destLab.L;

   484 

   485             // Check the mid range in 20% after MinL

   486             if (L > (MinL + 0.2 * (MaxL - MinL))) {

   487 

   488                 // Is close enough?

   489                 if (fabs(L - l) > 4.0) {

   490 

   491                     // Too far away, profile is buggy!

   492                     NearlyStraightMidRange = FALSE;

   493                     break;

   494                 }

   495             }

   496         }

   497     }

   498     else {

   499         // Check is always performed for perceptual and saturation intents

   500         NearlyStraightMidRange = FALSE;

   501     }

   502 

   503 

   504     // If no furter checking is needed, we are done

   505     if (NearlyStraightMidRange) {

   506 

   507         cmsLab2XYZ(NULL, BlackPoint, &InitialLab);

   508         cmsDeleteTransform(hRoundTrip);

   509         return TRUE;

   510     }

   511 

   512     // The round-trip curve normally looks like a nearly constant section at the black point,

   513     // with a corner and a nearly straight line to the white point.

   514 

   515     // STEP 4

   516     // =======

   517 

   518     // find the black point using the least squares error quadratic curve fitting

   519 

   520     if (Intent == INTENT_RELATIVE_COLORIMETRIC) {

   521         lo = 0.1;

   522         hi = 0.5;

   523     }

   524     else {

   525 

   526         // Perceptual and saturation

   527         lo = 0.03;

   528         hi = 0.25;

   529     }

   530 

   531     // Capture points for the fitting.

   532     n = 0;

   533     for (l=0; l <= 100; l++) {

   534 

   535         cmsFloat64Number ff;

   536 

   537         Lab.L = (cmsFloat64Number) l;

   538         Lab.a = InitialLab.a;

   539         Lab.b = InitialLab.b;

   540 

   541         cmsDoTransform(hRoundTrip, &Lab, &destLab, 1);

   542 

   543         ff = (destLab.L - MinL)/(MaxL - MinL);

   544 

   545         if (ff >= lo && ff < hi) {

   546 

   547             x[n] = Lab.L;

   548             y[n] = ff;

   549             n++;

   550         }

   551 

   552     }

   553 

   554     // This part is not on the Adobe paper, but I found is necessary for getting any result.

   555 

   556     if (IsMonotonic(n, y)) {

   557 

   558         // Monotonic means lower point is stil valid

   559         cmsLab2XYZ(NULL, BlackPoint, &InitialLab);

   560         cmsDeleteTransform(hRoundTrip);

   561         return TRUE;

   562     }

   563 

   564     // No suitable points, regret and use safer algorithm

   565     if (n == 0) {

   566         cmsDeleteTransform(hRoundTrip);

   567         return cmsDetectBlackPoint(BlackPoint, hProfile, Intent, dwFlags);

   568     }

   569 

   570 

   571     NonMonoMin = 100;

   572     NonMonoIndx = 0;

   573     for (i=0; i < n; i++) {

   574 

   575         if (y[i] < NonMonoMin) {

   576             NonMonoIndx = i;

   577             NonMonoMin = y[i];

   578         }

   579     }

   580 

   581     Lab.L = x[NonMonoIndx];

   582 

   583     // fit and get the vertex of quadratic curve

   584     Lab.L = RootOfLeastSquaresFitQuadraticCurve(n, x, y);

   585 

   586     if (Lab.L < 0.0 || Lab.L > 50.0) { // clip to zero L* if the vertex is negative

   587         Lab.L = 0;

   588     }

   589 

   590     Lab.a = InitialLab.a;

   591     Lab.b = InitialLab.b;

   592 

   593     cmsLab2XYZ(NULL, BlackPoint, &Lab);

   594 

   595     cmsDeleteTransform(hRoundTrip);

   596     return TRUE;

   597 }