From 839a5b09fdef1922091506185c12280e4343c9f9 Mon Sep 17 00:00:00 2001 From: Marti Maria Date: Fri, 23 Mar 2012 17:28:15 +0100 Subject: Added Adobe Black point detection algorithm --- ChangeLog | 3 +- include/lcms2.h | 1 + src/cmscnvrt.c | 2 +- src/cmsps2.c | 2 +- src/cmssamp.c | 304 +++++++++++++++++++++++++++++++++++++++++++++- src/cmstypes.c | 47 ++++--- src/lcms2.def | 3 +- utils/delphi/lcms2dll.pas | 3 + 8 files changed, 342 insertions(+), 23 deletions(-) diff --git a/ChangeLog b/ChangeLog index 55b40ea..394a3b8 100644 --- a/ChangeLog +++ b/ChangeLog @@ -60,4 +60,5 @@ Added compatibilty with Argyll's CGATS parser Fixed a bug in the named color devicelink generation Fixed uint64 to work in systems without long long native type Added performance improvements from several contributors, mostly artifex -Fixed a bug in black preservation checking \ No newline at end of file +Fixed a bug in black preservation checking +Added black point detection algorithm from Adobe paper diff --git a/include/lcms2.h b/include/lcms2.h index c2bbf9c..11965eb 100644 --- a/include/lcms2.h +++ b/include/lcms2.h @@ -1756,6 +1756,7 @@ CMSAPI cmsBool CMSEXPORT cmsGDBCheckPoint(cmsHANDLE hGBD, const cmsCIEL // Estimate the black point CMSAPI cmsBool CMSEXPORT cmsDetectBlackPoint(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfile, cmsUInt32Number Intent, cmsUInt32Number dwFlags); +CMSAPI cmsBool CMSEXPORT cmsDetectDestinationBlackPoint(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfile, cmsUInt32Number Intent, cmsUInt32Number dwFlags); // Estimate total area coverage CMSAPI cmsFloat64Number CMSEXPORT cmsDetectTAC(cmsHPROFILE hProfile); diff --git a/src/cmscnvrt.c b/src/cmscnvrt.c index 647a60e..c36a7d8 100644 --- a/src/cmscnvrt.c +++ b/src/cmscnvrt.c @@ -345,7 +345,7 @@ cmsBool ComputeConversion(int i, cmsHPROFILE hProfiles[], cmsCIEXYZ BlackPointIn, BlackPointOut; cmsDetectBlackPoint(&BlackPointIn, hProfiles[i-1], Intent, 0); - cmsDetectBlackPoint(&BlackPointOut, hProfiles[i], Intent, 0); + cmsDetectDestinationBlackPoint(&BlackPointOut, hProfiles[i], Intent, 0); // If black points are equal, then do nothing if (BlackPointIn.X != BlackPointOut.X || diff --git a/src/cmsps2.c b/src/cmsps2.c index 6b0a73c..5561926 100644 --- a/src/cmsps2.c +++ b/src/cmsps2.c @@ -948,7 +948,7 @@ cmsFloat64Number* GetPtrToMatrix(const cmsStage* mpe) // Does create CSA based on matrix-shaper. Allowed types are gray and RGB based static - int WriteInputMatrixShaper(cmsIOHANDLER* m, cmsHPROFILE hProfile, cmsStage* Matrix, cmsStage* Shaper) +int WriteInputMatrixShaper(cmsIOHANDLER* m, cmsHPROFILE hProfile, cmsStage* Matrix, cmsStage* Shaper) { cmsColorSpaceSignature ColorSpace; int rc; diff --git a/src/cmssamp.c b/src/cmssamp.c index 090d96d..a918c66 100644 --- a/src/cmssamp.c +++ b/src/cmssamp.c @@ -187,7 +187,6 @@ cmsBool BlackPointUsingPerceptualBlack(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfi // just that. There is a special flag for using black point tag, but turned // off by default because it is bogus on most profiles. The detection algorithm // involves to turn BP to neutral and to use only L component. - cmsBool CMSEXPORT cmsDetectBlackPoint(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfile, cmsUInt32Number Intent, cmsUInt32Number dwFlags) { @@ -264,3 +263,306 @@ cmsBool CMSEXPORT cmsDetectBlackPoint(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfil } + +// --------------------------------------------------------------------------------------------------------- + +// Least Squares Fit of a Quadratic Curve to Data +// http://www.personal.psu.edu/jhm/f90/lectures/lsq2.html + +static +cmsFloat64Number VertexOfLeastSquaresFitQuadraticCurve(int n, cmsFloat64Number x[], cmsFloat64Number y[]) +{ + double sum_x = 0, sum_x2 = 0, sum_x3 = 0, sum_x4 = 0; + double sum_y = 0, sum_yx = 0, sum_yx2 = 0; + double disc; + int i; + cmsMAT3 m; + cmsVEC3 v, res; + + if (n < 4) return 0; + + for (i=0; i < n; i++) { + + double xn = x[i]; + double yn = y[i]; + + sum_x += xn; + sum_x2 += xn*xn; + sum_x3 += xn*xn*xn; + sum_x4 += xn*xn*xn*xn; + + sum_y += yn; + sum_yx += yn*xn; + sum_yx2 += yn*xn*xn; + } + + _cmsVEC3init(&m.v[0], n, sum_x, sum_x2); + _cmsVEC3init(&m.v[1], sum_x, sum_x2, sum_x3); + _cmsVEC3init(&m.v[2], sum_x2, sum_x3, sum_x4); + + _cmsVEC3init(&v, sum_y, sum_yx, sum_yx2); + + if (!_cmsMAT3solve(&res, &m, &v)) return 0; + + // y = t x2 + u x + c + // x = ( - u + Sqrt( u^2 - 4 t c ) ) / ( 2 t ) + disc = res.n[1]*res.n[1] - 4.0 * res.n[0] * res.n[2]; + if (disc < 0) return -1; + + return ( -1.0 * res.n[1] + sqrt( disc )) / (2.0 * res.n[0]); +} + +static +cmsBool IsMonotonic(int n, const cmsFloat64Number Table[]) +{ + int i; + cmsFloat64Number last; + + last = Table[n-1]; + + for (i = n-2; i >= 0; --i) { + + if (Table[i] > last) + + return FALSE; + else + last = Table[i]; + + } + + return TRUE; +} + +// Calculates the black point of a destination profile. +// This algorithm comes from the Adobe paper disclosing its black point compensation method. +cmsBool CMSEXPORT cmsDetectDestinationBlackPoint(cmsCIEXYZ* BlackPoint, cmsHPROFILE hProfile, cmsUInt32Number Intent, cmsUInt32Number dwFlags) +{ + cmsColorSpaceSignature ColorSpace; + cmsHTRANSFORM hRoundTrip = NULL; + cmsCIELab InitialLab, destLab, Lab; + + cmsFloat64Number MinL, MaxL; + cmsBool NearlyStraightMidRange = FALSE; + cmsFloat64Number L; + cmsFloat64Number x[101], y[101]; + cmsFloat64Number lo, hi, NonMonoMin; + int n, l, i, NonMonoIndx; + + + // Make sure intent is adequate + if (Intent != INTENT_PERCEPTUAL && + Intent != INTENT_RELATIVE_COLORIMETRIC && + Intent != INTENT_SATURATION) { + BlackPoint -> X = BlackPoint ->Y = BlackPoint -> Z = 0.0; + return FALSE; + } + + + // v4 + perceptual & saturation intents does have its own black point, and it is + // well specified enough to use it. Black point tag is deprecated in V4. + if ((cmsGetEncodedICCversion(hProfile) >= 0x4000000) && + (Intent == INTENT_PERCEPTUAL || Intent == INTENT_SATURATION)) { + + // Matrix shaper share MRC & perceptual intents + if (cmsIsMatrixShaper(hProfile)) + return BlackPointAsDarkerColorant(hProfile, INTENT_RELATIVE_COLORIMETRIC, BlackPoint, 0); + + // Get Perceptual black out of v4 profiles. That is fixed for perceptual & saturation intents + BlackPoint -> X = cmsPERCEPTUAL_BLACK_X; + BlackPoint -> Y = cmsPERCEPTUAL_BLACK_Y; + BlackPoint -> Z = cmsPERCEPTUAL_BLACK_Z; + return TRUE; + } + + + // Check if the profile is lut based and gray, rgb or cmyk (7.2 in Adobe's document) + ColorSpace = cmsGetColorSpace(hProfile); + if (!cmsIsCLUT(hProfile, Intent, LCMS_USED_AS_OUTPUT ) || + (ColorSpace != cmsSigGrayData && + ColorSpace != cmsSigRgbData && + ColorSpace != cmsSigCmykData)) { + + // In this case, handle as input case + return cmsDetectBlackPoint(BlackPoint, hProfile, Intent, dwFlags); + } + + // It is one of the valid cases!, presto chargo hocus pocus, go for the Adobe magic + + // Step 1 + // ====== + + // Set a first guess, that should work on good profiles. + if (Intent == INTENT_RELATIVE_COLORIMETRIC) { + + cmsCIEXYZ IniXYZ; + + // calculate initial Lab as source black point + if (!cmsDetectBlackPoint(&IniXYZ, hProfile, Intent, dwFlags)) { + return FALSE; + } + + // convert the XYZ to lab + cmsXYZ2Lab(NULL, &InitialLab, &IniXYZ); + + } else { + + // set the initial Lab to zero, that should be the black point for perceptual and saturation + InitialLab.L = 0; + InitialLab.a = 0; + InitialLab.b = 0; + } + + + // Step 2 + // ====== + + // Create a roundtrip. Define a Transform BT for all x in L*a*b* + hRoundTrip = CreateRoundtripXForm(hProfile, Intent); + if (hRoundTrip == NULL) return FALSE; + + // Calculate Min L* + Lab = InitialLab; + Lab.L = 0; + cmsDoTransform(hRoundTrip, &Lab, &destLab, 1); + MinL = destLab.L; + + // Calculate Max L* + Lab = InitialLab; + Lab.L = 100; + cmsDoTransform(hRoundTrip, &Lab, &destLab, 1); + MaxL = destLab.L; + + // Step 3 + // ====== + + // check if quadratic estimation needs to be done. + if (Intent == INTENT_RELATIVE_COLORIMETRIC) { + + // Conceptually, this code tests how close the source l and converted L are to one another in the mid-range + // of the values. If the converted ramp of L values is close enough to a straight line y=x, then InitialLab + // is good enough to be the DestinationBlackPoint, + NearlyStraightMidRange = TRUE; + + for (l=0; l <= 100; l++) { + + Lab.L = l; + Lab.a = InitialLab.a; + Lab.b = InitialLab.b; + + cmsDoTransform(hRoundTrip, &Lab, &destLab, 1); + + L = destLab.L; + + // Check the mid range in 20% after MinL + if (L > (MinL + 0.2 * (MaxL - MinL))) { + + // Is close enough? + if (fabs(L - l) > 4.0) { + + // Too far away, profile is buggy! + NearlyStraightMidRange = FALSE; + break; + } + } + } + } + else { + // Check is always performed for perceptual and saturation intents + NearlyStraightMidRange = FALSE; + } + + + // If no furter checking is needed, we are done + if (NearlyStraightMidRange) { + + cmsLab2XYZ(NULL, BlackPoint, &InitialLab); + cmsDeleteTransform(hRoundTrip); + return TRUE; + } + + // The round-trip curve normally looks like a nearly constant section at the black point, + // with a corner and a nearly straight line to the white point. + + // STEP 4 + // ======= + + // find the black point using the least squares error quadratic curve fitting + + if (Intent == INTENT_RELATIVE_COLORIMETRIC) { + lo = 0.1; + hi = 0.5; + } + else { + + // Perceptual and saturation + lo = 0.03; + hi = 0.25; + } + + // Capture points for the fitting. + n = 0; + for (l=0; l <= 100; l++) { + + cmsFloat64Number ff; + + Lab.L = (cmsFloat64Number) l; + Lab.a = InitialLab.a; + Lab.b = InitialLab.b; + + cmsDoTransform(hRoundTrip, &Lab, &destLab, 1); + + ff = (destLab.L - MinL)/(MaxL - MinL); + + if (ff >= lo && ff < hi) { + + x[n] = Lab.L; + y[n] = ff; + n++; + } + + } + + // This part is not on the Adobe paper, but I found is necessary for getting any result. + + if (IsMonotonic(n, y)) { + + // Monotonic means lower point is stil valid + cmsLab2XYZ(NULL, BlackPoint, &InitialLab); + cmsDeleteTransform(hRoundTrip); + return TRUE; + } + + // No suitable points, regret and use safer algorithm + if (n == 0) { + cmsDeleteTransform(hRoundTrip); + return cmsDetectBlackPoint(BlackPoint, hProfile, Intent, dwFlags); + } + + + NonMonoMin = 100; + NonMonoIndx = 0; + for (i=0; i < n; i++) { + + if (y[i] < NonMonoMin) { + NonMonoIndx = i; + NonMonoMin = y[i]; + } + } + + Lab.L = x[NonMonoIndx]; + + // fit and get the vertex of quadratic curve + Lab.L = VertexOfLeastSquaresFitQuadraticCurve(n, x, y); + + if (Lab.L < 0.0 || Lab.L > 50.0) { // clip to zero L* if the vertex is negative + Lab.L = 0; + } + + Lab.a = InitialLab.a; + Lab.b = InitialLab.b; + + cmsLab2XYZ(NULL, BlackPoint, &Lab); + + cmsDeleteTransform(hRoundTrip); + return TRUE; +} diff --git a/src/cmstypes.c b/src/cmstypes.c index 1df9bd6..7a8c5c8 100644 --- a/src/cmstypes.c +++ b/src/cmstypes.c @@ -125,7 +125,7 @@ cmsBool _cmsWriteWCharArray(cmsIOHANDLER* io, cmsUInt32Number n, const wchar_t* cmsUInt32Number i; _cmsAssert(io != NULL); - _cmsAssert(Array != NULL); + _cmsAssert(!(Array == NULL && n > 0)); for (i=0; i < n; i++) { if (!_cmsWriteUInt16Number(io, (cmsUInt16Number) Array[i])) return FALSE; @@ -1464,13 +1464,19 @@ void *Type_MLU_Read(struct _cms_typehandler_struct* self, cmsIOHANDLER* io, cmsU // Now read the remaining of tag and fill all strings. Substract the directory SizeOfTag = (LargestPosition * sizeof(wchar_t)) / sizeof(cmsUInt16Number); + if (SizeOfTag == 0) + { + Block = NULL; + NumOfWchar = 0; - Block = (wchar_t*) _cmsMalloc(self ->ContextID, SizeOfTag); - if (Block == NULL) goto Error; - - NumOfWchar = SizeOfTag / sizeof(wchar_t); - - if (!_cmsReadWCharArray(io, NumOfWchar, Block)) goto Error; + } + else + { + Block = (wchar_t*) _cmsMalloc(self ->ContextID, SizeOfTag); + if (Block == NULL) goto Error; + NumOfWchar = SizeOfTag / sizeof(wchar_t); + if (!_cmsReadWCharArray(io, NumOfWchar, Block)) goto Error; + } mlu ->MemPool = Block; mlu ->PoolSize = SizeOfTag; @@ -2535,7 +2541,8 @@ cmsBool WriteSetOfCurves(struct _cms_typehandler_struct* self, cmsIOHANDLER* io, // If this is a table-based curve, use curve type even on V4 CurrentType = Type; - if ((Curves[i] ->nSegments == 0)||((Curves[i]->nSegments == 2) && (Curves[i] ->Segments[1].Type == 0))) + if ((Curves[i] ->nSegments == 0)|| + ((Curves[i]->nSegments == 2) && (Curves[i] ->Segments[1].Type == 0)) ) CurrentType = cmsSigCurveType; else if (Curves[i] ->Segments[0].Type < 0) @@ -3217,28 +3224,32 @@ void *Type_ProfileSequenceDesc_Read(struct _cms_typehandler_struct* self, cmsIOH cmsPSEQDESC* sec = &OutSeq -> seq[i]; - if (!_cmsReadUInt32Number(io, &sec ->deviceMfg)) return NULL; - if (SizeOfTag < sizeof(cmsUInt32Number)) return NULL; + if (!_cmsReadUInt32Number(io, &sec ->deviceMfg)) goto Error; + if (SizeOfTag < sizeof(cmsUInt32Number)) goto Error; SizeOfTag -= sizeof(cmsUInt32Number); - if (!_cmsReadUInt32Number(io, &sec ->deviceModel)) return NULL; - if (SizeOfTag < sizeof(cmsUInt32Number)) return NULL; + if (!_cmsReadUInt32Number(io, &sec ->deviceModel)) goto Error; + if (SizeOfTag < sizeof(cmsUInt32Number)) goto Error; SizeOfTag -= sizeof(cmsUInt32Number); - if (!_cmsReadUInt64Number(io, &sec ->attributes)) return NULL; - if (SizeOfTag < sizeof(cmsUInt32Number)) return NULL; + if (!_cmsReadUInt64Number(io, &sec ->attributes)) goto Error; + if (SizeOfTag < sizeof(cmsUInt32Number)) goto Error; SizeOfTag -= sizeof(cmsUInt64Number); - if (!_cmsReadUInt32Number(io, (cmsUInt32Number *)&sec ->technology)) return NULL; - if (SizeOfTag < sizeof(cmsUInt32Number)) return NULL; + if (!_cmsReadUInt32Number(io, (cmsUInt32Number *)&sec ->technology)) goto Error; + if (SizeOfTag < sizeof(cmsUInt32Number)) goto Error; SizeOfTag -= sizeof(cmsUInt32Number); - if (!ReadEmbeddedText(self, io, &sec ->Manufacturer, SizeOfTag)) return NULL; - if (!ReadEmbeddedText(self, io, &sec ->Model, SizeOfTag)) return NULL; + if (!ReadEmbeddedText(self, io, &sec ->Manufacturer, SizeOfTag)) goto Error; + if (!ReadEmbeddedText(self, io, &sec ->Model, SizeOfTag)) goto Error; } *nItems = 1; return OutSeq; + +Error: + cmsFreeProfileSequenceDescription(OutSeq); + return NULL; } diff --git a/src/lcms2.def b/src/lcms2.def index f4d1371..64240fe 100644 --- a/src/lcms2.def +++ b/src/lcms2.def @@ -62,7 +62,8 @@ _cmsDecodeDateTimeNumber = _cmsDecodeDateTimeNumber _cmsDefaultICCintents = _cmsDefaultICCintents cmsDeleteTransform = cmsDeleteTransform cmsDeltaE = cmsDeltaE -cmsDetectBlackPoint = cmsDetectBlackPoint +cmsDetectBlackPoint = cmsDetectBlackPoint +cmsDetectDestinationBlackPoint = cmsDetectDestinationBlackPoint cmsDetectTAC = cmsDetectTAC cmsDesaturateLab = cmsDesaturateLab cmsDoTransform = cmsDoTransform diff --git a/utils/delphi/lcms2dll.pas b/utils/delphi/lcms2dll.pas index e632631..19e3f87 100755 --- a/utils/delphi/lcms2dll.pas +++ b/utils/delphi/lcms2dll.pas @@ -1615,6 +1615,8 @@ FUNCTION cmsGDBCheckPoint(hGBD: cmsHANDLE; Lab: LPcmsCIELab): cmsBool; StdCall; // Estimate the black point FUNCTION cmsDetectBlackPoint( BlackPoint: LPcmsCIEXYZ; hProfile: cmsHPROFILE; Intent: cmsUInt32Number; dwFlags: cmsUInt32Number): cmsBool; StdCall; +FUNCTION cmsDetectDestinationBlackPoint( BlackPoint: LPcmsCIEXYZ; hProfile: cmsHPROFILE; Intent: cmsUInt32Number; dwFlags: cmsUInt32Number): cmsBool; StdCall; + // Estimate total area coverage FUNCTION cmsDetectTAC(hProfile: cmsHPROFILE): cmsFloat64Number; StdCall; @@ -2100,6 +2102,7 @@ FUNCTION cmsGDBCompute(hGDB: cmsHANDLE; dwFlags: cmsUInt32Number): cmsBool; Std FUNCTION cmsGDBCheckPoint(hGBD: cmsHANDLE; Lab: LPcmsCIELab): cmsBool; StdCall; external 'lcms2.dll'; FUNCTION cmsDetectBlackPoint( BlackPoint: LPcmsCIEXYZ; hProfile: cmsHPROFILE; Intent: cmsUInt32Number; dwFlags: cmsUInt32Number): cmsBool; StdCall; external 'lcms2.dll'; +FUNCTION cmsDetectDestinationBlackPoint( BlackPoint: LPcmsCIEXYZ; hProfile: cmsHPROFILE; Intent: cmsUInt32Number; dwFlags: cmsUInt32Number): cmsBool; StdCall; external 'lcms2.dll'; FUNCTION cmsDetectTAC(hProfile: cmsHPROFILE): cmsFloat64Number; StdCall; external 'lcms2.dll'; -- cgit v1.2.1