path: root/src/libs/lprof/cmsreg.cpp

/* */
/*  Little cms - profiler construction set */
/*  Copyright (C) 1998-2001 Marti Maria <[email protected]> */
/* */
/* THIS SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, */
/* EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY */
/* WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. */
/* */
/* IN NO EVENT SHALL MARTI MARIA BE LIABLE FOR ANY SPECIAL, INCIDENTAL, */
/* INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
/* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, */
/* WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF */
/* LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE */
/* OF THIS SOFTWARE. */
/* */
/* This file is free software; you can redistribute it and/or modify it */
/* under the terms of the GNU General Public License as published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, but */
/* WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU */
/* General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software */
/* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */
/* */
/* As a special exception to the GNU General Public License, if you */
/* distribute this file as part of a program that contains a */
/* configuration script generated by Autoconf, you may include it under */
/* the same distribution terms that you use for the rest of that program. */
/* */
/* Version 1.09a */


#include "lcmsprf.h"


/* There are three kinds of lies: */
/* */
/*   * lies */
/*   * damn lies */
/*   * statistics */
/* */
/*                    -Some Wag */
/* */
/* */
/* This module handles multiple linear regression stuff */



/* A measurement of error

typedef struct {

        double SSE;             // The error sum of squares
        double MSE;             // The error mean sum of squares
        double SSR;             // The regression sum of squares
        double MSR;             // The regression mean sum of squares
        double SSTO;            // Total sum of squares
        double F;               // The Fisher-F value (MSR / MSE)
        double R2;              // Proportion of variability explained by the regression
                                // (root is Pearson correlation coefficient)

        double R2adj;           // The adjusted coefficient of multiple determination.
                                // R2-adjusted or R2adj. This is calculated as
                                // R2adj = 1 - (1-R2)(N-n-1)/(N-1)
                                // and used as multiple correlation coefficient
                                // (really, it should be square root)

    } MLRSTATISTICS, FAR* LPMLRSTATISTICS;

*/


int  cdecl cmsxRegressionCreateMatrix(LPMEASUREMENT m, SETOFPATCHES Allowed, int nterms,
                                       int ColorSpace,
                                       LPMATN* lpMat, LPMLRSTATISTICS Stat);

BOOL cdecl cmsxRegressionRGB2Lab(double r, double g, double b,
                                       LPMATN tfm, LPcmsCIELab Lab);

BOOL cdecl cmsxRegressionRGB2XYZ(double r, double g, double b,
                                       LPMATN tfm, LPcmsCIEXYZ XYZ);


/* -------------------------------------------------------------- Implementation */

/* #define DEBUG 1 */


/* Multiple linear regression. Also keep track of error. */
/* Returns false if something goes wrong, or true if all Ok. */

static
BOOL MultipleLinearRegression(const LPMATN xi,                /* Dependent variable */
                              const LPMATN y,                 /* Independent variable */
                              int nvar,                       /* Number of samples */
                              int npar,                       /* Number of parameters (terms) */
                              double* coeff,                  /* Returned coefficients */
                              LPMATN vcb,                     /* Variance-covariance array */
                              double *tvl,                    /* T-Values */
                              LPMLRSTATISTICS ans)            /* The returned statistics */
{
    LPMATN bt, xt, a, xy, yt, b;
    double  sum;
    LPMATN temp1, temp2;
    int i;


    /* |xt| = |xi| T */
    xt = MATNtranspose(xi);
    if (xt == NULL) return false;


    /* |a| = |xt|* |xi| */
    a = MATNmult(xt, xi);
    if (a == NULL) return false;


    /* |xy| = |xy| * |y| */
    xy = MATNmult (xt, y);
    if (xy == NULL) return false;


    /* solve system |a|*|xy| = 0 */
    if (!MATNsolve(a, xy)) return false;

    /* b will hold coefficients */
    b = MATNalloc (xy->Rows, 1);
    if (b == NULL) return false;

    for (i = 0; i < npar; i++)
        b->Values[i][0] = xy->Values[i][0];

    /* Store a copy for later user */
    for (i = 0; i < npar; i++)
        coeff[i] = b->Values[i][0];

    /* Error analysis. */

    /* SSE and MSE. */
    temp1 = MATNalloc (1,1);
    if ((temp1->Values[0][0] = MATNcross(y)) == 0) return false;

    /* |bt| = |b| T */
    bt = MATNtranspose (b);
    if (bt == NULL) return false;

    /*  |yt| = |bt| * |xt| */
    yt = MATNmult (bt, xt);
    if (yt == NULL) return false;


    /* |temp2| = |yt|* |y| */
    temp2 = MATNmult (yt, y);
    if (temp2 == NULL) return false;

    /* SSE, MSE */
    ans->SSE = temp1 -> Values[0][0] - temp2 -> Values[0][0];
    ans->MSE = ans->SSE / (double) (nvar - npar);

    /* SSTO */
    sum = 0;
    for (i=0; i < nvar; i++)
                sum += y->Values[i][0];

    sum *= sum / (double) nvar;
    ans->SSTO = temp1->Values[0][0] - sum;

    /* SSR, MSR, and Fisher-F */
    ans->SSR = temp2->Values[0][0] - sum;
    ans->MSR = ans->SSR / (double) (npar - 1);
    ans->F   = ans->MSR / ans->MSE;

    /* Correlation coefficients. */
    ans->R2 = ans->SSR/ans->SSTO;
    ans->R2adj = 1.0 - (ans->SSE/ans->SSTO)*((nvar-1.)/(nvar-npar));

    /* Variance-covariance matrix */
    /* */
    /*     In RGB->Lab, for example: */
    /* */
    /*          Var(R)    Cov(R,G)  Cov(R,B) */
    /* |vcb| =  Cov(R,G)  Var(G)    Cov(G,B) */
    /*          Cov(R,B)  Cov(G,B)  Var(B) */
    /* */

    MATNscalar(a, ans->MSE, vcb);

    /* Determine the T-values */

    for (i=0; i < npar; i++) {

        temp1->Values[0][0] = fabs(vcb->Values[i][0]);
        if ( temp1->Values[0][0] == 0)
            tvl[i] = 0;             /* This should never happen */
        else
            tvl[i] = b->Values[i][0] / sqrt(temp1->Values[0][0]);
    }


    /* Ok, done */

    MATNfree(a); MATNfree(xy); MATNfree(yt); MATNfree(b);
    MATNfree(temp1); MATNfree(temp2); MATNfree(bt); MATNfree(xt);


    return true;
}



/* Does create (so, it allocates) the regression matrix, */
/* keeping track of error as well. */

static
BOOL CreateRegressionMatrix(const LPMATN Input, const LPMATN Output,
                            LPMATN* ptrMatrix, LPMLRSTATISTICS maxErrorMeas)
{
    double* coef;
    double* tval;
    LPMATN ivar, dvar, vcov;
    MLRSTATISTICS ErrorMeas, PeakErrorMeas;
    int i, j, nIn, nOut, NumOfPatches;

        nIn          = Input  -> Cols;
        nOut         = Output -> Cols;
        NumOfPatches = Input  -> Rows;

        /* Checkpoint */
        if (Output -> Rows != NumOfPatches) {

                cmsSignalError(LCMS_ERRC_ABORTED, "(internal) Regression matrix mismatch");
                return false;
        }

    coef = (double*) malloc(nIn * sizeof(double));
    if (coef == NULL) return false;

    tval = (double*) malloc(nIn * sizeof(double));
    if (tval == NULL) {
        free(coef);
        return false;
    }

    ivar = MATNalloc(NumOfPatches, nIn);
    dvar = MATNalloc(NumOfPatches, 1);

    /* Copy In to ivar, */
    for (i = 0; i < NumOfPatches; i++) {

        for (j = 0; j < nIn; j++)
            ivar->Values[i][j] = Input->Values[i][j];
    }

    /* This is the (symmetric) Covariance matrix */
    vcov = MATNalloc(nIn, nIn);

    /* This is the regression matrix */
    *ptrMatrix = MATNalloc(nIn, nOut);

    PeakErrorMeas.R2adj = 0;
    for (j = 0; j < nOut; ++j)
    {
        for (i = 0; i < NumOfPatches; ++i)
            dvar->Values[i][0] = Output->Values[i][j];

        if (MultipleLinearRegression(ivar, dvar, NumOfPatches, nIn, coef, vcov, tval, &ErrorMeas)) {

            /* Ok so far... store values */
            for (i = 0; i < nIn; i++)
                (*ptrMatrix)->Values[i][j] = coef[i];
        }
        else {
                    /* Boo... got error. Discard whole point. */
                    MATNfree(ivar); MATNfree(dvar); MATNfree(vcov);
                    if (coef) free(coef);
                    if (tval) free(tval);
                    MATNfree(*ptrMatrix); *ptrMatrix = NULL;
                    return false;
        }

        /* Did this colorant got higer error? If so, this is  */
        /* the peak of all pixel */

        if(fabs(ErrorMeas.R2adj) > fabs(PeakErrorMeas.R2adj))
            PeakErrorMeas = ErrorMeas;
    }

    /* This is the peak error on all components */
    *maxErrorMeas = PeakErrorMeas;


#ifdef DEBUG
        MATNprintf("Variance-Covariance", vcov);
        printf("R2adj: %g, F: %g\n", PeakErrorMeas.R2adj, PeakErrorMeas.F);
#endif

    /* Free stuff. */
    MATNfree(ivar); MATNfree(dvar); MATNfree(vcov);
    if (coef) free(coef);
    if (tval) free(tval);

    return true;
}


/* Does compute the term of regression based on inputs. */

static
double Term(int n, double r, double g, double b)
{
        
        switch (n) {

        /* 0 */
        case 0 : return 255.0;          /* 0  0  0 */

        /* 1 */
        case 1 : return r;              /* 1  0  0 */
        case 2 : return g;              /* 0  1  0 */
        case 3 : return b;              /* 0  0  1 */

        /* 2 */
        case 4 : return r * g;          /* 1  1  0 */
        case 5 : return r * b;          /* 1  0  1 */
        case 6 : return g * b;          /* 0  1  1 */
        case 7 : return r * r;          /* 2  0  0 */
        case 8 : return g * g;          /* 0  2  0 */
        case 9 : return b * b;          /* 0  0  2 */

        /* 3 */
        case 10: return r * g * b;      /*  1  1  1 */
        case 11: return r * r * r;      /*  3  0  0 */
        case 12: return g * g * g;      /*  0  3  0 */
        case 13: return b * b * b;      /*  0  0  3 */
        case 14: return r * g * g;      /*  1  2  0 */
        case 15: return r * r * g;      /*  2  1  0 */
        case 16: return g * g * b;      /*  0  2  1          */
        case 17: return b * r * r;      /*  2  0  1 */
        case 18: return b * b * r;      /*  1  0  2 */

        /* 4  */

        case 19: return r * r * g * g;          /* 2  2  0 */
        case 20: return g * g * b * b;          /* 0  2  2 */
        case 21: return r * r * b * b;          /* 2  0  2 */
        case 22: return r * r * g * b;          /* 2  1  1 */
        case 23: return r * g * g * b;          /* 1  2  1    */
        case 24: return r * g * b * b;          /* 1  1  2   */
        case 25: return r * r * r * g;          /* 3  1  0 */
        case 26: return r * r * r * b;          /* 3  0  1 */
        case 27: return r * g * g * g;          /* 1  3  0 */
        case 28: return g * g * g * b;          /* 0  3  1 */
        case 29: return r * b * b * b;          /* 1  0  3 */
        case 30: return g * b * b * b;          /* 0  1  3 */
        case 31: return r * r * r * r;          /* 4  0  0 */
        case 32: return g * g * g * g;          /* 0  4  0 */
        case 33: return b * b * b * b;          /* 0  0  4 */

        /* 5 */

        case 34: return r * r * g * g * b;      /* 2  2  1  */
        case 35: return r * g * g * b * b;      /* 1  2  2 */
        case 36: return r * r * g * b * b;      /* 2  1  2 */
        case 37: return r * r * r * g * g;      /* 3  2  0 */
        case 38: return r * r * r * g * b;      /* 3  1  1 */
        case 39: return r * r * r * b * b;      /* 3  0  2 */
        case 40: return g * g * g * b * b;      /* 0  3  2 */
        case 41: return r * r * g * g * g;      /* 2  3  0 */
        case 42: return r * g * g * g * b;      /* 1  3  1 */
        case 43: return r * r * b * b * b;      /* 2  0  3 */
        case 44: return g * g * b * b * b;      /* 0  2  3 */
        case 45: return r * g * b * b * b;      /* 1  1  3 */
        case 46: return r * r * r * r * g;      /* 4  1  0 */
        case 47: return r * r * r * r * b;      /* 4  0  1 */
        case 48: return r * g * g * g * g;      /* 1  4  0 */
        case 49: return g * g * g * g * b;      /* 0  4  1 */
        case 50: return r * b * b * b * b;      /* 1  0  4 */
        case 51: return g * b * b * b * b;      /* 0  1  4 */
        case 52: return r * r * r * r * r;      /* 5  0  0 */
        case 53: return g * g * g * g * g;      /* 0  5  0 */
        case 54: return b * b * b * b * b;      /* 0  0  5 */

        
        default: return 0;
        }
}



int cmsxRegressionCreateMatrix(LPMEASUREMENT m, SETOFPATCHES Allowed, int nterms,
                                                int ColorSpace,
                                                LPMATN* lpMat, LPMLRSTATISTICS Stat)
{
        LPMATN Input, Output;
        int nCollected = cmsxPCollCountSet(m, Allowed);
        int i, j, n, rc;

        /* We are going always 3 -> 3 for now.... */

        Input  = MATNalloc(nCollected, nterms);
        Output = MATNalloc(nCollected, 3);

        /* Set independent terms */

        for (n = i = 0; i < m -> nPatches; i++)
        {
         if (Allowed[i]) {

                LPPATCH p = m -> Patches + i;

                for (j=0; j < nterms; j++)
                        Input -> Values[n][j] = Term(j, p -> Colorant.RGB[0], p -> Colorant.RGB[1], p->Colorant.RGB[2]);

                switch (ColorSpace) {

                case PT_Lab:

                Output-> Values[n][0] = p -> Lab.L;
                Output-> Values[n][1] = p -> Lab.a;
                Output-> Values[n][2] = p -> Lab.b;
                break;

                case PT_XYZ:
                Output-> Values[n][0] = p -> XYZ.X;
                Output-> Values[n][1] = p -> XYZ.Y;
                Output-> Values[n][2] = p -> XYZ.Z;
                break;


                default:
                        cmsSignalError(LCMS_ERRC_ABORTED, "Invalid colorspace");
                }

        n++;
        }
    }


    /* Apply multiple linear regression */

    if (*lpMat) MATNfree(*lpMat);
    rc = CreateRegressionMatrix(Input, Output, lpMat, Stat);

    /* Free variables */

    MATNfree(Input);
    MATNfree(Output);


#ifdef DEBUG
        if (rc == true)
            MATNprintf("tfm", *lpMat);
#endif

        return rc;
}


/* Convert a RGB triplet to Lab by using regression matrix */

BOOL cmsxRegressionRGB2Lab(double r, double g, double b, LPMATN tfm, LPcmsCIELab Lab)
{
    LPMATN inVec, outVec;
        int i;

    inVec = MATNalloc(1, tfm->Rows);
    if (inVec == NULL)
            return false;

        /* Put terms */
        for (i=0; i < tfm->Rows; i++)
                inVec -> Values[0][i] = Term(i, r, g, b);

        /* Across regression matrix */
        outVec = MATNmult(inVec, tfm);

        /* Store result */
        if (outVec != NULL) {

        Lab->L = outVec->Values[0][0];
        Lab->a = outVec->Values[0][1];
        Lab->b = outVec->Values[0][2];
        MATNfree(outVec);
    }

    MATNfree(inVec);
    return true;
}


/* Convert a RGB triplet to XYX by using regression matrix */

BOOL cmsxRegressionRGB2XYZ(double r, double g, double b, LPMATN tfm, LPcmsCIEXYZ XYZ)
{
    LPMATN inVec, outVec;
        int i;

    inVec = MATNalloc(1, tfm->Rows);
    if (inVec == NULL)
            return false;

        /* Put terms */
        for (i=0; i < tfm->Rows; i++)
                inVec -> Values[0][i] = Term(i, r, g, b);

        /* Across regression matrix */
        outVec = MATNmult(inVec, tfm);

        /* Store result */
        if (outVec != NULL) {

        XYZ->X = outVec->Values[0][0];
        XYZ->Y = outVec->Values[0][1];
        XYZ->Z = outVec->Values[0][2];
        MATNfree(outVec);
    }

    MATNfree(inVec);
    return true;
}


/* Convert a RGB triplet to XYX by using regression matrix */

BOOL cmsxRegressionXYZ2RGB(LPcmsCIEXYZ XYZ, LPMATN tfm, double RGB[3])
{
    LPMATN inVec, outVec;
        int i;

    inVec = MATNalloc(1, tfm->Rows);
    if (inVec == NULL)
            return false;

        /* Put terms */
        for (i=0; i < tfm->Rows; i++)
                inVec -> Values[0][i] = Term(i, XYZ->X, XYZ->Y, XYZ->Z);

        /* Across regression matrix */
        outVec = MATNmult(inVec, tfm);

        /* Store result */
        if (outVec != NULL) {

        RGB[0] = outVec->Values[0][0];
        RGB[1] = outVec->Values[0][1];
        RGB[2] = outVec->Values[0][2];
        MATNfree(outVec);
    }

    MATNfree(inVec);
    return true;
}