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
59
#include <glib-object.h>
61
#include <libgimpmath/gimpmath.h>
64
/* cbrt() is a GNU extension */
65
#define cbrt(x) (pow(x, 1.0/3.0))
73
SANITY: emits warnings when passed non-sane colours (and usually
74
corrects them) -- useful when debugging.
76
APPROX: speeds up the conversion from RGB to the colourspace by
77
assuming that the RGB values passed in are integral and definitely
80
SRGB: assumes that the RGB values being passed in (and out) are
81
destined for an sRGB-alike display device (a typical modern monitor)
82
-- if you change this then you'll probably want to change ASSUMED_GAMMA,
83
the phosphor colours and the white point definition.
92
#define ASSUMED_GAMMA (2.2F)
94
/*#define ASSUMED_GAMMA (2.591F)*/
95
#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;
309
ffunc(const double t)
317
return (7.787F * t + 16.0F/116.0F);
323
ffunc_inv(const double t)
331
return ((t - 16.0F/116.0F) / 7.787F);
337
xyz_to_lab (double *inx,
343
const double X = *inx;
344
const double Y = *iny;
345
const double Z = *inz;
351
L = (116.0F * cbrt(Y)) - 16.0F;
361
g_printerr (" <eek1>%f \007",(float)L);
366
g_printerr (" <eek2>%f \007",(float)L);
376
a = 500.0F * (ffunc(X/xnn) - ffuncY);
377
b = 200.0F * (ffuncY - ffunc(Z/znn));
386
lab_to_xyz (double *inl,
392
const double L = *inl;
393
const double a = *ina;
394
const double b = *inb;
398
P = Y = (L + 16.0F) / 116.0F;
404
P = 7.787F * Y + 16.0F/116.0F;
407
X = (P + a / 500.0F);
408
X = xnn * ffunc_inv(X);
409
Z = (P - b / 200.0F);
410
Z = znn * ffunc_inv(Z);
416
g_printerr ("{badX %f {%f,%f,%f}}",X,L,a,b);
422
g_printerr ("{badY %f}",Y);
428
g_printerr ("{badZ %f}",Z);
440
/* call this before using the CPercep function */
444
static gboolean initialized = FALSE;
454
cpercep_rgb_to_space (double inr,
463
/* ADM extra sanity */
464
if ((inr) > 255.0F ||
473
inr = powtable[(int)inr];
474
ing = powtable[(int)ing];
475
inb = powtable[(int)inb];
478
/* sRGB gamma curve */
479
if (inr <= (0.03928F * 255.0F))
480
inr = inr / (255.0F * 12.92F);
482
inr = pow( (inr + (0.055F * 255.0F)) / (1.055F * 255.0F), 2.4F);
484
if (ing <= (0.03928F * 255.0F))
485
ing = ing / (255.0F * 12.92F);
487
ing = pow( (ing + (0.055F * 255.0F)) / (1.055F * 255.0F), 2.4F);
489
if (inb <= (0.03928F * 255.0F))
490
inb = inb / (255.0F * 12.92F);
492
inb = pow( (inb + (0.055F * 255.0F)) / (1.055F * 255.0F), 2.4F);
494
/* pure gamma function */
495
inr = pow((inr)/255.0F, ASSUMED_GAMMA);
496
ing = pow((ing)/255.0F, ASSUMED_GAMMA);
497
inb = pow((inb)/255.0F, ASSUMED_GAMMA);
502
/* ADM extra sanity */
516
rgb_to_xyz(&inr, &ing, &inb);
519
if (inr < 0.0F || ing < 0.0F || inb < 0.0F)
521
g_printerr (" [BAD2 XYZ: %f,%f,%f]\007 ",
526
xyz_to_lab(&inr, &ing, &inb);
535
cpercep_space_to_rgb (double inr,
542
lab_to_xyz(&inr, &ing, &inb);
545
if (inr<-0.0F || ing<-0.0F || inb<-0.0F)
547
g_printerr (" [BAD1 XYZ: %f,%f,%f]\007 ",
552
xyz_to_rgb(&inr, &ing, &inb);
554
/* yes, essential. :( */
555
inr = CLAMP(inr,0.0F,1.0F);
556
ing = CLAMP(ing,0.0F,1.0F);
557
inb = CLAMP(inb,0.0F,1.0F);
560
if (inr <= 0.0030402477F)
561
inr = inr * (12.92F * 255.0F);
563
inr = pow(inr, 1.0F/2.4F) * (1.055F * 255.0F) - (0.055F * 255.0F);
565
if (ing <= 0.0030402477F)
566
ing = ing * (12.92F * 255.0F);
568
ing = pow(ing, 1.0F/2.4F) * (1.055F * 255.0F) - (0.055F * 255.0F);
570
if (inb <= 0.0030402477F)
571
inb = inb * (12.92F * 255.0F);
573
inb = pow(inb, 1.0F/2.4F) * (1.055F * 255.0F) - (0.055F * 255.0F);
575
inr = 255.0F * pow(inr, REV_GAMMA);
576
ing = 255.0F * pow(ing, REV_GAMMA);
577
inb = 255.0F * pow(inb, REV_GAMMA);
587
/* EXPERIMENTAL SECTION */
590
xscaler(const double start, const double end,
591
const double me, const double him)
593
return start + ((end-start) * him) / (me + him);
598
mix_colours (const double L1, const double a1, const double b1,
599
const double L2, const double a2, const double b2,
600
double *rtnL, double *rtna, double *rtnb,
601
double mass1, double mass2)
606
*rtnL = xscaler (L1, L2, mass1, mass2);
607
*rtna = xscaler (a1, a2, mass1, mass2);
608
*rtnb = xscaler (b1, b2, mass1, mass2);
615
w1 = mass1 * (L1*L1*L1);
616
w2 = mass2 * (L2*L2*L2);
619
*rtnL = xscaler (L1, L2, mass1, mass2);
627
/* g_printerr ("\007OUCH. "); */
632
*rtna = xscaler(a1, a2, w1, w2);
633
*rtnb = xscaler(b1, b2, w1, w2);
637
#endif /* EXPERIMENTAL SECTION */
added
removed
app/base/cpercep.c
