1 //--------------------------------------------------------------------------------- 2 // 3 // Little Color Management System 4 // Copyright (c) 1998-2016 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 #define cmsmin(a, b) (((a) < (b)) ? (a) : (b)) 31 #define cmsmax(a, b) (((a) > (b)) ? (a) : (b)) 32 33 // This file contains routines for resampling and LUT optimization, black point detection 34 // and black preservation. 35 36 // Black point detection ------------------------------------------------------------------------- 37 38 39 // PCS -> PCS round trip transform, always uses relative intent on the device -> pcs 40 static 41 cmsHTRANSFORM CreateRoundtripXForm(cmsHPROFILE hProfile, cmsUInt32Number nIntent) 42 { 43 cmsContext ContextID = cmsGetProfileContextID(hProfile); 44 cmsHPROFILE hLab = cmsCreateLab4ProfileTHR(ContextID, NULL); 45 cmsHTRANSFORM xform; 46 cmsBool BPC[4] = { FALSE, FALSE, FALSE, FALSE }; 47 cmsFloat64Number States[4] = { 1.0, 1.0, 1.0, 1.0 }; 48 cmsHPROFILE hProfiles[4]; 49 cmsUInt32Number Intents[4]; 50 51 hProfiles[0] = hLab; hProfiles[1] = hProfile; hProfiles[2] = hProfile; hProfiles[3] = hLab; 52 Intents[0] = INTENT_RELATIVE_COLORIMETRIC; Intents[1] = nIntent; Intents[2] = INTENT_RELATIVE_COLORIMETRIC; Intents[3] = INTENT_RELATIVE_COLORIMETRIC; 53 54 xform = cmsCreateExtendedTransform(ContextID, 4, hProfiles, BPC, Intents, 55 States, NULL, 0, TYPE_Lab_DBL, TYPE_Lab_DBL, cmsFLAGS_NOCACHE|cmsFLAGS_NOOPTIMIZE); 56 57 cmsCloseProfile(hLab); 58 return xform; 59 } 60 61 // Use darker colorants to obtain black point. This works in the relative colorimetric intent and 62 // assumes more ink results in darker colors. No ink limit is assumed. 63 static 64 cmsBool BlackPointAsDarkerColorant(cmsHPROFILE hInput, 65 cmsUInt32Number Intent, 66 cmsCIEXYZ* BlackPoint, 67 cmsUInt32Number dwFlags) 68 { 69 cmsUInt16Number *Black; 70 cmsHTRANSFORM xform; 71 cmsColorSpaceSignature Space; 72 cmsUInt32Number nChannels; 73 cmsUInt32Number dwFormat; 74 cmsHPROFILE hLab; 75 cmsCIELab Lab; 76 cmsCIEXYZ BlackXYZ; 77 cmsContext ContextID = cmsGetProfileContextID(hInput); 78 79 // If the profile does not support input direction, assume Black point 0 80 if (!cmsIsIntentSupported(hInput, Intent, LCMS_USED_AS_INPUT)) { 81 82 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 83 return FALSE; 84 } 85 86 // Create a formatter which has n channels and floating point 87 dwFormat = cmsFormatterForColorspaceOfProfile(hInput, 2, FALSE); 88 89 // Try to get black by using black colorant 90 Space = cmsGetColorSpace(hInput); 91 92 // This function returns darker colorant in 16 bits for several spaces 93 if (!_cmsEndPointsBySpace(Space, NULL, &Black, &nChannels)) { 94 95 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 96 return FALSE; 97 } 98 99 if (nChannels != T_CHANNELS(dwFormat)) { 100 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 101 return FALSE; 102 } 103 104 // Lab will be used as the output space, but lab2 will avoid recursion 105 hLab = cmsCreateLab2ProfileTHR(ContextID, NULL); 106 if (hLab == NULL) { 107 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 108 return FALSE; 109 } 110 111 // Create the transform 112 xform = cmsCreateTransformTHR(ContextID, hInput, dwFormat, 113 hLab, TYPE_Lab_DBL, Intent, cmsFLAGS_NOOPTIMIZE|cmsFLAGS_NOCACHE); 114 cmsCloseProfile(hLab); 115 116 if (xform == NULL) { 117 118 // Something went wrong. Get rid of open resources and return zero as black 119 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 120 return FALSE; 121 } 122 123 // Convert black to Lab 124 cmsDoTransform(xform, Black, &Lab, 1); 125 126 // Force it to be neutral, clip to max. L* of 50 127 Lab.a = Lab.b = 0; 128 if (Lab.L > 50) Lab.L = 50; 129 130 // Free the resources 131 cmsDeleteTransform(xform); 132 133 // Convert from Lab (which is now clipped) to XYZ. 134 cmsLab2XYZ(NULL, &BlackXYZ, &Lab); 135 136 if (BlackPoint != NULL) 137 *BlackPoint = BlackXYZ; 138 139 return TRUE; 140 141 cmsUNUSED_PARAMETER(dwFlags); 142 } 143 144 // Get a black point of output CMYK profile, discounting any ink-limiting embedded 145 // in the profile. For doing that, we use perceptual intent in input direction: 146 // Lab (0, 0, 0) -> [Perceptual] Profile -> CMYK -> [Rel. colorimetric] Profile -> Lab 147 static 148 cmsBool BlackPointUsingPerceptualBlack(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfile) 149 { 150 cmsHTRANSFORM hRoundTrip; 151 cmsCIELab LabIn, LabOut; 152 cmsCIEXYZ BlackXYZ; 153 154 // Is the intent supported by the profile? 155 if (!cmsIsIntentSupported(hProfile, INTENT_PERCEPTUAL, LCMS_USED_AS_INPUT)) { 156 157 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 158 return TRUE; 159 } 160 161 hRoundTrip = CreateRoundtripXForm(hProfile, INTENT_PERCEPTUAL); 162 if (hRoundTrip == NULL) { 163 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 164 return FALSE; 165 } 166 167 LabIn.L = LabIn.a = LabIn.b = 0; 168 cmsDoTransform(hRoundTrip, &LabIn, &LabOut, 1); 169 170 // Clip Lab to reasonable limits 171 if (LabOut.L > 50) LabOut.L = 50; 172 LabOut.a = LabOut.b = 0; 173 174 cmsDeleteTransform(hRoundTrip); 175 176 // Convert it to XYZ 177 cmsLab2XYZ(NULL, &BlackXYZ, &LabOut); 178 179 if (BlackPoint != NULL) 180 *BlackPoint = BlackXYZ; 181 182 return TRUE; 183 } 184 185 // This function shouldn't exist at all -- there is such quantity of broken 186 // profiles on black point tag, that we must somehow fix chromaticity to 187 // avoid huge tint when doing Black point compensation. This function does 188 // just that. There is a special flag for using black point tag, but turned 189 // off by default because it is bogus on most profiles. The detection algorithm 190 // involves to turn BP to neutral and to use only L component. 191 cmsBool CMSEXPORT cmsDetectBlackPoint(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfile, cmsUInt32Number Intent, cmsUInt32Number dwFlags) 192 { 193 cmsProfileClassSignature devClass; 194 195 // Make sure the device class is adequate 196 devClass = cmsGetDeviceClass(hProfile); 197 if (devClass == cmsSigLinkClass || 198 devClass == cmsSigAbstractClass || 199 devClass == cmsSigNamedColorClass) { 200 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 201 return FALSE; 202 } 203 204 // Make sure intent is adequate 205 if (Intent != INTENT_PERCEPTUAL && 206 Intent != INTENT_RELATIVE_COLORIMETRIC && 207 Intent != INTENT_SATURATION) { 208 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 209 return FALSE; 210 } 211 212 // v4 + perceptual & saturation intents does have its own black point, and it is 213 // well specified enough to use it. Black point tag is deprecated in V4. 214 if ((cmsGetEncodedICCversion(hProfile) >= 0x4000000) && 215 (Intent == INTENT_PERCEPTUAL || Intent == INTENT_SATURATION)) { 216 217 // Matrix shaper share MRC & perceptual intents 218 if (cmsIsMatrixShaper(hProfile)) 219 return BlackPointAsDarkerColorant(hProfile, INTENT_RELATIVE_COLORIMETRIC, BlackPoint, 0); 220 221 // Get Perceptual black out of v4 profiles. That is fixed for perceptual & saturation intents 222 BlackPoint -> X = cmsPERCEPTUAL_BLACK_X; 223 BlackPoint -> Y = cmsPERCEPTUAL_BLACK_Y; 224 BlackPoint -> Z = cmsPERCEPTUAL_BLACK_Z; 225 226 return TRUE; 227 } 228 229 230 #ifdef CMS_USE_PROFILE_BLACK_POINT_TAG 231 232 // v2, v4 rel/abs colorimetric 233 if (cmsIsTag(hProfile, cmsSigMediaBlackPointTag) && 234 Intent == INTENT_RELATIVE_COLORIMETRIC) { 235 236 cmsCIEXYZ *BlackPtr, BlackXYZ, UntrustedBlackPoint, TrustedBlackPoint, MediaWhite; 237 cmsCIELab Lab; 238 239 // If black point is specified, then use it, 240 241 BlackPtr = cmsReadTag(hProfile, cmsSigMediaBlackPointTag); 242 if (BlackPtr != NULL) { 243 244 BlackXYZ = *BlackPtr; 245 _cmsReadMediaWhitePoint(&MediaWhite, hProfile); 246 247 // Black point is absolute XYZ, so adapt to D50 to get PCS value 248 cmsAdaptToIlluminant(&UntrustedBlackPoint, &MediaWhite, cmsD50_XYZ(), &BlackXYZ); 249 250 // Force a=b=0 to get rid of any chroma 251 cmsXYZ2Lab(NULL, &Lab, &UntrustedBlackPoint); 252 Lab.a = Lab.b = 0; 253 if (Lab.L > 50) Lab.L = 50; // Clip to L* <= 50 254 cmsLab2XYZ(NULL, &TrustedBlackPoint, &Lab); 255 256 if (BlackPoint != NULL) 257 *BlackPoint = TrustedBlackPoint; 258 259 return TRUE; 260 } 261 } 262 #endif 263 264 // That is about v2 profiles. 265 266 // If output profile, discount ink-limiting and that's all 267 if (Intent == INTENT_RELATIVE_COLORIMETRIC && 268 (cmsGetDeviceClass(hProfile) == cmsSigOutputClass) && 269 (cmsGetColorSpace(hProfile) == cmsSigCmykData)) 270 return BlackPointUsingPerceptualBlack(BlackPoint, hProfile); 271 272 // Nope, compute BP using current intent. 273 return BlackPointAsDarkerColorant(hProfile, Intent, BlackPoint, dwFlags); 274 } 275 276 277 278 // --------------------------------------------------------------------------------------------------------- 279 280 // Least Squares Fit of a Quadratic Curve to Data 281 // http://www.personal.psu.edu/jhm/f90/lectures/lsq2.html 282 283 static 284 cmsFloat64Number RootOfLeastSquaresFitQuadraticCurve(int n, cmsFloat64Number x[], cmsFloat64Number y[]) 285 { 286 double sum_x = 0, sum_x2 = 0, sum_x3 = 0, sum_x4 = 0; 287 double sum_y = 0, sum_yx = 0, sum_yx2 = 0; 288 double d, a, b, c; 289 int i; 290 cmsMAT3 m; 291 cmsVEC3 v, res; 292 293 if (n < 4) return 0; 294 295 for (i=0; i < n; i++) { 296 297 double xn = x[i]; 298 double yn = y[i]; 299 300 sum_x += xn; 301 sum_x2 += xn*xn; 302 sum_x3 += xn*xn*xn; 303 sum_x4 += xn*xn*xn*xn; 304 305 sum_y += yn; 306 sum_yx += yn*xn; 307 sum_yx2 += yn*xn*xn; 308 } 309 310 _cmsVEC3init(&m.v[0], n, sum_x, sum_x2); 311 _cmsVEC3init(&m.v[1], sum_x, sum_x2, sum_x3); 312 _cmsVEC3init(&m.v[2], sum_x2, sum_x3, sum_x4); 313 314 _cmsVEC3init(&v, sum_y, sum_yx, sum_yx2); 315 316 if (!_cmsMAT3solve(&res, &m, &v)) return 0; 317 318 319 a = res.n[2]; 320 b = res.n[1]; 321 c = res.n[0]; 322 323 if (fabs(a) < 1.0E-10) { 324 325 return cmsmin(0, cmsmax(50, -c/b )); 326 } 327 else { 328 329 d = b*b - 4.0 * a * c; 330 if (d <= 0) { 331 return 0; 332 } 333 else { 334 335 double rt = (-b + sqrt(d)) / (2.0 * a); 336 337 return cmsmax(0, cmsmin(50, rt)); 338 } 339 } 340 341 } 342 343 344 345 // Calculates the black point of a destination profile. 346 // This algorithm comes from the Adobe paper disclosing its black point compensation method. 347 cmsBool CMSEXPORT cmsDetectDestinationBlackPoint(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfile, cmsUInt32Number Intent, cmsUInt32Number dwFlags) 348 { 349 cmsColorSpaceSignature ColorSpace; 350 cmsHTRANSFORM hRoundTrip = NULL; 351 cmsCIELab InitialLab, destLab, Lab; 352 cmsFloat64Number inRamp[256], outRamp[256]; 353 cmsFloat64Number MinL, MaxL; 354 cmsBool NearlyStraightMidrange = TRUE; 355 cmsFloat64Number yRamp[256]; 356 cmsFloat64Number x[256], y[256]; 357 cmsFloat64Number lo, hi; 358 int n, l; 359 cmsProfileClassSignature devClass; 360 361 // Make sure the device class is adequate 362 devClass = cmsGetDeviceClass(hProfile); 363 if (devClass == cmsSigLinkClass || 364 devClass == cmsSigAbstractClass || 365 devClass == cmsSigNamedColorClass) { 366 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 367 return FALSE; 368 } 369 370 // Make sure intent is adequate 371 if (Intent != INTENT_PERCEPTUAL && 372 Intent != INTENT_RELATIVE_COLORIMETRIC && 373 Intent != INTENT_SATURATION) { 374 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 375 return FALSE; 376 } 377 378 379 // v4 + perceptual & saturation intents does have its own black point, and it is 380 // well specified enough to use it. Black point tag is deprecated in V4. 381 if ((cmsGetEncodedICCversion(hProfile) >= 0x4000000) && 382 (Intent == INTENT_PERCEPTUAL || Intent == INTENT_SATURATION)) { 383 384 // Matrix shaper share MRC & perceptual intents 385 if (cmsIsMatrixShaper(hProfile)) 386 return BlackPointAsDarkerColorant(hProfile, INTENT_RELATIVE_COLORIMETRIC, BlackPoint, 0); 387 388 // Get Perceptual black out of v4 profiles. That is fixed for perceptual & saturation intents 389 BlackPoint -> X = cmsPERCEPTUAL_BLACK_X; 390 BlackPoint -> Y = cmsPERCEPTUAL_BLACK_Y; 391 BlackPoint -> Z = cmsPERCEPTUAL_BLACK_Z; 392 return TRUE; 393 } 394 395 396 // Check if the profile is lut based and gray, rgb or cmyk (7.2 in Adobe's document) 397 ColorSpace = cmsGetColorSpace(hProfile); 398 if (!cmsIsCLUT(hProfile, Intent, LCMS_USED_AS_OUTPUT ) || 399 (ColorSpace != cmsSigGrayData && 400 ColorSpace != cmsSigRgbData && 401 ColorSpace != cmsSigCmykData)) { 402 403 // In this case, handle as input case 404 return cmsDetectBlackPoint(BlackPoint, hProfile, Intent, dwFlags); 405 } 406 407 // It is one of the valid cases!, use Adobe algorithm 408 409 410 // Set a first guess, that should work on good profiles. 411 if (Intent == INTENT_RELATIVE_COLORIMETRIC) { 412 413 cmsCIEXYZ IniXYZ; 414 415 // calculate initial Lab as source black point 416 if (!cmsDetectBlackPoint(&IniXYZ, hProfile, Intent, dwFlags)) { 417 return FALSE; 418 } 419 420 // convert the XYZ to lab 421 cmsXYZ2Lab(NULL, &InitialLab, &IniXYZ); 422 423 } else { 424 425 // set the initial Lab to zero, that should be the black point for perceptual and saturation 426 InitialLab.L = 0; 427 InitialLab.a = 0; 428 InitialLab.b = 0; 429 } 430 431 432 // Step 2 433 // ====== 434 435 // Create a roundtrip. Define a Transform BT for all x in L*a*b* 436 hRoundTrip = CreateRoundtripXForm(hProfile, Intent); 437 if (hRoundTrip == NULL) return FALSE; 438 439 // Compute ramps 440 441 for (l=0; l < 256; l++) { 442 443 Lab.L = (cmsFloat64Number) (l * 100.0) / 255.0; 444 Lab.a = cmsmin(50, cmsmax(-50, InitialLab.a)); 445 Lab.b = cmsmin(50, cmsmax(-50, InitialLab.b)); 446 447 cmsDoTransform(hRoundTrip, &Lab, &destLab, 1); 448 449 inRamp[l] = Lab.L; 450 outRamp[l] = destLab.L; 451 } 452 453 // Make monotonic 454 for (l = 254; l > 0; --l) { 455 outRamp[l] = cmsmin(outRamp[l], outRamp[l+1]); 456 } 457 458 // Check 459 if (! (outRamp[0] < outRamp[255])) { 460 461 cmsDeleteTransform(hRoundTrip); 462 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 463 return FALSE; 464 } 465 466 467 // Test for mid range straight (only on relative colorimetric) 468 NearlyStraightMidrange = TRUE; 469 MinL = outRamp[0]; MaxL = outRamp[255]; 470 if (Intent == INTENT_RELATIVE_COLORIMETRIC) { 471 472 for (l=0; l < 256; l++) { 473 474 if (! ((inRamp[l] <= MinL + 0.2 * (MaxL - MinL) ) || 475 (fabs(inRamp[l] - outRamp[l]) < 4.0 ))) 476 NearlyStraightMidrange = FALSE; 477 } 478 479 // If the mid range is straight (as determined above) then the 480 // DestinationBlackPoint shall be the same as initialLab. 481 // Otherwise, the DestinationBlackPoint shall be determined 482 // using curve fitting. 483 if (NearlyStraightMidrange) { 484 485 cmsLab2XYZ(NULL, BlackPoint, &InitialLab); 486 cmsDeleteTransform(hRoundTrip); 487 return TRUE; 488 } 489 } 490 491 492 // curve fitting: The round-trip curve normally looks like a nearly constant section at the black point, 493 // with a corner and a nearly straight line to the white point. 494 for (l=0; l < 256; l++) { 495 496 yRamp[l] = (outRamp[l] - MinL) / (MaxL - MinL); 497 } 498 499 // find the black point using the least squares error quadratic curve fitting 500 if (Intent == INTENT_RELATIVE_COLORIMETRIC) { 501 lo = 0.1; 502 hi = 0.5; 503 } 504 else { 505 506 // Perceptual and saturation 507 lo = 0.03; 508 hi = 0.25; 509 } 510 511 // Capture shadow points for the fitting. 512 n = 0; 513 for (l=0; l < 256; l++) { 514 515 cmsFloat64Number ff = yRamp[l]; 516 517 if (ff >= lo && ff < hi) { 518 x[n] = inRamp[l]; 519 y[n] = yRamp[l]; 520 n++; 521 } 522 } 523 524 525 // No suitable points 526 if (n < 3 ) { 527 cmsDeleteTransform(hRoundTrip); 528 BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; 529 return FALSE; 530 } 531 532 533 // fit and get the vertex of quadratic curve 534 Lab.L = RootOfLeastSquaresFitQuadraticCurve(n, x, y); 535 536 if (Lab.L < 0.0) { // clip to zero L* if the vertex is negative 537 Lab.L = 0; 538 } 539 540 Lab.a = InitialLab.a; 541 Lab.b = InitialLab.b; 542 543 cmsLab2XYZ(NULL, BlackPoint, &Lab); 544 545 cmsDeleteTransform(hRoundTrip); 546 return TRUE; 547 } 548