diff options
Diffstat (limited to 'src/libs/lprof/cmsreg.cpp')
-rw-r--r-- | src/libs/lprof/cmsreg.cpp | 558 |
1 files changed, 558 insertions, 0 deletions
diff --git a/src/libs/lprof/cmsreg.cpp b/src/libs/lprof/cmsreg.cpp new file mode 100644 index 00000000..4e66a9ef --- /dev/null +++ b/src/libs/lprof/cmsreg.cpp @@ -0,0 +1,558 @@ +/* */ +/* 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; +} + |