1/*
2 * Copyright 2006, Haiku. All rights reserved.
3 * Distributed under the terms of the MIT License.
4 *
5 * Authors:
6 *		Stephan Aßmus <superstippi@gmx.de>
7 */
8
9#include <math.h>
10#include <stdio.h>
11
12#include "support.h"
13
14#include "lab_convert.h"
15
16#define GAMMA_ZERO_ENTRIES 256
17#define GAMMA_ENTRIES 10240
18#define GAMMA_MAX_ENTRIES 256
19#define GAMMA_TOTAL_ENTRIES GAMMA_ZERO_ENTRIES + GAMMA_ENTRIES + GAMMA_MAX_ENTRIES
20
21// init_gamma_table
22uint8*
23init_gamma_table()
24{
25	uint8* table = new uint8[GAMMA_TOTAL_ENTRIES];
26	for (int32 i = 0; i < GAMMA_ZERO_ENTRIES; i++)
27		table[i] = 0;
28	for (int32 i = 0; i < GAMMA_ENTRIES; i++)
29		table[i + GAMMA_ZERO_ENTRIES] = (uint8)(pow((float)i / (float)(GAMMA_ENTRIES - 1), 0.4) * 255.0 + 0.5);
30	for (int32 i = 0; i < GAMMA_MAX_ENTRIES; i++)
31		table[i + GAMMA_ZERO_ENTRIES + GAMMA_ENTRIES] = 255;
32	return table;
33}
34
35// init_linear_table
36float*
37init_linear_table()
38{
39	float* table = new float[256];
40	for (int32 i = 0; i < 256; i++)
41		table[i] = pow((float)i / 255.0, 2.5);
42	return table;
43}
44
45// conversion from RGB (0...255) to linear and normalized RGB (0...1)
46/*static uint8* gammaTable = init_gamma_table();*/
47static float* linearTable = init_linear_table();
48
49// matrix entries: XYZ -> RGB (709 RGB, D65 Whitepoint)
50/*const float Rx = 3.240479;
51const float Ry = -1.537150;
52const float Rz = -0.498535;
53const float Gx = -0.969256;
54const float Gy = 1.875992;
55const float Gz = 0.041556;
56const float Bx = 0.055648;
57const float By = -0.204043;
58const float Bz = 1.057311;*/
59
60// matrix entries: XYZ -> sRGB (D65 Whitepoint)
61const float Rx = 3.24071;
62const float Ry = -1.53726;
63const float Rz = -0.498571;
64const float Gx = -0.969258;
65const float Gy = 1.87599;
66const float Gz = 0.0415557;
67const float Bx = 0.0556352;
68const float By = -0.203996;
69const float Bz = 1.05707;
70
71
72// matrix entries scaled: XYZ(0...1) -> RGB(0...255)
73const float RX = Rx * 255;
74const float RY = Ry * 255;
75const float RZ = Rz * 255;
76const float GX = Gx * 255;
77const float GY = Gy * 255;
78const float GZ = Gz * 255;
79const float BX = Bx * 255;
80const float BY = By * 255;
81const float BZ = Bz * 255;
82
83// matrix etries: RGB -> XYZ
84/*const float Xr = 0.412453;
85const float Xg = 0.357580;
86const float Xb = 0.189423;
87const float Yr = 0.212671;
88const float Yg = 0.715160;
89const float Yb = 0.072169;
90const float Zr = 0.019334;
91const float Zg = 0.119193;
92const float Zb = 0.950227;*/
93
94// matrix etries: sRGB -> XYZ (D65 Whitepoint)
95const float Xr = 0.412424;
96const float Xg = 0.357579;
97const float Xb = 0.180464;
98const float Yr = 0.212656;
99const float Yg = 0.715158;
100const float Yb = 0.0721856;
101const float Zr = 0.0193324;
102const float Zg = 0.119193;
103const float Zb = 0.950444;
104
105// matrix entries scaled: RGB(0...255) -> XYZ(0...1)
106const float XR = Xr / 255;
107const float XG = Xg / 255;
108const float XB = Xb / 255;
109const float YR = Yr / 255;
110const float YG = Yg / 255;
111const float YB = Yb / 255;
112const float ZR = Zr / 255;
113const float ZG = Zg / 255;
114const float ZB = Zb / 255;
115
116// white point
117const float Xn = Xr + Xg + Xb;
118const float Yn = Yr + Yg + Yb;
119const float Zn = Zr + Zg + Zb;
120
121// matrix entries scaled and white point normalized: RGB(0...255) -> XYZ(0...1/WP)
122const float XRn = Xr / Xn;
123const float XGn = Xg / Xn;
124const float XBn = Xb / Xn;
125const float YRn = Yr / Yn;
126const float YGn = Yg / Yn;
127const float YBn = Yb / Yn;
128const float ZRn = Zr / Zn;
129const float ZGn = Zg / Zn;
130const float ZBn = Zb / Zn;
131
132// lab2rgb
133void
134lab2rgb(float L, float a, float b, uint8& R, uint8& G, uint8& B)
135{
136	float P = (L + 16) / 116;
137	float Pa500 = P + a / 500;
138	float Pb200 = P - b / 200;
139	float X = Xn * Pa500 * Pa500 * Pa500;
140	float Y = Yn * P * P * P;
141	float Z = Zn * Pb200 * Pb200 * Pb200;
142/*	float linearR = max_c(0.0, min_c(1.0, Rx * X + Ry * Y + Rz * Z));
143	float linearG = max_c(0.0, min_c(1.0, Gx * X + Gy * Y + Gz * Z));
144	float linearB = max_c(0.0, min_c(1.0, Bx * X + By * Y + Bz * Z));
145	R = (uint8)(pow(linearR, 0.4) * 255.0 + 0.5);
146	G = (uint8)(pow(linearG, 0.4) * 255.0 + 0.5);
147	B = (uint8)(pow(linearB, 0.4) * 255.0 + 0.5);*/
148	float linearR = Rx * X + Ry * Y + Rz * Z;
149	float linearG = Gx * X + Gy * Y + Gz * Z;
150	float linearB = Bx * X + By * Y + Bz * Z;
151	R = (uint8)constrain_int32_0_255((int32)(pow(linearR, 0.4) * 255.0 + 0.5));
152	G = (uint8)constrain_int32_0_255((int32)(pow(linearG, 0.4) * 255.0 + 0.5));
153	B = (uint8)constrain_int32_0_255((int32)(pow(linearB, 0.4) * 255.0 + 0.5));
154/*	float linearR = Rx * X + Ry * Y + Rz * Z;
155	float linearG = Gx * X + Gy * Y + Gz * Z;
156	float linearB = Bx * X + By * Y + Bz * Z;
157	R = gammaTable[(uint32)(linearR * (GAMMA_ENTRIES - 1) + 0.5) + GAMMA_ZERO_ENTRIES];
158	G = gammaTable[(uint32)(linearG * (GAMMA_ENTRIES - 1) + 0.5) + GAMMA_ZERO_ENTRIES];
159	B = gammaTable[(uint32)(linearB * (GAMMA_ENTRIES - 1) + 0.5) + GAMMA_ZERO_ENTRIES];*/
160}
161
162inline float
163f(float t)
164{
165	if (t > 0.008856)
166		return pow(t, 1.0 / 3.0);
167	return 7.787 * t + 16.0 / 116;
168}
169
170// rgb2lab
171void
172rgb2lab(uint8 R, uint8 G, uint8 B, float& L, float& a, float& b)
173{
174/*	float linearR = pow((float)R / 255.0, 2.5);
175	float linearG = pow((float)G / 255.0, 2.5);
176	float linearB = pow((float)B / 255.0, 2.5);*/
177	float linearR = linearTable[R];
178	float linearG = linearTable[G];
179	float linearB = linearTable[B];
180	float Xq = XRn * linearR + XGn * linearG + XBn * linearB; // == X/Xn
181	float Yq = YRn * linearR + YGn * linearG + YBn * linearB; // == Y/Yn
182	float Zq = ZRn * linearR + ZGn * linearG + ZBn * linearB; // == Z/Zn
183	if (Yq > 0.008856)
184		L = 116.0 * pow(Yq, 1.0 / 3.0) - 16;
185	else
186		L = 903.3 * Yq;
187	a = 500 * (f(Xq) - f(Yq));
188	b = 200 * (f(Yq) - f(Zq));
189}
190
191// replace_luminance
192void
193replace_luminance(uint8& r1, uint8& g1, uint8& b1, uint8 r2, uint8 g2, uint8 b2)
194{
195	float CIEL1, CIEa1, CIEb1, CIEL2;//, CIEa2, CIEb2;
196	rgb2lab(r1, g1, b1, CIEL1, CIEa1, CIEb1);
197//	rgb2lab(r2, g2, b2, CIEL2, CIEa2, CIEb2);
198	CIEL2 = ((linearTable[r2] + linearTable[g2] + linearTable[b2]) / 3.0) * 100.0;
199	lab2rgb(CIEL2, CIEa1, CIEb1, r1, g1, b1);
200}
201