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
LPGAMMATABLE cdecl cmsxEstimateGamma(LPSAMPLEDCURVE X, LPSAMPLEDCURVE Y, int nResultingPoints);
42
void             cdecl cmsxCompleteLabOfPatches(LPMEASUREMENT m, SETOFPATCHES Valids, int Medium);
43
 
44
void cdecl cmsxComputeLinearizationTables(LPMEASUREMENT m,
45
                                    int ColorSpace,
46
                                    LPGAMMATABLE Lin[3],
47
                                                                        int nResultingPoints,
48
                                                                        int Medium);
49
 
50
 
51
void cdecl cmsxApplyLinearizationTable(double In[3],
52
                                       LPGAMMATABLE Gamma[3],
53
                                       double Out[3]);
54
 
55
void cdecl cmsxApplyLinearizationGamma(WORD In[3], LPGAMMATABLE Gamma[3], WORD Out[3]);
56
 
57
 
58
 
59
// ------------------------------------------------------------- Implementation
60
 
61
 
62
#define EPSILON 0.00005
63
#define LEVENBERG_MARQUARDT_ITERATE_MAX  150
64
 
65
// In order to track linearization tables, we use following procedure
66
//
67
// We first assume R', G' and B' does exhibit a non-linear behaviour
68
// that can be separated for each channel as Yr(R'), Yg(G'), Yb(B')
69
// This is the shaper step
70
//
71
//  R = Lr(R')
72
//  G = Lg(G')
73
//  B = Lb(B')   (0.0)
74
//
75
// After this step, RGB is converted to XYZ by a matrix multiplication
76
//
77
//  |X|       |R|
78
//  |Y| = [M]|G|
79
//  |Z|       |B|    (1.0)
80
//
81
// In order to extract Lr,Lg,Lb tables, we are interested only on Y part 
82
// 
83
//  Y = (m1 * R + m2 * G + m3 * B)  (1.1)
84
//
85
// The total intensity for maximum RGB = (1, 1, 1) should be 1,
86
//
87
//  1 = m1 * 1 + m2 * 1 + m3 * 1, so
88
//
89
//  m1 + m2 + m3 = 1.0          (1.2)
90
//
91
// We now impose that for neutral (gray) patches, RGB components must be equal
92
//
93
//      R = G = B = Gray 
94
// 
95
// So, substituting in (1.1):
96
// 
97
//  Y = (m1 + m2 + m3) Gray
98
// 
99
// and for (1.2), (m1+m2+m3) = 1, so
100
//
101
//  Y = Gray = Lr(R') = Lg(G') = Lb(B')
102
//
103
// That is, after prelinearization, RGB of gray patches should give 
104
// same values for R, G and B. And this value is Y.
105
//
106
//
107
 
108
 
109
static
110
LPSAMPLEDCURVE NormalizeTo(LPSAMPLEDCURVE X, double N, BOOL lAddEndPoint)
111
{
112
                int i, nItems;
113
                LPSAMPLEDCURVE XNorm;
114
 
115
                nItems = X ->nItems;
116
                if (lAddEndPoint) nItems++;
117
 
118
                XNorm = cmsAllocSampledCurve(nItems);
119
 
120
                for (i=0; i < X ->nItems; i++) {
121
 
122
                        XNorm ->Values[i] = X ->Values[i] / N;                         
123
                }
124
 
125
                if (lAddEndPoint)
126
                        XNorm -> Values[X ->nItems] = 1.0;
127
 
128
                return XNorm;  
129
}
130
 
131
 
132
//
133
// ------------------------------------------------------------------------------
134
//
135
// Our Monitor model. We assume gamma has a general expression of
136
//
137
//  Fn(x) = (Gain * x + offset) ^ gamma | for x >= 0
138
//      Fn(x) = 0                                                       | for x < 0
139
//      
140
//      First partial derivatives are
141
//
142
//  dFn/dGamma  = Fn * ln(Base)
143
//  dFn/dGain   = gamma * x * ((Gain * x + Offset) ^ (gamma -1))
144
//      dFn/dOffset = gamma * ((Gain * x + Offset) ^ (gamma -1))
145
//
146
 
147
static
148
void GammaGainOffsetFn(double x, double *a, double *y, double *dyda, int na)
149
{
150
    double Gamma,Gain,Offset;
151
        double Base;
152
 
153
    Gamma  = a[0];
154
    Gain   = a[1];
155
    Offset = a[2];
156
 
157
    Base = Gain * x + Offset;  
158
 
159
        if (Base < 0) {
160
 
161
                Base = 0.0;
162
                *y = 0.0;    
163
                dyda[0] = 0.0;
164
                dyda[1] = 0.0;
165
                dyda[2] = 0.0;
166
 
167
 
168
        } else {
169
 
170
 
171
                // The function itself
172
                *y = pow(Base, Gamma);
173
 
174
                // dyda[0] is partial derivative across Gamma
175
                dyda[0] = *y * log(Base);
176
 
177
                // dyda[1] is partial derivative across gain
178
                dyda[1] = (x * Gamma) * pow(Base, Gamma-1.0);
179
 
180
                // dyda[2] is partial derivative across offset
181
                dyda[2] =     Gamma * pow(Base, Gamma-1.0);
182
        }
183
}
184
 
185
 
186
// Fit curve to our gamma-gain-offset model.
187
 
188
static
189
BOOL OneTry(LPSAMPLEDCURVE XNorm, LPSAMPLEDCURVE YNorm, double a[])
190
{
191
        LCMSHANDLE h;
192
        double ChiSq, OldChiSq;        
193
        int i;
194
        BOOL Status = TRUE;
195
 
196
                // initial guesses
197
 
198
            a[0] = 3.0;                 // gamma
199
            a[1] = 4.0;                 // gain
200
            a[2] = 6.0;                 // offset
201
                a[3] = 0.0;                     // Thereshold
202
                a[4] = 0.0;                     // Black
203
 
204
 
205
                // Significance = 0.02 gives good results
206
 
207
                h = cmsxLevenbergMarquardtInit(XNorm, YNorm,  0.02, a, 3, GammaGainOffsetFn);                                                  
208
                if (h == NULL) return FALSE;
209
 
210
 
211
                OldChiSq = cmsxLevenbergMarquardtChiSq(h);
212
 
213
                for(i = 0; i < LEVENBERG_MARQUARDT_ITERATE_MAX; i++) {
214
 
215
                        if (!cmsxLevenbergMarquardtIterate(h)) {
216
                                Status = FALSE;
217
                                break;
218
                        }
219
 
220
                        ChiSq = cmsxLevenbergMarquardtChiSq(h);
221
 
222
                        if(OldChiSq != ChiSq && (OldChiSq - ChiSq) < EPSILON)
223
                                break;
224
 
225
                        OldChiSq = ChiSq;
226
                }
227
 
228
                cmsxLevenbergMarquardtFree(h);
229
 
230
                return Status;
231
}
232
 
233
// Tries to fit gamma as per IEC 61966-2.1 using Levenberg-Marquardt method
234
// 
235
// Y = (aX + b)^Gamma | X >= d
236
// Y = cX             | X < d
237
 
238
LPGAMMATABLE cmsxEstimateGamma(LPSAMPLEDCURVE X, LPSAMPLEDCURVE Y, int nResultingPoints)
239
{
240
        double a[5];
241
        LPSAMPLEDCURVE XNorm, YNorm;   
242
        double e, Max;
243
 
244
 
245
                // Coarse approximation, to find maximum. 
246
            // We have only a portion of curve. It is likely
247
            // maximum will not fall on exactly 100.
248
 
249
                if (!OneTry(X, Y, a))
250
                        return FALSE;
251
 
252
                // Got parameters. Compute maximum.
253
                e = a[1]* 255.0 + a[2];
254
                if (e < 0) return FALSE;
255
                Max = pow(e, a[0]);
256
 
257
 
258
                // Normalize values to maximum
259
                XNorm = NormalizeTo(X, 255.0, FALSE);
260
                YNorm = NormalizeTo(Y, Max, FALSE);
261
 
262
                // Do the final fitting 
263
                if (!OneTry(XNorm, YNorm, a))
264
                                return FALSE;
265
 
266
                // Type 3 = IEC 61966-2.1 (sRGB)
267
        // Y = (aX + b)^Gamma | X >= d
268
        // Y = cX             | X < d
269
                return  cmsBuildParametricGamma(nResultingPoints, 3, a);
270
}
271
 
272
 
273
 
274
 
275
 
276
// A dumb bubble sort
277
 
278
static
279
void Bubble(LPSAMPLEDCURVE C, LPSAMPLEDCURVE L)
280
{
281
#define SWAP(a, b)      { tmp = (a); (a) = (b); (b) = tmp; }
282
 
283
        BOOL lSwapped;
284
        int i, nItems;
285
        double tmp;
286
 
287
                nItems = C -> nItems;
288
        do {
289
                lSwapped = FALSE;
290
 
291
                for (i= 0; i <  nItems - 1; i++) {
292
 
293
                        if (C->Values[i] > C->Values[i+1]) {
294
 
295
                                SWAP(C->Values[i], C->Values[i+1]);
296
                                                                SWAP(L->Values[i], L->Values[i+1]);                            
297
                                lSwapped = TRUE;
298
                        }
299
                }
300
 
301
        } while (lSwapped);
302
 
303
#undef SWAP
304
}
305
 
306
 
307
 
308
// Check for monotonicity. Force it if is not the case.
309
 
310
static
311
void CheckForMonotonicSampledCurve(LPSAMPLEDCURVE t)
312
{
313
    int n = t ->nItems;
314
    int i;
315
        double last;
316
 
317
    last = t ->Values[n-1];
318
    for (i = n-2; i >= 0; --i) {
319
 
320
        if (t ->Values[i] > last)
321
 
322
                t ->Values[i] = last;
323
        else
324
                last = t ->Values[i];
325
 
326
    }
327
 
328
}
329
 
330
// The main gamma inferer. Tries first by gamma-gain-offset, 
331
// if not proper reverts to curve guessing.
332
 
333
static
334
LPGAMMATABLE BuildGammaTable(LPSAMPLEDCURVE C, LPSAMPLEDCURVE L, int nResultingPoints)
335
{
336
        LPSAMPLEDCURVE Cw, Lw, Cn, Ln;
337
        LPSAMPLEDCURVE out;
338
        LPGAMMATABLE Result;
339
        double Lmax, Lend, Cmax;
340
 
341
        // Try to see if it can be fitted       
342
        Result = cmsxEstimateGamma(C, L, nResultingPoints);  
343
        if (Result)
344
                return Result;
345
 
346
 
347
        // No... build curve from scratch. Since we have not
348
        // endpoints, a coarse linear extrapolation should be
349
        // applied in order to get the expected maximum.
350
 
351
        Cw = cmsDupSampledCurve(C);
352
        Lw = cmsDupSampledCurve(L);
353
 
354
        Bubble(Cw, Lw);
355
 
356
    // Get endpoint
357
        Lmax = Lw->Values[Lw ->nItems - 1];
358
        Cmax = Cw->Values[Cw ->nItems - 1];
359
 
360
        // Linearly extrapolate
361
        Lend = (255 * Lmax) / Cmax;
362
 
363
        Ln = NormalizeTo(Lw, Lend, TRUE);
364
        Cn = NormalizeTo(Cw, 255.0, TRUE);
365
 
366
        cmsFreeSampledCurve(Cw);
367
        cmsFreeSampledCurve(Lw);
368
 
369
        // Add endpoint
370
        out = cmsJoinSampledCurves(Cn, Ln,  nResultingPoints);
371
 
372
        cmsFreeSampledCurve(Cn);
373
        cmsFreeSampledCurve(Ln);
374
 
375
        CheckForMonotonicSampledCurve(out);
376
 
377
        cmsSmoothSampledCurve(out, nResultingPoints*4.);
378
        cmsClampSampledCurve(out, 0, 1.0);                     
379
 
380
        Result = cmsConvertSampledCurveToGamma(out, 1.0);    
381
 
382
        cmsFreeSampledCurve(out);
383
        return Result;
384
}
385
 
386
 
387
 
388
 
389
void cmsxCompleteLabOfPatches(LPMEASUREMENT m, SETOFPATCHES Valids, int Medium)
390
{
391
        LPPATCH White; 
392
        cmsCIEXYZ WhiteXYZ;
393
        int i;
394
 
395
        if (Medium == MEDIUM_REFLECTIVE_D50)
396
        {
397
                WhiteXYZ.X = D50X * 100.;
398
                WhiteXYZ.Y = D50Y * 100.;
399
                WhiteXYZ.Z = D50Z * 100.;
400
        }
401
        else {
402
 
403
                White = cmsxPCollFindWhite(m, Valids, NULL);
404
                if (!White) return;
405
 
406
                WhiteXYZ = White ->XYZ;
407
        }
408
 
409
        // For all patches with XYZ and without Lab, add Lab values.
410
        // Transmissive profiles does need to locate its own white
411
        // point for device gray. Reflective does use D50
412
 
413
         for (i=0; i < m -> nPatches; i++) {
414
 
415
       if (Valids[i]) {
416
 
417
        LPPATCH p = m -> Patches + i;
418
 
419
                if ((p ->dwFlags & PATCH_HAS_XYZ) &&
420
                        (!(p ->dwFlags & PATCH_HAS_Lab) || (Medium == MEDIUM_TRANSMISSIVE))) {
421
 
422
                        cmsXYZ2Lab(&WhiteXYZ, &p->Lab, &p->XYZ);
423
                        p -> dwFlags |= PATCH_HAS_Lab;
424
                }
425
           }
426
         }
427
}
428
 
429
 
430
// Compute linearization tables, trying to fit in a pure 
431
// exponential gamma. If gamma cannot be accurately infered,
432
// then does build a smooth, monotonic curve that does the job.
433
 
434
void cmsxComputeLinearizationTables(LPMEASUREMENT m,
435
                                    int ColorSpace,
436
                                    LPGAMMATABLE Lin[3],
437
                                                                        int nResultingPoints,
438
                                                                        int Medium)                                    
439
 
440
{
441
    LPSAMPLEDCURVE R, G, B, L;
442
        LPGAMMATABLE gr, gg, gb;
443
    SETOFPATCHES Neutrals;
444
    int nGrays;
445
    int i;
446
 
447
        // We need Lab for grays.
448
        cmsxCompleteLabOfPatches(m, m->Allowed, Medium);
449
 
450
    // Add neutrals, normalize to max
451
    Neutrals = cmsxPCollBuildSet(m, FALSE);
452
    cmsxPCollPatchesNearNeutral(m, m ->Allowed, 15, Neutrals); 
453
 
454
    nGrays = cmsxPCollCountSet(m, Neutrals);
455
 
456
    R = cmsAllocSampledCurve(nGrays);  
457
    G = cmsAllocSampledCurve(nGrays);
458
    B = cmsAllocSampledCurve(nGrays);
459
    L = cmsAllocSampledCurve(nGrays);
460
 
461
        nGrays = 0;
462
 
463
    // Collect patches   
464
    for (i=0; i < m -> nPatches; i++) {
465
 
466
                if (Neutrals[i]) {
467
 
468
                        LPPATCH gr = m -> Patches + i;
469
 
470
 
471
                        R -> Values[nGrays] = gr -> Colorant.RGB[0];
472
                        G -> Values[nGrays] = gr -> Colorant.RGB[1];
473
                        B -> Values[nGrays] = gr -> Colorant.RGB[2];
474
                                                L -> Values[nGrays] = gr -> XYZ.Y;
475
 
476
                        nGrays++;
477
                }
478
 
479
    }
480
 
481
 
482
        gr = BuildGammaTable(R, L, nResultingPoints);
483
        gg = BuildGammaTable(G, L, nResultingPoints);
484
        gb = BuildGammaTable(B, L, nResultingPoints);
485
 
486
        cmsFreeSampledCurve(R);
487
        cmsFreeSampledCurve(G);
488
        cmsFreeSampledCurve(B);
489
        cmsFreeSampledCurve(L);
490
 
491
        if (ColorSpace == PT_Lab) {
492
 
493
                LPGAMMATABLE Gamma3 = cmsBuildGamma(nResultingPoints, 3.0);
494
 
495
                Lin[0] = cmsJoinGammaEx(gr, Gamma3, nResultingPoints);
496
                Lin[1] = cmsJoinGammaEx(gg, Gamma3, nResultingPoints);
497
                Lin[2] = cmsJoinGammaEx(gb, Gamma3, nResultingPoints);
498
 
499
                cmsFreeGamma(gr); cmsFreeGamma(gg); cmsFreeGamma(gb);
500
                cmsFreeGamma(Gamma3);
501
        }
502
        else {
503
 
504
 
505
                LPGAMMATABLE Gamma1 = cmsBuildGamma(nResultingPoints, 1.0);
506
 
507
                Lin[0] = cmsJoinGammaEx(gr, Gamma1, nResultingPoints);
508
                Lin[1] = cmsJoinGammaEx(gg, Gamma1, nResultingPoints);
509
                Lin[2] = cmsJoinGammaEx(gb, Gamma1, nResultingPoints);
510
 
511
                cmsFreeGamma(gr); cmsFreeGamma(gg); cmsFreeGamma(gb);
512
                cmsFreeGamma(Gamma1);
513
 
514
        }
515
 
516
}
517
 
518
 
519
 
520
// Apply linearization. WORD encoded version
521
 
522
void cmsxApplyLinearizationGamma(WORD In[3], LPGAMMATABLE Gamma[3], WORD Out[3])
523
{
524
        L16PARAMS Lut16;
525
 
526
        cmsCalcL16Params(Gamma[0] -> nEntries, &Lut16);
527
 
528
        Out[0] = cmsLinearInterpLUT16(In[0], Gamma[0] -> GammaTable, &Lut16);
529
        Out[1] = cmsLinearInterpLUT16(In[1], Gamma[1] -> GammaTable, &Lut16);
530
        Out[2] = cmsLinearInterpLUT16(In[2], Gamma[2] -> GammaTable, &Lut16);
531
 
532
 
533
}
534
 
535
 
536
 
537
// Apply linearization. double version
538
 
539
void cmsxApplyLinearizationTable(double In[3], LPGAMMATABLE Gamma[3], double Out[3])
540
{
541
        WORD rw, gw, bw;
542
        double rd, gd, bd;
543
        L16PARAMS Lut16;
544
 
545
 
546
        cmsCalcL16Params(Gamma[0] -> nEntries, &Lut16);
547
 
548
                rw = (WORD) floor(_cmsxSaturate255To65535(In[0]) + .5);
549
                gw = (WORD) floor(_cmsxSaturate255To65535(In[1]) + .5);
550
                bw = (WORD) floor(_cmsxSaturate255To65535(In[2]) + .5);
551
 
552
        rd = cmsLinearInterpLUT16(rw , Gamma[0] -> GammaTable, &Lut16);
553
        gd = cmsLinearInterpLUT16(gw,  Gamma[1] -> GammaTable, &Lut16);
554
        bd = cmsLinearInterpLUT16(bw,  Gamma[2] -> GammaTable, &Lut16);
555
 
556
        Out[0] = _cmsxSaturate65535To255(rd);            // back to 0..255
557
        Out[1] = _cmsxSaturate65535To255(gd);
558
        Out[2] = _cmsxSaturate65535To255(bd);
559
}
560