/* */ /* Little cms - profiler construction set */ /* Copyright (C) 1998-2001 Marti Maria */ /* */ /* 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; }