1 /* See LICENSE file for copyright and license details. */
2 #include <math.h>
3 
4 #define D65_XYY_X 0.312726871026564878786047074755
5 #define D65_XYY_Y 0.329023206641284038376227272238
6 
7 #define D65_XYZ_X (D65_XYY_X / D65_XYY_Y)
8 #define D65_XYZ_Z (1 / D65_XYY_Y - 1 - D65_XYZ_X)
9 
10 static inline double
srgb_encode(double t)11 srgb_encode(double t)
12 {
13 	double sign = 1;
14 	if (t < 0) {
15 		t = -t;
16 		sign = -1;
17 	}
18 	t = t <= 0.0031306684425217108 ? 12.92 * t : 1.055 * pow(t, 1 / 2.4) - 0.055;
19 	return t * sign;
20 }
21 
22 static inline double
srgb_decode(double t)23 srgb_decode(double t)
24 {
25 	double sign = 1;
26 	if (t < 0) {
27 		t = -t;
28 		sign = -1;
29 	}
30 	t = t <= 0.0031306684425217108 * 12.92 ? t / 12.92 : pow((t + 0.055) / 1.055, 2.4);
31 	return t * sign;
32 }
33 
34 static inline void
yuv_to_srgb(double y,double u,double v,double * r,double * g,double * b)35 yuv_to_srgb(double y, double u, double v, double *r, double *g, double *b)
36 {
37 #define MULTIPLY(CY, CU, CV)  ((CY) * y + (CU) * u + (CV) * v)
38 	*r = MULTIPLY(1,  0.00028328010485821202317155420580263580632163211703,  1.14070449590558520291949662350816652178764343261719);
39 	*g = MULTIPLY(1, -0.39630886669497211727275498560629785060882568359375, -0.58107364288228224857846271333983168005943298339844);
40 	*b = MULTIPLY(1,  2.03990003507541306504435851820744574069976806640625,  0.00017179031692307700847528739718228507626918144524);
41 #undef MULTIPLY
42 }
43 
44 static inline void
srgb_to_yuv(double r,double g,double b,double * y,double * u,double * v)45 srgb_to_yuv(double r, double g, double b, double *y, double *u, double *v)
46 {
47 #define MULTIPLY(CR, CG, CB) ((CR) * r + (CG) * g + (CB) * b)
48 	*y = MULTIPLY(0.299, 0.587, 0.114);
49 	*u = MULTIPLY(-0.14662756598240470062854967636667424812912940979004,
50 		      -0.28771586836102963635752871596196200698614120483398,
51 		       0.43434343434343436474165400795754976570606231689453);
52 	*v = MULTIPLY( 0.61456892577224520035628074765554629266262054443359,
53 		      -0.51452282157676354490405401520547457039356231689453,
54 		      -0.10004610419548178035231700278018251992762088775635);
55 #undef MULTIPLY
56 }
57 
58 static inline void
ciexyz_to_srgb(double x,double y,double z,double * r,double * g,double * b)59 ciexyz_to_srgb(double x, double y, double z, double *r, double *g, double *b)
60 {
61 #define MULTIPLY(CX, CY, CZ)  ((CX) * x + (CY) * y + (CZ) * z)
62 	*r = MULTIPLY(3.240446254647737500675930277794, -1.537134761820080575134284117667, -0.498530193022728718155178739835);
63 	*g = MULTIPLY(-0.969266606244679751469561779231, 1.876011959788370209167851498933, 0.041556042214430065351304932619);
64 	*b = MULTIPLY(0.055643503564352832235773149705, -0.204026179735960239147729566866, 1.057226567722703292062647051353);
65 #undef MULTIPLY
66 }
67 
68 static inline void
srgb_to_ciexyz(double r,double g,double b,double * x,double * y,double * z)69 srgb_to_ciexyz(double r, double g, double b, double *x, double *y, double *z)
70 {
71 #define MULTIPLY(CR, CG, CB) ((CR) * r + (CG) * g + (CB) * b)
72 	*x = MULTIPLY(0.412457445582367576708548995157, 0.357575865245515878143578447634, 0.180437247826399665973085006954);
73 	*y = MULTIPLY(0.212673370378408277403536885686, 0.715151730491031756287156895269, 0.072174899130559869164791564344);
74 	*z = MULTIPLY(0.019333942761673460208893260415, 0.119191955081838593666354597644, 0.950302838552371742508739771438);
75 #undef MULTIPLY
76 }
77 
78 static inline void
scaled_yuv_to_ciexyz(double y,double u,double v,double * xp,double * yp,double * zp)79 scaled_yuv_to_ciexyz(double y, double u, double v, double *xp, double *yp, double *zp)
80 {
81 #define MULTIPLY(CY, CU, CV) ((CY) * y + (CU) * u + (CV) * v)
82 	*xp = MULTIPLY( 0.00001450325106667098632156481796684488472237717360,
83 		        0.00000345586790639342739093228633329157872822179343,
84 		        0.00000400923398630552893485111398685916128670214675);
85 	*yp = MULTIPLY( 0.00001525902189669641837040624243737596543724066578,
86 		       -0.00000207722814409390653614547427030512238843584782,
87 		       -0.00000263898607692305410302407824019166326934282552);
88 	*zp = MULTIPLY( 0.00001661446153041708825425643025752719950105529279,
89 		        0.00002885925752619118069149627137104374696718878113,
90 			-0.00000071781086875769179526501342566979779746816348);
91 #undef MULTIPLY
92 }
93 
94 static inline void
ciexyz_to_scaled_yuv(double x,double y,double z,double * yp,double * up,double * vp)95 ciexyz_to_scaled_yuv(double x, double y, double z, double *yp, double *up, double *vp)
96 {
97 #define MULTIPLY(CX, CY, CZ) ((CX) * x + (CY) * y + (CZ) * z)
98 	*yp = MULTIPLY(  26625.38231027395886485464870929718017578125,
99 		         40524.0090949436053051613271236419677734375,
100 		          -271.5313105642117079696618020534515380859375);
101 	*up = MULTIPLY( -11278.3751445417292416095733642578125,
102 		        -26409.91773157499847002327442169189453125,
103 		         34100.5706543184860493056476116180419921875);
104 	*vp = MULTIPLY( 162829.60100012840121053159236907958984375,
105 		       -123829.313212639070115983486175537109375,
106 			-28411.65702312920984695665538311004638671875);
107 #undef MULTIPLY
108 }
109