Subversion Repositories Scribus

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2454 craig 1
//
2
//  Little cms - profiler construction set
3
//  Copyright (C) 1998-2001 Marti Maria
4
//
5
// THIS SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND,
6
// EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
7
// WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
8
//
9
// IN NO EVENT SHALL MARTI MARIA BE LIABLE FOR ANY SPECIAL, INCIDENTAL,
10
// INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND,
11
// OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
12
// WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF
13
// LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
14
// OF THIS SOFTWARE.
15
//
16
// This file is free software; you can redistribute it and/or modify it
17
// under the terms of the GNU General Public License as published by
18
// the Free Software Foundation; either version 2 of the License, or
19
// (at your option) any later version.
20
//
21
// This program is distributed in the hope that it will be useful, but
22
// WITHOUT ANY WARRANTY; without even the implied warranty of
23
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
24
// General Public License for more details.
25
//
26
// You should have received a copy of the GNU General Public License
27
// along with this program; if not, write to the Free Software
28
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
29
//
30
// As a special exception to the GNU General Public License, if you
31
// distribute this file as part of a program that contains a
32
// configuration script generated by Autoconf, you may include it under
33
// the same distribution terms that you use for the rest of that program.
34
//
35
// Version 1.09a
36
 
37
 
38
#include "lcmsprf.h"
39
 
40
 
41
// There are three kinds of lies:
42
//
43
//   * lies
44
//   * damn lies
45
//   * statistics
46
//
47
//                    -Some Wag
48
//
49
//
50
// This module handles multiple linear regression stuff
51
 
52
 
53
/*
54
// A measurement of error
55
 
56
typedef struct {
57
 
58
        double SSE;             // The error sum of squares
59
        double MSE;             // The error mean sum of squares
60
        double SSR;             // The regression sum of squares
61
        double MSR;             // The regression mean sum of squares
62
        double SSTO;            // Total sum of squares
63
        double F;               // The Fisher-F value (MSR / MSE)
64
        double R2;              // Proportion of variability explained by the regression
65
                                // (root is Pearson correlation coefficient)
66
 
67
        double R2adj;           // The adjusted coefficient of multiple determination.
68
                                // R2-adjusted or R2adj. This is calculated as
69
                                // R2adj = 1 - (1-R2)(N-n-1)/(N-1)
70
                                // and used as multiple correlation coefficient
71
                                // (really, it should be square root)
72
 
73
    } MLRSTATISTICS, FAR* LPMLRSTATISTICS;
74
 
75
*/
76
 
77
 
78
int  cdecl cmsxRegressionCreateMatrix(LPMEASUREMENT m, SETOFPATCHES Allowed, int nterms,
79
                                       int ColorSpace,
80
                                       LPMATN* lpMat, LPMLRSTATISTICS Stat);
81
 
82
BOOL cdecl cmsxRegressionRGB2Lab(double r, double g, double b,
83
                                       LPMATN tfm, LPcmsCIELab Lab);
84
 
85
BOOL cdecl cmsxRegressionRGB2XYZ(double r, double g, double b,
86
                                       LPMATN tfm, LPcmsCIEXYZ XYZ);
87
 
88
 
89
// -------------------------------------------------------------- Implementation
90
 
91
// #define DEBUG 1
92
 
93
 
94
// Multiple linear regression. Also keep track of error.
95
// Returns FALSE if something goes wrong, or TRUE if all Ok.
96
 
97
static
98
BOOL MultipleLinearRegression(const LPMATN xi,                // Dependent variable
99
                              const LPMATN y,                 // Independent variable
100
                              int nvar,                       // Number of samples
101
                              int npar,                       // Number of parameters (terms)
102
                              double* coeff,                  // Returned coefficients
103
                              LPMATN vcb,                     // Variance-covariance array
104
                              double *tvl,                    // T-Values
105
                              LPMLRSTATISTICS ans)            // The returned statistics
106
{
107
    LPMATN bt, xt, a, xy, yt, b;
108
    double  sum;
109
    LPMATN temp1, temp2;
110
    int i;
111
 
112
 
113
    // |xt| = |xi| T
114
    xt = MATNtranspose(xi);
115
    if (xt == NULL) return FALSE;
116
 
117
 
118
    // |a| = |xt|* |xi|
119
    a = MATNmult(xt, xi);
120
    if (a == NULL) return FALSE;
121
 
122
 
123
    // |xy| = |xy| * |y|
124
    xy = MATNmult (xt, y);
125
    if (xy == NULL) return FALSE;
126
 
127
 
128
    // solve system |a|*|xy| = 0
129
    if (!MATNsolve(a, xy)) return FALSE;
130
 
131
    // b will hold coefficients
132
    b = MATNalloc (xy->Rows, 1);
133
    if (b == NULL) return FALSE;
134
 
135
    for (i = 0; i < npar; i++)
136
        b->Values[i][0] = xy->Values[i][0];
137
 
138
    // Store a copy for later user
139
    for (i = 0; i < npar; i++)
140
        coeff[i] = b->Values[i][0];
141
 
142
    // Error analysis.
143
 
144
    // SSE and MSE.
145
    temp1 = MATNalloc (1,1);
146
    if ((temp1->Values[0][0] = MATNcross(y)) == 0) return FALSE;
147
 
148
    // |bt| = |b| T
149
    bt = MATNtranspose (b);
150
    if (bt == FALSE) return FALSE;
151
 
152
    //  |yt| = |bt| * |xt|
153
    yt = MATNmult (bt, xt);
154
    if (yt == NULL) return FALSE;
155
 
156
 
157
    // |temp2| = |yt|* |y|
158
    temp2 = MATNmult (yt, y);
159
    if (temp2 == NULL) return FALSE;
160
 
161
    // SSE, MSE
162
    ans->SSE = temp1 -> Values[0][0] - temp2 -> Values[0][0];
163
    ans->MSE = ans->SSE / (double) (nvar - npar);
164
 
165
    // SSTO
166
    sum = 0;
167
    for (i=0; i < nvar; i++)
168
                sum += y->Values[i][0];
169
 
170
    sum *= sum / (double) nvar;
171
    ans->SSTO = temp1->Values[0][0] - sum;
172
 
173
    // SSR, MSR, and Fisher-F
174
    ans->SSR = temp2->Values[0][0] - sum;
175
    ans->MSR = ans->SSR / (double) (npar - 1);
176
    ans->F   = ans->MSR / ans->MSE;
177
 
178
    // Correlation coefficients.
179
    ans->R2 = ans->SSR/ans->SSTO;
180
    ans->R2adj = 1.0 - (ans->SSE/ans->SSTO)*((nvar-1.)/(nvar-npar));
181
 
182
    // Variance-covariance matrix
183
    //
184
    //     In RGB->Lab, for example:
185
    //
186
    //          Var(R)    Cov(R,G)  Cov(R,B)
187
    // |vcb| =  Cov(R,G)  Var(G)    Cov(G,B)
188
    //          Cov(R,B)  Cov(G,B)  Var(B)
189
    //
190
 
191
    MATNscalar(a, ans->MSE, vcb);
192
 
193
    // Determine the T-values
194
 
195
    for (i=0; i < npar; i++) {
196
 
197
        temp1->Values[0][0] = fabs(vcb->Values[i][0]);
198
        if ( temp1->Values[0][0] == 0)
199
            tvl[i] = 0;             // This should never happen
200
        else
201
            tvl[i] = b->Values[i][0] / sqrt(temp1->Values[0][0]);
202
    }
203
 
204
 
205
    // Ok, done
206
 
207
    MATNfree(a); MATNfree(xy); MATNfree(yt); MATNfree(b);
208
    MATNfree(temp1); MATNfree(temp2); MATNfree(bt); MATNfree(xt);
209
 
210
 
211
    return TRUE;
212
}
213
 
214
 
215
 
216
// Does create (so, it allocates) the regression matrix,
217
// keeping track of error as well.
218
 
219
static
220
BOOL CreateRegressionMatrix(const LPMATN Input, const LPMATN Output,
221
                            LPMATN* ptrMatrix, LPMLRSTATISTICS maxErrorMeas)
222
{
223
    double* coef;
224
    double* tval;
225
    LPMATN ivar, dvar, vcov;
226
    MLRSTATISTICS ErrorMeas, PeakErrorMeas;
227
    int i, j, nIn, nOut, NumOfPatches;
228
 
229
        nIn          = Input  -> Cols;
230
        nOut         = Output -> Cols;
231
        NumOfPatches = Input  -> Rows;
232
 
233
        // Checkpoint
234
        if (Output -> Rows != NumOfPatches) {
235
 
236
                cmsSignalError(LCMS_ERRC_ABORTED, "(internal) Regression matrix mismatch");
237
                return FALSE;
238
        }
239
 
240
    coef = (double*) malloc(nIn * sizeof(double));
241
    if (coef == NULL) return FALSE;
242
 
243
    tval = (double*) malloc(nIn * sizeof(double));
244
    if (tval == NULL) {
245
        free(coef);
246
        return FALSE;
247
    }
248
 
249
    ivar = MATNalloc(NumOfPatches, nIn);
250
    dvar = MATNalloc(NumOfPatches, 1);
251
 
252
    // Copy In to ivar,
253
    for (i = 0; i < NumOfPatches; i++) {
254
 
255
        for (j = 0; j < nIn; j++)
256
            ivar->Values[i][j] = Input->Values[i][j];
257
    }
258
 
259
    // This is the (symmetric) Covariance matrix
260
    vcov = MATNalloc(nIn, nIn);
261
 
262
    // This is the regression matrix
263
    *ptrMatrix = MATNalloc(nIn, nOut);
264
 
265
    PeakErrorMeas.R2adj = 0;
266
    for (j = 0; j < nOut; ++j)
267
    {
268
        for (i = 0; i < NumOfPatches; ++i)
269
            dvar->Values[i][0] = Output->Values[i][j];
270
 
271
        if (MultipleLinearRegression(ivar, dvar, NumOfPatches, nIn, coef, vcov, tval, &ErrorMeas)) {
272
 
273
            // Ok so far... store values
274
            for (i = 0; i < nIn; i++)
275
                (*ptrMatrix)->Values[i][j] = coef[i];
276
        }
277
        else {
278
                    // Boo... got error. Discard whole point.
279
                    MATNfree(ivar); MATNfree(dvar); MATNfree(vcov);
280
                    if (coef) free(coef);
281
                    if (tval) free(tval);
282
                    MATNfree(*ptrMatrix); *ptrMatrix = NULL;
283
                    return FALSE;
284
        }
285
 
286
        // Did this colorant got higer error? If so, this is 
287
        // the peak of all pixel
288
 
289
        if(fabs(ErrorMeas.R2adj) > fabs(PeakErrorMeas.R2adj))
290
            PeakErrorMeas = ErrorMeas;
291
    }
292
 
293
    // This is the peak error on all components
294
    *maxErrorMeas = PeakErrorMeas;
295
 
296
 
297
#ifdef DEBUG
298
        MATNprintf("Variance-Covariance", vcov);
299
        printf("R2adj: %g, F: %g\n", PeakErrorMeas.R2adj, PeakErrorMeas.F);
300
#endif
301
 
302
    // Free stuff.
303
    MATNfree(ivar); MATNfree(dvar); MATNfree(vcov);
304
    if (coef) free(coef);
305
    if (tval) free(tval);
306
 
307
    return TRUE;
308
}
309
 
310
 
311
// Does compute the term of regression based on inputs.
312
 
313
static
314
double Term(int n, double r, double g, double b)
315
{
316
 
317
        switch (n) {
318
 
319
        // 0
320
        case 0 : return 255.0;          // 0  0  0
321
 
322
        // 1
323
        case 1 : return r;              // 1  0  0
324
        case 2 : return g;              // 0  1  0
325
        case 3 : return b;              // 0  0  1
326
 
327
        // 2
328
        case 4 : return r * g;          // 1  1  0
329
        case 5 : return r * b;          // 1  0  1
330
        case 6 : return g * b;          // 0  1  1
331
        case 7 : return r * r;          // 2  0  0
332
        case 8 : return g * g;          // 0  2  0
333
        case 9 : return b * b;          // 0  0  2
334
 
335
        // 3
336
        case 10: return r * g * b;      //  1  1  1
337
        case 11: return r * r * r;      //  3  0  0
338
        case 12: return g * g * g;      //  0  3  0
339
        case 13: return b * b * b;      //  0  0  3
340
        case 14: return r * g * g;      //  1  2  0
341
        case 15: return r * r * g;      //  2  1  0
342
        case 16: return g * g * b;      //  0  2  1         
343
        case 17: return b * r * r;      //  2  0  1
344
        case 18: return b * b * r;      //  1  0  2
345
 
346
        // 4 
347
 
348
        case 19: return r * r * g * g;          // 2  2  0
349
        case 20: return g * g * b * b;          // 0  2  2
350
        case 21: return r * r * b * b;          // 2  0  2
351
        case 22: return r * r * g * b;          // 2  1  1
352
        case 23: return r * g * g * b;          // 1  2  1   
353
        case 24: return r * g * b * b;          // 1  1  2  
354
        case 25: return r * r * r * g;          // 3  1  0
355
        case 26: return r * r * r * b;          // 3  0  1
356
        case 27: return r * g * g * g;          // 1  3  0
357
        case 28: return g * g * g * b;          // 0  3  1
358
        case 29: return r * b * b * b;          // 1  0  3
359
        case 30: return g * b * b * b;          // 0  1  3
360
        case 31: return r * r * r * r;          // 4  0  0
361
        case 32: return g * g * g * g;          // 0  4  0
362
        case 33: return b * b * b * b;          // 0  0  4
363
 
364
        // 5
365
 
366
        case 34: return r * r * g * g * b;      // 2  2  1 
367
        case 35: return r * g * g * b * b;      // 1  2  2
368
        case 36: return r * r * g * b * b;      // 2  1  2
369
        case 37: return r * r * r * g * g;      // 3  2  0
370
        case 38: return r * r * r * g * b;      // 3  1  1
371
        case 39: return r * r * r * b * b;      // 3  0  2
372
        case 40: return g * g * g * b * b;      // 0  3  2
373
        case 41: return r * r * g * g * g;      // 2  3  0
374
        case 42: return r * g * g * g * b;      // 1  3  1
375
        case 43: return r * r * b * b * b;      // 2  0  3
376
        case 44: return g * g * b * b * b;      // 0  2  3
377
        case 45: return r * g * b * b * b;      // 1  1  3
378
        case 46: return r * r * r * r * g;      // 4  1  0
379
        case 47: return r * r * r * r * b;      // 4  0  1
380
        case 48: return r * g * g * g * g;      // 1  4  0
381
        case 49: return g * g * g * g * b;      // 0  4  1
382
        case 50: return r * b * b * b * b;      // 1  0  4
383
        case 51: return g * b * b * b * b;      // 0  1  4
384
        case 52: return r * r * r * r * r;      // 5  0  0
385
        case 53: return g * g * g * g * g;      // 0  5  0
386
        case 54: return b * b * b * b * b;      // 0  0  5
387
 
388
 
389
        default: return 0;
390
        }
391
}
392
 
393
 
394
 
395
int cmsxRegressionCreateMatrix(LPMEASUREMENT m, SETOFPATCHES Allowed, int nterms,
396
                                                int ColorSpace,
397
                                                LPMATN* lpMat, LPMLRSTATISTICS Stat)
398
{
399
        LPMATN Input, Output;
400
        int nCollected = cmsxPCollCountSet(m, Allowed);
401
        int i, j, n, rc;
402
 
403
        // We are going always 3 -> 3 for now....
404
 
405
        Input  = MATNalloc(nCollected, nterms);
406
        Output = MATNalloc(nCollected, 3);
407
 
408
        // Set independent terms
409
 
410
        for (n = i = 0; i < m -> nPatches; i++)
411
        {
412
         if (Allowed[i]) {
413
 
414
                LPPATCH p = m -> Patches + i;
415
 
416
                for (j=0; j < nterms; j++)
417
                        Input -> Values[n][j] = Term(j, p -> Colorant.RGB[0], p -> Colorant.RGB[1], p->Colorant.RGB[2]);
418
 
419
                switch (ColorSpace) {
420
 
421
                case PT_Lab:
422
 
423
                Output-> Values[n][0] = p -> Lab.L;
424
                Output-> Values[n][1] = p -> Lab.a;
425
                Output-> Values[n][2] = p -> Lab.b;
426
                break;
427
 
428
                case PT_XYZ:
429
                Output-> Values[n][0] = p -> XYZ.X;
430
                Output-> Values[n][1] = p -> XYZ.Y;
431
                Output-> Values[n][2] = p -> XYZ.Z;
432
                break;
433
 
434
 
435
                default:
436
                        cmsSignalError(LCMS_ERRC_ABORTED, "Invalid colorspace");
437
                }
438
 
439
        n++;
440
        }
441
    }
442
 
443
 
444
    // Apply multiple linear regression
445
 
446
    if (*lpMat) MATNfree(*lpMat);
447
    rc = CreateRegressionMatrix(Input, Output, lpMat, Stat);
448
 
449
    // Free variables
450
 
451
    MATNfree(Input);
452
    MATNfree(Output);
453
 
454
 
455
#ifdef DEBUG
456
        if (rc == TRUE)
457
            MATNprintf("tfm", *lpMat);
458
#endif
459
 
460
        return rc;
461
}
462
 
463
 
464
// Convert a RGB triplet to Lab by using regression matrix
465
 
466
BOOL cmsxRegressionRGB2Lab(double r, double g, double b, LPMATN tfm, LPcmsCIELab Lab)
467
{
468
    LPMATN inVec, outVec;
469
        int i;
470
 
471
    inVec = MATNalloc(1, tfm->Rows);
472
    if (inVec == NULL)
473
            return FALSE;
474
 
475
        // Put terms
476
        for (i=0; i < tfm->Rows; i++)
477
                inVec -> Values[0][i] = Term(i, r, g, b);
478
 
479
        // Across regression matrix
480
        outVec = MATNmult(inVec, tfm);
481
 
482
        // Store result
483
        if (outVec != NULL) {
484
 
485
        Lab->L = outVec->Values[0][0];
486
        Lab->a = outVec->Values[0][1];
487
        Lab->b = outVec->Values[0][2];
488
        MATNfree(outVec);
489
    }
490
 
491
    MATNfree(inVec);
492
    return TRUE;
493
}
494
 
495
 
496
// Convert a RGB triplet to XYX by using regression matrix
497
 
498
BOOL cmsxRegressionRGB2XYZ(double r, double g, double b, LPMATN tfm, LPcmsCIEXYZ XYZ)
499
{
500
    LPMATN inVec, outVec;
501
        int i;
502
 
503
    inVec = MATNalloc(1, tfm->Rows);
504
    if (inVec == NULL)
505
            return FALSE;
506
 
507
        // Put terms
508
        for (i=0; i < tfm->Rows; i++)
509
                inVec -> Values[0][i] = Term(i, r, g, b);
510
 
511
        // Across regression matrix
512
        outVec = MATNmult(inVec, tfm);
513
 
514
        // Store result
515
        if (outVec != NULL) {
516
 
517
        XYZ->X = outVec->Values[0][0];
518
        XYZ->Y = outVec->Values[0][1];
519
        XYZ->Z = outVec->Values[0][2];
520
        MATNfree(outVec);
521
    }
522
 
523
    MATNfree(inVec);
524
    return TRUE;
525
}
526
 
527
 
528
// Convert a RGB triplet to XYX by using regression matrix
529
 
530
BOOL cmsxRegressionXYZ2RGB(LPcmsCIEXYZ XYZ, LPMATN tfm, double RGB[3])
531
{
532
    LPMATN inVec, outVec;
533
        int i;
534
 
535
    inVec = MATNalloc(1, tfm->Rows);
536
    if (inVec == NULL)
537
            return FALSE;
538
 
539
        // Put terms
540
        for (i=0; i < tfm->Rows; i++)
541
                inVec -> Values[0][i] = Term(i, XYZ->X, XYZ->Y, XYZ->Z);
542
 
543
        // Across regression matrix
544
        outVec = MATNmult(inVec, tfm);
545
 
546
        // Store result
547
        if (outVec != NULL) {
548
 
549
        RGB[0] = outVec->Values[0][0];
550
        RGB[1] = outVec->Values[0][1];
551
        RGB[2] = outVec->Values[0][2];
552
        MATNfree(outVec);
553
    }
554
 
555
    MATNfree(inVec);
556
    return TRUE;
557
}
558