2
Copyright (C) 1999-2002 Adam D. Moss (the "Author"). All Rights Reserved.
4
Permission is hereby granted, free of charge, to any person obtaining a copy
5
of this software and associated documentation files (the "Software"), to deal
6
in the Software without restriction, including without limitation the rights
7
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8
copies of the Software, and to permit persons to whom the Software is fur-
9
nished to do so, subject to the following conditions:
11
The above copyright notice and this permission notice shall be included in
12
all copies or substantial portions of the Software.
14
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FIT-
16
NESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17
AUTHOR BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
18
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CON-
19
NECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21
Except as contained in this notice, the name of the Author of the
22
Software shall not be used in advertising or otherwise to promote the sale,
23
use or other dealings in this Software without prior written authorization
28
cpercep.c: The CPercep Functions v0.9: 2002-02-10
29
Adam D. Moss: adam@gimp.org <http://www.foxbox.org/adam/code/cpercep/>
31
This code module concerns itself with conversion from a hard-coded
32
RGB colour space (sRGB by default) to CIE L*a*b* and back again with
33
(primarily) precision and (secondarily) speed, oriented largely
34
towards the purposes of quantifying the PERCEPTUAL difference between
35
two arbitrary RGB colours with a minimum of fuss.
37
Motivation One: The author is disheartened at the amount of graphics
38
processing software around which uses weighted or non-weighted
39
Euclidean distance between co-ordinates within a (poorly-defined) RGB
40
space as the basis of what should really be an estimate of perceptual
41
difference to the human eye. Certainly it's fast to do it that way,
42
but please think carefully about whether your particular application
43
should be tolerating sloppy results for the sake of real-time response.
45
Motivation Two: Lack of tested, re-usable and free code available
46
for this purpose. The difficulty in finding something similar to
47
CPercep with a free license motivated this project; I hope that this
48
code also serves to illustrate how to perform the
49
R'G'B'->XYZ->L*a*b*->XYZ->R'G'B' transformation correctly since I
50
was distressed to note how many of the equations and code snippets
51
on the net were omitting the reverse transform and/or were using
52
incorrectly-derived or just plain wrong constants.
54
TODO: document functions, rename erroneously-named arguments
61
/* cbrt() is a GNU extension */
62
#define cbrt(x) (pow(x, 1.0/3.0))
65
#ifdef GIMP_COMPILATION
66
#include <glib.h> /* to get working 'inline' */
71
SANITY: emits warnings when passed non-sane colours (and usually
72
corrects them) -- useful when debugging.
74
APPROX: speeds up the conversion from RGB to the colourspace by
75
assuming that the RGB values passed in are integral and definitely
78
SRGB: assumes that the RGB values being passed in (and out) are
79
destined for an sRGB-alike display device (a typical modern monitor)
80
-- if you change this then you'll probably want to change ASSUMED_GAMMA,
81
the phosphor colours and the white point definition.
93
#define ASSUMED_GAMMA (2.2F)
95
/*#define ASSUMED_GAMMA (2.591F)*/
96
#define ASSUMED_GAMMA (2.2F)
98
#define REV_GAMMA ((1.0F / ASSUMED_GAMMA))
101
/* define characteristics of the source RGB space (and the space
102
within which we try to behave linearly). */
104
/* Phosphor colours: */
106
/* sRGB/HDTV phosphor colours */
107
static const double pxr = 0.64F;
108
static const double pyr = 0.33F;
109
static const double pxg = 0.30F;
110
static const double pyg = 0.60F;
111
static const double pxb = 0.15F;
112
static const double pyb = 0.06F;
116
/* D65 (6500K) (recommended but not a common display default) */
117
static const double lxn = 0.312713F;
118
static const double lyn = 0.329016F;
121
/*static const double lxn = 0.3457F; */
122
/*static const double lyn = 0.3585F; */
125
/*static const double lxn = 0.3324F; */
126
/*static const double lyn = 0.3474F; */
128
/* D93 (9300K) (a common monitor default, but poor colour reproduction) */
129
/* static const double lxn = 0.2848F; */
130
/* static const double lyn = 0.2932F; */
132
/* illum E (normalized) */
133
/*static const double lxn = 1.0/3.0F; */
134
/*static const double lyn = 1.0/3.0F; */
136
/* illum C (average sunlight) */
137
/*static const double lxn = 0.3101F; */
138
/*static const double lyn = 0.3162F; */
140
/* illum B (direct sunlight) */
141
/*static const double lxn = 0.3484F; */
142
/*static const double lyn = 0.3516F; */
144
/* illum A (tungsten lamp) */
145
/*static const double lxn = 0.4476F; */
146
/*static const double lyn = 0.4074F; */
149
static const double LRAMP = 7.99959199F;
152
static double xnn, znn;
154
static double powtable[256];
158
#define CLAMP(x,l,u) ((x)<(l)?(l):((x)>(u)?(u):(x)))
163
init_powtable(const double gamma)
168
/* pure gamma function */
169
for (i=0; i<256; i++)
171
powtable[i] = pow((i)/255.0F, gamma);
174
/* sRGB gamma curve */
175
for (i=0; i<11 /* 0.03928 * 255 */; i++)
177
powtable[i] = (i) / (255.0F * 12.92F);
181
powtable[i] = pow( (((i) / 255.0F) + 0.055F) / 1.055F, 2.4F);
187
typedef double CMatrix[3][3];
188
typedef double CVector[3];
190
static CMatrix Mrgb_to_xyz, Mxyz_to_rgb;
193
Minvert (CMatrix src, CMatrix dest)
197
dest[0][0] = src[1][1] * src[2][2] - src[1][2] * src[2][1];
198
dest[0][1] = src[0][2] * src[2][1] - src[0][1] * src[2][2];
199
dest[0][2] = src[0][1] * src[1][2] - src[0][2] * src[1][1];
200
dest[1][0] = src[1][2] * src[2][0] - src[1][0] * src[2][2];
201
dest[1][1] = src[0][0] * src[2][2] - src[0][2] * src[2][0];
202
dest[1][2] = src[0][2] * src[1][0] - src[0][0] * src[1][2];
203
dest[2][0] = src[1][0] * src[2][1] - src[1][1] * src[2][0];
204
dest[2][1] = src[0][1] * src[2][0] - src[0][0] * src[2][1];
205
dest[2][2] = src[0][0] * src[1][1] - src[0][1] * src[1][0];
208
src[0][0] * dest[0][0] +
209
src[0][1] * dest[1][0] +
210
src[0][2] * dest[2][0];
215
g_printerr ("\n\007 XXXX det: %f\n", det);
237
init_powtable (ASSUMED_GAMMA);
240
/* ynn taken as 1.0 */
241
znn = (1.0F - (lxn + lyn)) / lyn;
253
MRC[2][0] = 1.0F - (pxr + pyr);
254
MRC[2][1] = 1.0F - (pxg + pyg);
255
MRC[2][2] = 1.0F - (pxb + pyb);
259
C1 = MRCi[0][0]*xnn + MRCi[0][1] + MRCi[0][2]*znn;
260
C2 = MRCi[1][0]*xnn + MRCi[1][1] + MRCi[1][2]*znn;
261
C3 = MRCi[2][0]*xnn + MRCi[2][1] + MRCi[2][2]*znn;
263
Mrgb_to_xyz[0][0] = MRC[0][0] * C1;
264
Mrgb_to_xyz[0][1] = MRC[0][1] * C2;
265
Mrgb_to_xyz[0][2] = MRC[0][2] * C3;
266
Mrgb_to_xyz[1][0] = MRC[1][0] * C1;
267
Mrgb_to_xyz[1][1] = MRC[1][1] * C2;
268
Mrgb_to_xyz[1][2] = MRC[1][2] * C3;
269
Mrgb_to_xyz[2][0] = MRC[2][0] * C1;
270
Mrgb_to_xyz[2][1] = MRC[2][1] * C2;
271
Mrgb_to_xyz[2][2] = MRC[2][2] * C3;
273
Minvert (Mrgb_to_xyz, Mxyz_to_rgb);
279
xyz_to_rgb (double *inx_outr,
283
const double x = *inx_outr;
284
const double y = *iny_outg;
285
const double z = *inz_outb;
287
*inx_outr = Mxyz_to_rgb[0][0]*x + Mxyz_to_rgb[0][1]*y + Mxyz_to_rgb[0][2]*z;
288
*iny_outg = Mxyz_to_rgb[1][0]*x + Mxyz_to_rgb[1][1]*y + Mxyz_to_rgb[1][2]*z;
289
*inz_outb = Mxyz_to_rgb[2][0]*x + Mxyz_to_rgb[2][1]*y + Mxyz_to_rgb[2][2]*z;
294
rgb_to_xyz (double *inr_outx,
298
const double r = *inr_outx;
299
const double g = *ing_outy;
300
const double b = *inb_outz;
302
*inr_outx = Mrgb_to_xyz[0][0]*r + Mrgb_to_xyz[0][1]*g + Mrgb_to_xyz[0][2]*b;
303
*ing_outy = Mrgb_to_xyz[1][0]*r + Mrgb_to_xyz[1][1]*g + Mrgb_to_xyz[1][2]*b;
304
*inb_outz = Mrgb_to_xyz[2][0]*r + Mrgb_to_xyz[2][1]*g + Mrgb_to_xyz[2][2]*b;
308
/* call this before using the CPercep function */
310
cpercep_init_conversions(void)
318
ffunc(const double t)
326
return (7.787F * t + 16.0F/116.0F);
332
ffunc_inv(const double t)
340
return ((t - 16.0F/116.0F) / 7.787F);
346
xyz_to_lab (double *inx,
352
const double X = *inx;
353
const double Y = *iny;
354
const double Z = *inz;
360
L = (116.0F * cbrt(Y)) - 16.0F;
370
g_printerr (" <eek1>%f \007",(float)L);
375
g_printerr (" <eek2>%f \007",(float)L);
385
a = 500.0F * (ffunc(X/xnn) - ffuncY);
386
b = 200.0F * (ffuncY - ffunc(Z/znn));
395
lab_to_xyz (double *inl,
401
const double L = *inl;
402
const double a = *ina;
403
const double b = *inb;
407
P = Y = (L + 16.0F) / 116.0F;
413
P = 7.787F * Y + 16.0F/116.0F;
416
X = (P + a / 500.0F);
417
X = xnn * ffunc_inv(X);
418
Z = (P - b / 200.0F);
419
Z = znn * ffunc_inv(Z);
425
g_printerr ("{badX %f {%f,%f,%f}}",X,L,a,b);
431
g_printerr ("{badY %f}",Y);
437
g_printerr ("{badZ %f}",Z);
451
cpercep_rgb_to_space (double inr, double ing, double inb,
452
double* outr, double* outg, double* outb)
456
/* ADM extra sanity */
457
if ((inr) > 255.0F ||
466
inr = powtable[(int)inr];
467
ing = powtable[(int)ing];
468
inb = powtable[(int)inb];
471
/* sRGB gamma curve */
472
if (inr <= (0.03928F * 255.0F))
473
inr = inr / (255.0F * 12.92F);
475
inr = pow( (inr + (0.055F * 255.0F)) / (1.055F * 255.0F), 2.4F);
477
if (ing <= (0.03928F * 255.0F))
478
ing = ing / (255.0F * 12.92F);
480
ing = pow( (ing + (0.055F * 255.0F)) / (1.055F * 255.0F), 2.4F);
482
if (inb <= (0.03928F * 255.0F))
483
inb = inb / (255.0F * 12.92F);
485
inb = pow( (inb + (0.055F * 255.0F)) / (1.055F * 255.0F), 2.4F);
487
/* pure gamma function */
488
inr = pow((inr)/255.0F, ASSUMED_GAMMA);
489
ing = pow((ing)/255.0F, ASSUMED_GAMMA);
490
inb = pow((inb)/255.0F, ASSUMED_GAMMA);
495
/* ADM extra sanity */
509
rgb_to_xyz(&inr, &ing, &inb);
512
if (inr < 0.0F || ing < 0.0F || inb < 0.0F)
514
g_printerr (" [BAD2 XYZ: %f,%f,%f]\007 ",
519
xyz_to_lab(&inr, &ing, &inb);
528
cpercep_space_to_rgb (double inr, double ing, double inb,
529
double* outr, double* outg, double* outb)
531
lab_to_xyz(&inr, &ing, &inb);
534
if (inr<-0.0F || ing<-0.0F || inb<-0.0F)
536
g_printerr (" [BAD1 XYZ: %f,%f,%f]\007 ",
541
xyz_to_rgb(&inr, &ing, &inb);
543
/* yes, essential. :( */
544
inr = CLAMP(inr,0.0F,1.0F);
545
ing = CLAMP(ing,0.0F,1.0F);
546
inb = CLAMP(inb,0.0F,1.0F);
549
if (inr <= 0.0030402477F)
550
inr = inr * (12.92F * 255.0F);
552
inr = pow(inr, 1.0F/2.4F) * (1.055F * 255.0F) - (0.055F * 255.0F);
554
if (ing <= 0.0030402477F)
555
ing = ing * (12.92F * 255.0F);
557
ing = pow(ing, 1.0F/2.4F) * (1.055F * 255.0F) - (0.055F * 255.0F);
559
if (inb <= 0.0030402477F)
560
inb = inb * (12.92F * 255.0F);
562
inb = pow(inb, 1.0F/2.4F) * (1.055F * 255.0F) - (0.055F * 255.0F);
564
inr = 255.0F * pow(inr, REV_GAMMA);
565
ing = 255.0F * pow(ing, REV_GAMMA);
566
inb = 255.0F * pow(inb, REV_GAMMA);
576
/* EXPERIMENTAL SECTION */
579
xscaler(const double start, const double end,
580
const double me, const double him)
582
return start + ((end-start) * him) / (me + him);
587
mix_colours (const double L1, const double a1, const double b1,
588
const double L2, const double a2, const double b2,
589
double *rtnL, double *rtna, double *rtnb,
590
double mass1, double mass2)
595
*rtnL = xscaler (L1, L2, mass1, mass2);
596
*rtna = xscaler (a1, a2, mass1, mass2);
597
*rtnb = xscaler (b1, b2, mass1, mass2);
604
w1 = mass1 * (L1*L1*L1);
605
w2 = mass2 * (L2*L2*L2);
608
*rtnL = xscaler (L1, L2, mass1, mass2);
616
/* g_printerr ("\007OUCH. "); */
621
*rtna = xscaler(a1, a2, w1, w2);
622
*rtnb = xscaler(b1, b2, w1, w2);
626
#endif /* EXPERIMENTAL SECTION */
added
removed
app/core/cpercep.c
