Added Adobe Black point detection algorithm

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';