IccApplyBPC.cpp

Go to the documentation of this file.

/** @file
File:       IccApplyBPC.cpp
 
Contains:   Implementation of Black Point Compensation calculations.
 
Version:    V1
 
Copyright:  (c) see ICC Software License
*/
 
/*
* The ICC Software License, Version 0.2
*
*
* Copyright (c) 2003-2012 The International Color Consortium. All rights 
* reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
*
* 1. Redistributions of source code must retain the above copyright
*    notice, this list of conditions and the following disclaimer. 
*
* 2. Redistributions in binary form must reproduce the above copyright
*    notice, this list of conditions and the following disclaimer in
*    the documentation and/or other materials provided with the
*    distribution.
*
* 3. In the absence of prior written permission, the names "ICC" and "The
*    International Color Consortium" must not be used to imply that the
*    ICC organization endorses or promotes products derived from this
*    software.
*
*
* THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED.  IN NO EVENT SHALL THE INTERNATIONAL COLOR CONSORTIUM OR
* ITS CONTRIBUTING MEMBERS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
* USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
* ====================================================================
*
* This software consists of voluntary contributions made by many
* individuals on behalf of the The International Color Consortium. 
*
*
* Membership in the ICC is encouraged when this software is used for
* commercial purposes. 
*
*  
* For more information on The International Color Consortium, please
* see <http://www.color.org/>.
*  
* 
*/
 
////////////////////////////////////////////////////////////////////// 
// HISTORY:
//
// -Initial implementation by Rohit Patil 12-10-2008
//
//////////////////////////////////////////////////////////////////////
 
#include "IccApplyBPC.h"
#include <math.h>
 
#define IsSpacePCS(x) ((x)==icSigXYZData || (x)==icSigLabData)
 
/**
**************************************************************************
* Name: CIccApplyBPCHint::GetNewAdjustPCSXform
* 
* Purpose:
*  Returns a new CIccApplyBPC object. Returned object should be deleted 
*  by the caller.
* 
**************************************************************************
*/
IIccAdjustPCSXform* CIccApplyBPCHint::GetNewAdjustPCSXform() const
{
    return new CIccApplyBPC();
}
 
//////////////////////////////////////////////////////////////////////
// CIccApplyBPC utility functions
//////////////////////////////////////////////////////////////////////
 
// converts lab to pcs
void CIccApplyBPC::lab2pcs(icFloatNumber* pixel, const CIccProfile* pProfile) const
{
    switch (pProfile->m_Header.pcs)
    {
    case icSigLabData:
        icLabToPcs(pixel);
        break;
 
    case icSigXYZData:
        icLabtoXYZ(pixel);
        icXyzToPcs(pixel);
        break;
 
  default:
    break;
    }
}
 
// converts pcs to lab
void CIccApplyBPC::pcs2lab(icFloatNumber* pixel, const CIccProfile* pProfile) const
{
    switch (pProfile->m_Header.pcs)
    {
    case icSigLabData:
        icLabFromPcs(pixel);
        break;
 
    case icSigXYZData:
        icXyzFromPcs(pixel);
        icXYZtoLab(pixel);
        break;
 
  default:
    break;
    }
}
 
// calculates sum of product of x^j and y^k polynomials
icFloatNumber CIccApplyBPC::calcsum(icFloatNumber* x, icFloatNumber* y, int n, int j, int k) const
{
    icFloatNumber dSum = 0.0;
 
    int i;
    if (j && k) {
        for (i=0; i<n; i++) {
            dSum += (icFloatNumber)(pow(x[i], j)*pow(y[i], k));
        }
    }
    else if (j) {
        for (i=0; i<n; i++) {
            dSum += (icFloatNumber)pow(x[i], j);
        }
    }
    else if (k) {
        for (i=0; i<n; i++) {
            dSum += (icFloatNumber)pow(y[i], k);
        }
    }
    else {
        dSum = icFloatNumber(n);
    }
 
    return dSum;
}
 
// fits a quadratic curve through x,y points and returns the vertex of the parabola
icFloatNumber CIccApplyBPC::calcQuadraticVertex(icFloatNumber* x, icFloatNumber* y, int n) const
{
    icFloatNumber vert = 0.0;
    
    if (n>2) { // need at least three points to solve three linear equations
        icFloatNumber s00, s10, s20, s30, s40, s01, s11, s21, denom;
        s00 = calcsum(x, y, n, 0, 0);
        s10 = calcsum(x, y, n, 1, 0);
        s20 = calcsum(x, y, n, 2, 0);
        s30 = calcsum(x, y, n, 3, 0);
        s40 = calcsum(x, y, n, 4, 0);
        s01 = calcsum(x, y, n, 0, 1);
        s11 = calcsum(x, y, n, 1, 1);
        s21 = calcsum(x, y, n, 2, 1);
        denom = (icFloatNumber)(s00*s20*s40 - s10*s10*s40 - s00*s30*s30 + 2.0*s10*s20*s30 - s20*s20*s20);
        if (fabs(denom)>0.0) {
            // t and u are the coefficients of the quadratic equation y = tx^2 + ux + c
            // the three equations with 3 unknowns can be written as
            // [s40 s30 s20][t]   [s21]
            // [s30 s20 s10][u] = [s11]
            // [s20 s10 s00][c]   [s01]
            icFloatNumber t = (s01*s10*s30 - s11*s00*s30 - s01*s20*s20 + s11*s10*s20 + s21*s00*s20 - s21*s10*s10)/denom;
 
            icFloatNumber u = (s11*s00*s40 - s01*s10*s40 + s01*s20*s30 - s21*s00*s30 - s11*s20*s20 + s21*s10*s20)/denom;
 
            icFloatNumber c = (s01*s20*s40 - s11*s10*s40 - s01*s30*s30 + s11*s20*s30 + s21*s10*s30 - s21*s20*s20)/denom;
 
            // vertex is (-u + sqrt(u^2 - 4tc))/2t
            vert = (icFloatNumber)((-1.0 * u + sqrt(u*u - 4*t*c)) / (2.0 * t));
        }
    }
 
    return vert;
}
 
/**
**************************************************************************
* Name: CIccApplyBPC::CalculateFactors
* 
* Purpose:
*  This function does the suitable calculations to setup black point
*  compensation.
* 
* Args: 
*  pXform = pointer to the Xform object that calls this function
* 
* Return: 
*  true = all calculations done
*  false = an error occurred
**************************************************************************
*/
bool CIccApplyBPC::CalcFactors(const CIccProfile* pProfile, const CIccXform* pXform, icFloatNumber* Scale, icFloatNumber* Offset) const
{
    if (!pProfile || !pXform)
        return false;
 
    if (pXform->GetIntent()==icAbsoluteColorimetric) { // black point compensation not supported
        return false;
    }
 
    switch (pProfile->m_Header.deviceClass)
    { // These profile classes not supported
        case icSigLinkClass:
        case icSigAbstractClass:
        //case icSigColorSpaceClass:
        case icSigNamedColorClass:
            return false;
    default:
      break;
    }
 
    icFloatNumber XYZbp[3]; // storage for black point XYZ
 
    // calculate the black point
    if (!calcBlackPoint(pProfile, pXform, XYZbp)) {
        return false;
    }
 
    // calculate the scale and offset
    if (pXform->IsInput()) { // use PRM black as destination black
        Scale[0] = (icFloatNumber)((1.0 - icPerceptualRefBlackY)/(1.0 - XYZbp[1]));
    }
    else { // use PRM black as source black
        Scale[0] = (icFloatNumber)((1.0 - XYZbp[1])/(1.0 - icPerceptualRefBlackY));
    }
 
    Scale[1] = Scale[0];
    Scale[2] = Scale[0];
 
    Offset[0] = (icFloatNumber)((1.0 - Scale[0]) * icPerceptualRefWhiteX);
    Offset[1] = (icFloatNumber)((1.0 - Scale[1]) * icPerceptualRefWhiteY);
    Offset[2] = (icFloatNumber)((1.0 - Scale[2]) * icPerceptualRefWhiteZ);
 
    icXyzToPcs(Offset);
 
    return true;
}
 
/**
**************************************************************************
* Name: CIccApplyBPC::calcBlackPoint
* 
* Purpose:
*  Calculates the black point of a profile
* 
**************************************************************************
*/
bool CIccApplyBPC::calcBlackPoint(const CIccProfile* pProfile, const CIccXform* pXform, icFloatNumber* XYZb) const
{
    if (pXform->IsInput()) { // profile used as input/source profile
        return calcSrcBlackPoint(pProfile, pXform, XYZb);
    }
    else { // profile used as output profile
        return calcDstBlackPoint(pProfile, pXform, XYZb);
    }
 
    return true;
}
 
/**
**************************************************************************
* Name: CIccApplyBPC::calcSrcBlackPoint
* 
* Purpose:
*  Calculates the black point of a source profile
* 
**************************************************************************
*/
bool CIccApplyBPC::calcSrcBlackPoint(const CIccProfile* pProfile, const CIccXform* pXform, icFloatNumber* XYZb) const
{
    icFloatNumber Pixel[16];
    if ((pProfile->m_Header.colorSpace == icSigCmykData) && (pProfile->m_Header.deviceClass == icSigOutputClass)) {
 
        // calculate intermediate CMYK
        XYZb[0] = XYZb[1] = XYZb[2] = 0.0;
 
        // convert the Lab of 0,0,0 to relevant PCS
        lab2pcs(XYZb, pProfile);
 
        //convert the PCS value to CMYK
        if (!pixelXfm(Pixel, XYZb, pProfile->m_Header.pcs, icPerceptual, pProfile)) {
            return false;
        }
    }
    else {
        switch (pProfile->m_Header.colorSpace) {
                case icSigRgbData:
                    Pixel[0] = 0.0;
                    Pixel[1] = 0.0;
                    Pixel[2] = 0.0;
                    break;
 
                case icSigGrayData:
                    Pixel[0] = 0.0;
                    break;
 
                case icSigCmykData:
                case icSigCmyData:
                case icSig2colorData:
                case icSig3colorData:
                case icSig4colorData:
                case icSig5colorData:
                case icSig6colorData:
                case icSig7colorData:
                case icSig8colorData:
                case icSig9colorData:
                case icSig10colorData:
                case icSig11colorData:
                case icSig12colorData:
                case icSig13colorData:
                case icSig14colorData:
                case icSig15colorData:
                    {
                        icUInt32Number nSamples = icGetSpaceSamples(pProfile->m_Header.colorSpace);
                        for (icUInt32Number i=0; i<nSamples; i++) {
                            Pixel[i] = 1.0;
                        }
                    }
                    break;
 
                default:
                    return false;
        }
    }
 
    // convert the device value to PCS
    if (!pixelXfm(XYZb, Pixel, pProfile->m_Header.colorSpace, pXform->GetIntent(), pProfile)) {
        return false;
    }
 
    // convert PCS to Lab
    pcs2lab(XYZb, pProfile);
 
    // set a* b* to zero for cmyk profiles
    if (pProfile->m_Header.colorSpace == icSigCmykData) {
        XYZb[1] = XYZb[2] = 0.0;
    }
 
    // clip L* to 50
    if (XYZb[0]>50.0) {
        XYZb[0] = 50.0;
    }
 
    // convert Lab to XYZ
    icLabtoXYZ(XYZb);
    return true;
}
 
/**
**************************************************************************
* Name: CIccApplyBPC::calcDstBlackPoint
* 
* Purpose:
*  Calculates the black point of a destination profile
* 
**************************************************************************
*/
bool CIccApplyBPC::calcDstBlackPoint(const CIccProfile* pProfile, const CIccXform* pXform, icFloatNumber* XYZb) const
{
    icRenderingIntent nIntent = pXform->GetIntent();
    icFloatNumber Pixel[3];
    icFloatNumber pcsPixel[3];
 
    // check if the profile is lut based gray, rgb or cmyk
    if (pProfile->IsTagPresent(icSigBToA0Tag) && 
            (pProfile->m_Header.colorSpace==icSigGrayData || pProfile->m_Header.colorSpace==icSigRgbData || pProfile->m_Header.colorSpace==icSigCmykData))
    { // do the complicated and lengthy black point estimation
 
        // get the black transform
        CIccCmm* pCmm = getBlackXfm(nIntent, pProfile);
        if (!pCmm) {
            return false;
        }
 
        // set the initial Lab
        icFloatNumber iniLab[3] = {0.0, 0.0, 0.0};
 
        // calculate minL
        pcsPixel[0] = 0.0;
        pcsPixel[1] = iniLab[1];
        pcsPixel[2] = iniLab[2];
        lab2pcs(pcsPixel, pProfile);
        if (pCmm->Apply(Pixel, pcsPixel)!=icCmmStatOk) {
            delete pCmm;
            return false;
        }
        pcs2lab(Pixel, pProfile);
        icFloatNumber MinL = Pixel[0];
 
        // calculate MaxL
        pcsPixel[0] = 100.0;
        pcsPixel[1] = iniLab[1];
        pcsPixel[2] = iniLab[2];
        lab2pcs(pcsPixel, pProfile);
        if (pCmm->Apply(Pixel, pcsPixel)!=icCmmStatOk) {
            delete pCmm;
            return false;
        }
        pcs2lab(Pixel, pProfile);
        icFloatNumber MaxL = Pixel[0];
 
        // check if quadratic estimation needs to be done
        bool bStraightMidRange = false;
 
        // if the intent is relative
        if (nIntent==icRelativeColorimetric)
        { 
            // calculate initial Lab as source black point
            if (!calcSrcBlackPoint(pProfile, pXform, iniLab)) {
                delete pCmm;
                return false;
            }
 
            // convert the XYZ to lab
            icXYZtoLab(iniLab);
 
            // check mid range L* values
            icFloatNumber lcnt=0.0, roundtripL;
            bStraightMidRange = true;
            while (lcnt<100.1)
            {
                pcsPixel[0] = icFloatNumber(lcnt);
                pcsPixel[1] = iniLab[1];
                pcsPixel[2] = iniLab[2];
                lab2pcs(pcsPixel, pProfile);
                if (pCmm->Apply(Pixel, pcsPixel)!=icCmmStatOk) {
                    delete pCmm;
                    return false;
                }
                pcs2lab(Pixel, pProfile);
                roundtripL = Pixel[0];
 
                if (roundtripL>(MinL + 0.2 * (MaxL - MinL))) {
                    if (fabs(roundtripL - lcnt)>4.0) {
                        bStraightMidRange = false;
                        break;
                    }
                }
 
                lcnt += 1.0;
            }
        }
 
        // quadratic estimation is not needed
        if (bStraightMidRange) { // initial Lab is the destination black point
            XYZb[0] = iniLab[0];
            XYZb[1] = iniLab[1];
            XYZb[2] = iniLab[2];
            icLabtoXYZ(XYZb);
            delete pCmm;
            return true;
        }
 
        // find the black point using the least squares error quadratic curve fitting
 
        // calculate y values
        icFloatNumber x[101], y[101];
        icFloatNumber lo=0.03f, hi=0.25f;
        int i, n;
        if (nIntent==icRelativeColorimetric) {
            lo = 0.1f;
            hi = 0.5f;
        }
 
        for (i=0; i<101; i++) {
            x[i] = icFloatNumber(i);
            pcsPixel[0] = x[i];
            pcsPixel[1] = iniLab[1];
            pcsPixel[2] = iniLab[2];
            lab2pcs(pcsPixel, pProfile);
            if (pCmm->Apply(Pixel, pcsPixel)!=icCmmStatOk) {
                delete pCmm;
                return false;
            }
            pcs2lab(Pixel, pProfile);
            y[i] = (Pixel[0] - MinL)/(MaxL - MinL);
        }
 
        // check for y values in the range and rearrange
        n = 0;
        for (i=0; i<101; i++) {
            if (y[i]>=lo && y[i]<hi) {
                x[n] = x[i];
                y[n] = y[i];
                n++;
            }
        }
 
        if (!n) {
            delete pCmm;
            return false;
        }
 
        // fit and get the vertex of quadratic curve
        XYZb[0] = calcQuadraticVertex(x, y, n);
        if (XYZb[0]<0.0) { // clip to zero L* if the vertex is negative
            XYZb[0] = 0.0;
        }
        XYZb[1] = iniLab[1];
        XYZb[2] = iniLab[2];
        icLabtoXYZ(XYZb);
 
        delete pCmm;
    }
    else { // use the procedure for source black point
        return calcSrcBlackPoint(pProfile, pXform, XYZb);
    }
 
    return true;
}
 
/**
**************************************************************************
* Name: CIccApplyBPC::pixelXfm
* 
* Purpose:
*  Applies the specified transform to the source pixel
* 
**************************************************************************
*/
bool CIccApplyBPC::pixelXfm(icFloatNumber *DstPixel, icFloatNumber *SrcPixel, icColorSpaceSignature SrcSpace, 
                                                        icRenderingIntent nIntent, const CIccProfile *pProfile) const
{
    // create the cmm object
    CIccCmm cmm(SrcSpace, icSigUnknownData, !IsSpacePCS(SrcSpace));
 
    // first create a copy of the profile because the copy will be owned by the cmm
    CIccProfile* pICC = new CIccProfile(*pProfile);
    if (!pICC) return false;
 
    // add the xform
    if (cmm.AddXform(pICC, nIntent, icInterpTetrahedral, NULL, icXformLutColorimetric, pICC->m_Header.version >= icVersionNumberV5 ? false : true)!=icCmmStatOk) {
        delete pICC;
        return false;
    }
 
    // get the cmm ready to do transforms
    if (cmm.Begin()!=icCmmStatOk) {
        return false;
    }
 
    // Apply the pixel
    if (cmm.Apply(DstPixel, SrcPixel)!=icCmmStatOk) {
        return false;
    }
 
    return true;
}
 
/**
**************************************************************************
* Name: CIccApplyBPC::blackXfm
* 
* Purpose:
*  PCS -> PCS round trip transform, always uses relative intent on the device -> pcs transform
* 
**************************************************************************
*/
CIccCmm* CIccApplyBPC::getBlackXfm(icRenderingIntent nIntent, const CIccProfile *pProfile) const
{
    // create the cmm object
    CIccCmm* pCmm = new CIccCmm(pProfile->m_Header.pcs, icSigUnknownData, false);
    if (!pCmm) return NULL;
 
    // first create a copy of the profile because the copy will be owned by the cmm
    CIccProfile* pICC1 = new CIccProfile(*pProfile);
    if (!pICC1) {
        delete pCmm;
        return NULL;
    }
 
    // add the xform
    if (pCmm->AddXform(pICC1, nIntent, icInterpTetrahedral, NULL, icXformLutColor, pICC1->m_Header.version >= icVersionNumberV5 ? false : true)!=icCmmStatOk) {
        delete pICC1;
        delete pCmm;
        return NULL;
    }
 
    // create another copy of the profile because the copy will be owned by the cmm
    CIccProfile* pICC2 = new CIccProfile(*pProfile);
    if (!pICC2) {
        delete pCmm;
        return NULL;
    }
 
    // add the xform
    if (pCmm->AddXform(pICC2, icRelativeColorimetric, icInterpTetrahedral, NULL, icXformLutColor, pICC2->m_Header.version >= icVersionNumberV5 ? false : true)!=icCmmStatOk) { // uses the relative intent on the device to Lab side
        delete pICC2;
        delete pCmm;
        return NULL;
    }
 
    // get the cmm ready to do transforms
    if (pCmm->Begin()!=icCmmStatOk) {
        delete pCmm;
        return NULL;
    }
 
    return pCmm;
}