1 /* babl - dynamically extendable universal pixel conversion library.
2  * Copyright (C) 2005, 2014, 2019 Øyvind Kolås.
3  * Copyright (C) 2009, Martin Nordholts
4  * Copyright (C) 2014, 2019 Elle Stone
5  * Copyright (C) 2017, 2018 Red Hat, Inc.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 3 of the License, or (at your option) any later version.
11  *
12  * This library is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General
18  * Public License along with this library; if not, see
19  * <https://www.gnu.org/licenses/>.
20  */
21 
22 #include "config.h"
23 #include <math.h>
24 #include <stdint.h>
25 #include <string.h>
26 
27 #if defined(USE_SSE2)
28 #include <emmintrin.h>
29 #endif /* defined(USE_SSE2) */
30 
31 #include "babl-internal.h"
32 #include "extensions/util.h"
33 
34 #define DEGREES_PER_RADIAN (180 / 3.14159265358979323846)
35 #define RADIANS_PER_DEGREE (1 / DEGREES_PER_RADIAN)
36 
37 #define LAB_EPSILON       (216.0f / 24389.0f)
38 #define LAB_KAPPA         (24389.0f / 27.0f)
39 
40 /* The constants below hard-code the D50-adapted sRGB ICC profile
41  * reference white, aka the ICC profile D50 illuminant.
42  *
43  * In a properly ICC profile color-managed application, the profile
44  * illuminant values should be retrieved from the image's
45  * ICC profile's illuminant.
46  *
47  * At present, the ICC profile illuminant is always D50. This might
48  * change when the next version of the ICC specs is released.
49  *
50  * As encoded in an actual V2 or V4 ICC profile,
51  * the illuminant values are hexadecimal-rounded, as are the following
52  * hard-coded D50 ICC profile illuminant values:
53  */
54 
55 #define D50_WHITE_REF_X   0.964202880f
56 #define D50_WHITE_REF_Y   1.000000000f
57 #define D50_WHITE_REF_Z   0.824905400f
58 
59 #define NEAR_ZERO         0.0000000001f
60 #define near_zero(a)   ((a) < NEAR_ZERO && (a) > -NEAR_ZERO)
61 
62 #define D50_WHITE_REF_x   0.345702921222f
63 #define D50_WHITE_REF_y   0.358537532290f
64 
65 
66 static void types (void);
67 static void components (void);
68 static void models (void);
69 static void conversions (void);
70 static void formats (void);
71 
72 int init (void);
73 
74 int
init(void)75 init (void)
76 {
77   types ();
78   components ();
79   models ();
80   formats ();
81   conversions ();
82   return 0;
83 }
84 
85 static void
components(void)86 components (void)
87 {
88   babl_component_new ("CIE L", "doc", "Luminance, range 0.0-100.0 in float", NULL);
89   babl_component_new ("CIE a", "chroma", "doc", "chroma component 0.0 is no saturation", NULL);
90   babl_component_new ("CIE b", "chroma", "doc", "chroma component 0.0 is no saturation", NULL);
91   babl_component_new ("CIE C(ab)", "chroma", "doc", "chrominance/saturation", NULL);
92   babl_component_new ("CIE H(ab)", "chroma", "doc", "hue value range 0.0-360.0", NULL);
93   babl_component_new ("CIE X", NULL);
94   babl_component_new ("CIE Y", NULL);
95   babl_component_new ("CIE Z", NULL);
96   babl_component_new ("CIE x", NULL);
97   babl_component_new ("CIE y", NULL);
98   /* CIE 1976 UCS */
99   babl_component_new ("CIE u", NULL);
100   babl_component_new ("CIE v", NULL);
101 /*  babl_component_new ("CIE z", NULL);*/
102 }
103 
104 static void
models(void)105 models (void)
106 {
107   babl_model_new (
108     "name", "CIE Lab",
109     "doc", "CIE Lab color model, a perceptually uniform space, euclidian distance in this space represents delta E.",
110     babl_component ("CIE L"),
111     babl_component ("CIE a"),
112     babl_component ("CIE b"),
113     "CIE",
114     NULL);
115 
116   babl_model_new (
117     "name", "CIE Lab alpha",
118     "doc", "CIE Lab color model, with separate alpha",
119     babl_component ("CIE L"),
120     babl_component ("CIE a"),
121     babl_component ("CIE b"),
122     babl_component ("A"),
123     "CIE",
124     "alpha",
125     NULL);
126 
127   babl_model_new (
128     "name", "CIE LCH(ab)",
129     "doc", "CIE LCH color model, using cylindrical coordinates",
130     babl_component ("CIE L"),
131     babl_component ("CIE C(ab)"),
132     babl_component ("CIE H(ab)"),
133     "CIE",
134     NULL);
135 
136   babl_model_new (
137     "name", "CIE LCH(ab) alpha",
138     "doc", "CIE LCH color model, using cylindrical coordinates, with separate alpha",
139     babl_component ("CIE L"),
140     babl_component ("CIE C(ab)"),
141     babl_component ("CIE H(ab)"),
142     babl_component ("A"),
143     "CIE",
144     "alpha",
145     NULL);
146 
147   babl_model_new (
148     "name", "CIE XYZ",
149     babl_component ("CIE X"),
150     babl_component ("CIE Y"),
151     babl_component ("CIE Z"),
152     "CIE",
153     NULL);
154 
155   babl_model_new (
156     "name", "CIE XYZ alpha",
157     babl_component ("CIE X"),
158     babl_component ("CIE Y"),
159     babl_component ("CIE Z"),
160     babl_component ("A"),
161     "CIE",
162     "alpha",
163     NULL);
164 
165   babl_model_new (
166     "name", "CIE xyY",
167     "doc", "the coordinate system often used for drawing chromaticity diagrams. Y is luminance.",
168     babl_component ("CIE x"),
169     babl_component ("CIE y"),
170     babl_component ("CIE Y"),
171     "CIE",
172     NULL);
173 
174   babl_model_new (
175     "name", "CIE xyY alpha",
176     "doc", "the coordinate system often used for drawing chromaticity diagrams. Y is luminance, with separate alpha",
177     babl_component ("CIE x"),
178     babl_component ("CIE y"),
179     babl_component ("CIE Y"),
180     babl_component ("A"),
181     "CIE",
182     "alpha",
183     NULL);
184 
185 /* CIE 1976 UCS */
186   babl_model_new (
187     "name", "CIE Yuv",
188     "doc", "A newer more perceptually uniform space than xyY for chromaticity diagrams.",
189     babl_component ("CIE Y"),
190     babl_component ("CIE u"),
191     babl_component ("CIE v"),
192     NULL);
193 
194   babl_model_new (
195     "name", "CIE Yuv alpha",
196     "doc", "A newer more perceptually uniform space than xyY for chromaticity diagrams, with separate alpha.",
197     babl_component ("CIE Y"),
198     babl_component ("CIE u"),
199     babl_component ("CIE v"),
200     babl_component ("A"),
201     NULL);
202 }
203 
204 static void  rgbcie_init (void);
205 
206 /******** begin double RGB/CIE color space conversions ****************/
207 
208 static inline void  ab_to_CHab    (double  a,
209                                    double  b,
210                                    double *to_C,
211                                    double *to_H);
212 
213 static inline void  CHab_to_ab    (double  C,
214                                    double  H,
215                                    double *to_a,
216                                    double *to_b);
217 
218 static inline void XYZ_to_LAB     (double  X,
219                                    double  Y,
220                                    double  Z,
221                                    double *to_L,
222                                    double *to_a,
223                                    double *to_b
224                                    );
225 
226 static inline void LAB_to_XYZ     (double  L,
227                                    double  a,
228                                    double  b,
229                                    double *to_X,
230                                    double *to_Y,
231                                    double *to_Z
232                                    );
233 
234 static inline void XYZ_to_xyY     (double  X,
235                                    double  Y,
236                                    double  Z,
237                                    double *to_x,
238                                    double *to_y,
239                                    double *to_Y
240                                    );
241 
242 static inline void xyY_to_XYZ     (double  x,
243                                    double  y,
244                                    double  Y,
245                                    double *to_X,
246                                    double *to_Y,
247                                    double *to_Z
248                                    );
249 
250 /* CIE 1976 UCS */
251 static inline void XYZ_to_Yuv     (double X,
252                                    double Y,
253                                    double Z,
254                                    double *to_Y,
255                                    double *to_u,
256                                    double *to_v
257                                    );
258 
259 static inline void Yuv_to_XYZ     (double Y,
260                                    double u,
261                                    double v,
262                                    double *to_X,
263                                    double *to_Y,
264                                    double *to_Z
265                                    );
266 
267 static inline void
XYZ_to_LAB(double X,double Y,double Z,double * to_L,double * to_a,double * to_b)268 XYZ_to_LAB (double  X,
269             double  Y,
270             double  Z,
271             double *to_L,
272             double *to_a,
273             double *to_b)
274 {
275   double xr = X / D50_WHITE_REF_X;
276   double yr = Y / D50_WHITE_REF_Y;
277   double zr = Z / D50_WHITE_REF_Z;
278 
279   double fx = xr > LAB_EPSILON ? cbrt (xr) : (LAB_KAPPA * xr + 16.0) / 116.0;
280   double fy = yr > LAB_EPSILON ? cbrt (yr) : (LAB_KAPPA * yr + 16.0) / 116.0;
281   double fz = zr > LAB_EPSILON ? cbrt (zr) : (LAB_KAPPA * zr + 16.0) / 116.0;
282 
283   *to_L = 116.0 * fy - 16.0;
284   *to_a = 500.0 * (fx - fy);
285   *to_b = 200.0 * (fy - fz);
286 }
287 
288 static inline void
LAB_to_XYZ(double L,double a,double b,double * to_X,double * to_Y,double * to_Z)289 LAB_to_XYZ (double  L,
290             double  a,
291             double  b,
292             double *to_X,
293             double *to_Y,
294             double *to_Z)
295 {
296   double fy = (L + 16.0) / 116.0;
297   double fy_cubed = fy * fy * fy;
298 
299   double fx = fy + a / 500.0;
300   double fx_cubed = fx * fx * fx;
301 
302   double fz = fy - b / 200.0;
303   double fz_cubed = fz * fz * fz;
304 
305   double yr = L > LAB_KAPPA * LAB_EPSILON ? fy_cubed : L / LAB_KAPPA;
306   double xr = fx_cubed > LAB_EPSILON ? fx_cubed : (fx * 116.0 - 16.0) / LAB_KAPPA;
307   double zr = fz_cubed > LAB_EPSILON ? fz_cubed : (fz * 116.0 - 16.0) / LAB_KAPPA;
308 
309   *to_X = xr * D50_WHITE_REF_X;
310   *to_Y = yr * D50_WHITE_REF_Y;
311   *to_Z = zr * D50_WHITE_REF_Z;
312 }
313 
314 /* CIE xyY */
315 static inline void
XYZ_to_xyY(double X,double Y,double Z,double * to_x,double * to_y,double * to_Y)316 XYZ_to_xyY (double  X,
317             double  Y,
318             double  Z,
319             double *to_x,
320             double *to_y,
321             double *to_Y)
322 {
323    double sum = X + Y + Z;
324    if (near_zero (sum))
325 	{ *to_Y = 0.0;
326 	  *to_x = D50_WHITE_REF_x;
327 	  *to_y = D50_WHITE_REF_y;
328 	}
329 	else
330 	{
331 	*to_x = X / sum;
332 	*to_y = Y / sum;
333 	*to_Y = Y;
334     }
335 }
336 
337 static inline void
xyY_to_XYZ(double x,double y,double Y,double * to_X,double * to_Y,double * to_Z)338 xyY_to_XYZ (double  x,
339             double  y,
340             double  Y,
341             double *to_X,
342             double *to_Y,
343             double *to_Z)
344 {
345   if (near_zero (Y))
346     {
347 	  *to_X = 0.0;
348 	  *to_Y = 0.0;
349 	  *to_Z = 0.0;
350 	}
351   else
352 	{
353 	*to_X = (x * Y) / y;
354 	*to_Y = Y;
355 	*to_Z = ((1 - x - y) * Y) / y;
356     }
357 }
358 
359 
360 /* CIE 1976 UCS */
361 
362 /* Code for near-zero XYZ and RGB for CIE 1976 UCS patterned
363  * after code written by and used with kind permission of Graeme Gill.
364  * Graeme Gill's code is in "icc.c"
365  * downloadable from http://argyllcms.com/ */
366 
367 static inline void
XYZ_to_Yuv(double X,double Y,double Z,double * to_Y,double * to_u,double * to_v)368 XYZ_to_Yuv (double X,
369             double Y,
370             double Z,
371             double *to_Y,
372             double *to_u,
373             double *to_v)
374 {
375   double sum = X + (15.0 * Y) + (3.0 * Z);
376   if (near_zero (sum))
377     {
378         *to_Y = 0.0;
379         *to_u = 4.0/19.0;
380         *to_v = 9.0/19.0;
381     }
382         else
383     {
384         *to_Y = Y;
385         *to_u = (4.0 * X) / sum;
386         *to_v = (9.0 * Y) / sum;
387     }
388 }
389 
390 static inline void
Yuv_to_XYZ(double Y,double u,double v,double * to_X,double * to_Y,double * to_Z)391 Yuv_to_XYZ (double Y,
392             double u,
393             double v,
394             double *to_X,
395             double *to_Y,
396             double *to_Z)
397 {
398   if (near_zero (v))
399     {
400       *to_X = 0.0;
401       *to_Y = 0.0;
402       *to_Z = 0.0;
403     }
404   else
405     {
406       *to_X = ((9.0 * u * Y)/(4.0 * v));
407       *to_Y = Y;
408       *to_Z = -(((20.0 * v + 3.0 * u - 12.0) * Y)/(4.0 * v));
409     }
410 }
411 
412 /* rgb <-> XYZ */
413 
414 static void
rgba_to_xyz(const Babl * conversion,char * src,char * dst,long n)415 rgba_to_xyz (const Babl *conversion,
416              char       *src,
417              char       *dst,
418              long        n)
419 {
420   const Babl *space = babl_conversion_get_source_space (conversion);
421   while (n--)
422     {
423       double RGB[3]  = {((double *) src)[0],
424                         ((double *) src)[1],
425                         ((double *) src)[2] };
426       babl_space_to_xyz (space, RGB, (double*)dst);
427 
428       src += sizeof (double) * 4;
429       dst += sizeof (double) * 3;
430     }
431 }
432 
433 static void
xyz_to_rgba(const Babl * conversion,char * src,char * dst,long n)434 xyz_to_rgba (const Babl *conversion,
435              char       *src,
436              char       *dst,
437              long        n)
438 {
439   const Babl *space = babl_conversion_get_destination_space (conversion);
440   while (n--)
441     {
442       babl_space_from_xyz (space, (double*)src, (double*) dst);
443       ((double *) dst)[3] = 1.0;
444       src += sizeof (double) * 3;
445       dst += sizeof (double) * 4;
446     }
447 }
448 
449 static void
rgba_to_xyza(const Babl * conversion,char * src,char * dst,long n)450 rgba_to_xyza (const Babl *conversion,char *src,
451               char *dst,
452               long  n)
453 {
454   const Babl *space = babl_conversion_get_source_space (conversion);
455   while (n--)
456     {
457       babl_space_to_xyz (space, (double*)src, (double*)dst);
458       ((double *) dst)[3] = ((double *) src)[3];
459 
460       src += sizeof (double) * 4;
461       dst += sizeof (double) * 4;
462     }
463 }
464 
465 static void
xyza_to_rgba(const Babl * conversion,char * src,char * dst,long n)466 xyza_to_rgba (const Babl *conversion,char *src,
467               char *dst,
468               long  n)
469 {
470   const Babl *space = babl_conversion_get_destination_space (conversion);
471   while (n--)
472     {
473       babl_space_from_xyz (space, (double*)src, (double*) dst);
474       ((double *) dst)[3] = ((double *) src)[3];
475 
476       src += sizeof (double) * 4;
477       dst += sizeof (double) * 4;
478     }
479 }
480 
481 
482 /* rgb -> xyY */
483 static void
rgba_to_xyY(const Babl * conversion,char * src,char * dst,long n)484 rgba_to_xyY (const Babl *conversion,
485              char       *src,
486              char       *dst,
487              long        n)
488 {
489   const Babl *space = babl_conversion_get_source_space (conversion);
490   while (n--)
491     {
492       double XYZ[3], x, y, Y;
493 
494       babl_space_to_xyz (space, (double*)src, XYZ);
495       XYZ_to_xyY (XYZ[0], XYZ[1], XYZ[2], &x, &y, &Y);
496 
497       ((double *) dst)[0] = x;
498       ((double *) dst)[1] = y;
499       ((double *) dst)[2] = Y;
500 
501       src += sizeof (double) * 4;
502       dst += sizeof (double) * 3;
503     }
504 }
505 
506 static void
rgba_to_xyYa(const Babl * conversion,char * src,char * dst,long n)507 rgba_to_xyYa (const Babl *conversion,char *src,
508               char *dst,
509               long  n)
510 {
511   const Babl *space = babl_conversion_get_source_space (conversion);
512   while (n--)
513     {
514       double alpha = ((double *) src)[3];
515       double XYZ[3], x, y, Y;
516 
517       //convert RGB to XYZ
518       babl_space_to_xyz (space, (double*)src, XYZ);
519 
520       //convert XYZ to xyY
521       XYZ_to_xyY (XYZ[0], XYZ[1], XYZ[2], &x, &y, &Y);
522 
523       ((double *) dst)[0] = x;
524       ((double *) dst)[1] = y;
525       ((double *) dst)[2] = Y;
526       ((double *) dst)[3] = alpha;
527 
528       src += sizeof (double) * 4;
529       dst += sizeof (double) * 4;
530     }
531 }
532 
533 
534 
535 /* rgb -> Yuv */
536 static void
rgba_to_Yuv(const Babl * conversion,char * src,char * dst,long n)537 rgba_to_Yuv (const Babl *conversion,char *src,
538              char *dst,
539              long  n)
540 {
541   const Babl *space = babl_conversion_get_source_space (conversion);
542   while (n--)
543     {
544       double XYZ[3], Y, u, v;
545 
546       babl_space_to_xyz (space, (double*)src, XYZ);
547       XYZ_to_Yuv (XYZ[0], XYZ[1], XYZ[2], &Y, &u, &v);
548 
549       ((double *) dst)[0] = Y;
550       ((double *) dst)[1] = u;
551       ((double *) dst)[2] = v;
552 
553       src += sizeof (double) * 4;
554       dst += sizeof (double) * 3;
555     }
556 }
557 
558 static void
rgba_to_Yuva(const Babl * conversion,char * src,char * dst,long n)559 rgba_to_Yuva (const Babl *conversion,char *src,
560               char *dst,
561               long  n)
562 {
563   const Babl *space = babl_conversion_get_source_space (conversion);
564   while (n--)
565     {
566       double alpha = ((double *) src)[3];
567       double XYZ[3], Y, u, v;
568 
569       //convert RGB to XYZ
570       babl_space_to_xyz (space, (double*)src, XYZ);
571 
572       //convert XYZ to Yuv
573       XYZ_to_Yuv (XYZ[0], XYZ[1], XYZ[2], &Y, &u, &v);
574 
575       ((double *) dst)[0] = Y;
576       ((double *) dst)[1] = u;
577       ((double *) dst)[2] = v;
578       ((double *) dst)[3] = alpha;
579 
580       src += sizeof (double) * 4;
581       dst += sizeof (double) * 4;
582     }
583 }
584 
585 
586 /* rgbf <->xyYf */
587 static void
rgbaf_to_xyYaf(const Babl * conversion,float * src,float * dst,long samples)588 rgbaf_to_xyYaf (const Babl *conversion,
589                 float *src,
590                 float *dst,
591                 long   samples)
592 {
593   const Babl *space = babl_conversion_get_source_space (conversion);
594   float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
595   float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
596   float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
597   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
598   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
599   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
600   float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
601   float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
602   float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
603   long n = samples;
604 
605   while (n--)
606     {
607       float x, y, X, Y, Z, r, g, b, a;
608       r = src[0];
609       g = src[1];
610       b = src[2];
611       a = src[3];
612 
613       if (near_zero(r) && near_zero(g) && near_zero(b))
614         {
615           Y = 0.0f;
616           x = D50_WHITE_REF_x;
617           y = D50_WHITE_REF_y;
618         }
619       else
620         {
621           X = m_0_0 * r + m_0_1 * g + m_0_2 * b;
622           Y = m_1_0 * r + m_1_1 * g + m_1_2 * b;
623           Z = m_2_0 * r + m_2_1 * g + m_2_2 * b;
624 
625           x = X / (X + Y + Z);
626           y = Y / (X + Y + Z);
627         }
628 
629       dst[0] = x;
630       dst[1] = y;
631       dst[2] = Y;
632       dst[3] = a;
633 
634       src += 4;
635       dst += 4;
636     }
637 }
638 
639 static void
rgbf_to_xyYf(const Babl * conversion,float * src,float * dst,long samples)640 rgbf_to_xyYf (const Babl *conversion,float *src,
641               float *dst,
642               long   samples)
643 {
644   const Babl *space = babl_conversion_get_source_space (conversion);
645   float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
646   float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
647   float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
648   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
649   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
650   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
651   float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
652   float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
653   float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
654   long n = samples;
655 
656   while (n--)
657     {
658   float x, y, X, Y, Z, r, g, b;
659       r = src[0];
660       g = src[1];
661       b = src[2];
662 
663       if (near_zero(r) && near_zero(g) && near_zero(b))
664         {
665           Y = 0.0f;
666           x = D50_WHITE_REF_x;
667           y = D50_WHITE_REF_y;
668         }
669       else
670         {
671           X = m_0_0 * r + m_0_1 * g + m_0_2 * b;
672           Y = m_1_0 * r + m_1_1 * g + m_1_2 * b;
673           Z = m_2_0 * r + m_2_1 * g + m_2_2 * b;
674 
675           x = X / (X + Y + Z);
676           y = Y / (X + Y + Z);
677         }
678 
679       dst[0] = x;
680       dst[1] = y;
681       dst[2] = Y;
682 
683       src += 3;
684       dst += 3;
685     }
686 }
687 
688 
689 static void
rgbaf_to_xyYf(const Babl * conversion,float * src,float * dst,long samples)690 rgbaf_to_xyYf (const Babl *conversion,
691                float      *src,
692                float      *dst,
693                long        samples)
694 {
695   const Babl *space = babl_conversion_get_source_space (conversion);
696   float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
697   float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
698   float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
699   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
700   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
701   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
702   float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
703   float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
704   float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
705   long n = samples;
706 
707   while (n--)
708     {
709   float x, y, X, Y, Z, r, g, b;
710       r = src[0];
711       g = src[1];
712       b = src[2];
713 
714       if (near_zero(r) && near_zero(g) && near_zero(b))
715         {
716           Y = 0.0f;
717           x = D50_WHITE_REF_x;
718           y = D50_WHITE_REF_y;
719         }
720       else
721         {
722           X = m_0_0 * r + m_0_1 * g + m_0_2 * b;
723           Y = m_1_0 * r + m_1_1 * g + m_1_2 * b;
724           Z = m_2_0 * r + m_2_1 * g + m_2_2 * b;
725 
726           x = X / (X + Y + Z);
727           y = Y / (X + Y + Z);
728         }
729 
730       dst[0] = x;
731       dst[1] = y;
732       dst[2] = Y;
733 
734       src += 4;
735       dst += 3;
736     }
737 }
738 
739 
740 
741 /* rgbf -> Yuv */
742 static void
rgbaf_to_Yuvaf(const Babl * conversion,float * src,float * dst,long samples)743 rgbaf_to_Yuvaf (const Babl *conversion,
744                 float *src,
745                 float *dst,
746                 long   samples)
747 {
748   const Babl *space = babl_conversion_get_source_space (conversion);
749   float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
750   float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
751   float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
752   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
753   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
754   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
755   float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
756   float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
757   float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
758   long n = samples;
759 
760   while (n--)
761     {
762       float u, v, X, Y, Z, r, g, b, a, sum;
763       r = src[0];
764       g = src[1];
765       b = src[2];
766       a = src[3];
767 
768       if (near_zero(r) && near_zero(g) && near_zero(b))
769         {
770           Y = 0.0f;
771 		  u = 4.0/19.0;
772 		  v = 9.0/19.0;
773         }
774       else
775         {
776           X = m_0_0 * r + m_0_1 * g + m_0_2 * b;
777           Y = m_1_0 * r + m_1_1 * g + m_1_2 * b;
778           Z = m_2_0 * r + m_2_1 * g + m_2_2 * b;
779 
780 	      sum = (X + 15.0 * Y + 3.0 * Z);
781 	      u = (4.0 * X) / sum;
782 	      v = (9.0 * Y) / sum;
783         }
784 
785       dst[0] = Y;
786       dst[1] = u;
787       dst[2] = v;
788       dst[3] = a;
789 
790       src += 4;
791       dst += 4;
792     }
793 }
794 
795 static void
rgbf_to_Yuvf(const Babl * conversion,float * src,float * dst,long samples)796 rgbf_to_Yuvf (const Babl *conversion,float *src,
797               float *dst,
798               long   samples)
799 {
800   const Babl *space = babl_conversion_get_source_space (conversion);
801   float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
802   float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
803   float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
804   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
805   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
806   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
807   float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
808   float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
809   float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
810   long n = samples;
811 
812   while (n--)
813     {
814   float u, v, X, Y, Z, r, g, b, sum;
815       r = src[0];
816       g = src[1];
817       b = src[2];
818 
819       if (near_zero(r) && near_zero(g) && near_zero(b))
820         {
821           Y = 0.0f;
822 		  u = 4.0/19.0;
823 		  v = 9.0/19.0;
824         }
825       else
826         {
827           X = m_0_0 * r + m_0_1 * g + m_0_2 * b;
828           Y = m_1_0 * r + m_1_1 * g + m_1_2 * b;
829           Z = m_2_0 * r + m_2_1 * g + m_2_2 * b;
830 
831 	      sum = (X + 15.0 * Y + 3.0 * Z);
832 	      u = (4.0 * X) / sum;
833 	      v = (9.0 * Y) / sum;
834         }
835 
836       dst[0] = Y;
837       dst[1] = u;
838       dst[2] = v;
839 
840       src += 3;
841       dst += 3;
842     }
843 }
844 
845 
846 static void
rgbaf_to_Yuvf(const Babl * conversion,float * src,float * dst,long samples)847 rgbaf_to_Yuvf (const Babl *conversion,
848                float      *src,
849                float      *dst,
850                long        samples)
851 {
852   const Babl *space = babl_conversion_get_source_space (conversion);
853   float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
854   float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
855   float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
856   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
857   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
858   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
859   float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
860   float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
861   float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
862   long n = samples;
863 
864   while (n--)
865     {
866   float u, v, X, Y, Z, r, g, b, sum;
867       r = src[0];
868       g = src[1];
869       b = src[2];
870 
871       if (near_zero(r) && near_zero(g) && near_zero(b))
872         {
873           Y = 0.0f;
874 		  u = 4.0/19.0;
875 		  v = 9.0/19.0;
876         }
877       else
878         {
879           X = m_0_0 * r + m_0_1 * g + m_0_2 * b;
880           Y = m_1_0 * r + m_1_1 * g + m_1_2 * b;
881           Z = m_2_0 * r + m_2_1 * g + m_2_2 * b;
882 
883 	      sum = (X + 15.0 * Y + 3.0 * Z);
884 	      u = (4.0 * X) / sum;
885 	      v = (9.0 * Y) / sum;
886         }
887 
888       dst[0] = Y;
889       dst[1] = u;
890       dst[2] = v;
891 
892       src += 4;
893       dst += 3;
894     }
895 }
896 
897 
898 
899 
900 /* xyY -> rgb */
901 
902 static void
xyY_to_rgba(const Babl * conversion,char * src,char * dst,long n)903 xyY_to_rgba (const Babl *conversion,
904              char       *src,
905              char       *dst,
906              long        n)
907 {
908   const Babl *space = babl_conversion_get_destination_space (conversion);
909   while (n--)
910     {
911       double x = ((double *) src)[0];
912       double y = ((double *) src)[1];
913       double Y = ((double *) src)[2];
914 
915       double R, G, B, X, Z;
916 
917       //convert xyY to XYZ
918       xyY_to_XYZ (x, y, Y, &X, &Y, &Z);
919 
920       //convert XYZ to RGB
921       {
922         double XYZ[3]  = {X,Y,Z};
923         double RGB[3];
924         babl_space_from_xyz (space, XYZ, RGB);
925         R = RGB[0];
926         G = RGB[1];
927         B = RGB[2];
928       }
929 
930       ((double *) dst)[0] = R;
931       ((double *) dst)[1] = G;
932       ((double *) dst)[2] = B;
933       ((double *) dst)[3] = 1.0;
934 
935       src += sizeof (double) * 3;
936       dst += sizeof (double) * 4;
937     }
938 }
939 
940 
941 static void
xyYa_to_rgba(const Babl * conversion,char * src,char * dst,long n)942 xyYa_to_rgba (const Babl *conversion,char *src,
943               char *dst,
944               long  n)
945 {
946   const Babl *space = babl_conversion_get_destination_space (conversion);
947   while (n--)
948     {
949       double x     = ((double *) src)[0];
950       double y     = ((double *) src)[1];
951       double Y     = ((double *) src)[2];
952       double alpha = ((double *) src)[3];
953 
954       double X, Z;
955 
956       //convert xyY to XYZ
957       xyY_to_XYZ (x, y, Y, &X, &Y, &Z);
958 
959       {
960         //convert XYZ to RGB
961         double XYZ[3]  = {X,Y,Z};
962         babl_space_from_xyz (space, XYZ, (double*)dst);
963       }
964       ((double *) dst)[3] = alpha;
965 
966       src += sizeof (double) * 4;
967       dst += sizeof (double) * 4;
968     }
969 }
970 
971 
972 
973 /* Yuv -> rgb */
974 
975 static void
Yuv_to_rgba(const Babl * conversion,char * src,char * dst,long n)976 Yuv_to_rgba (const Babl *conversion,char *src,
977              char *dst,
978              long  n)
979 {
980   const Babl *space = babl_conversion_get_destination_space (conversion);
981   while (n--)
982     {
983       double Y = ((double *) src)[0];
984       double u = ((double *) src)[1];
985       double v = ((double *) src)[2];
986 
987       double R, G, B, X, Z;
988 
989       //convert Yuv to XYZ
990       Yuv_to_XYZ (Y, u, v, &X, &Y, &Z);
991 
992       //convert XYZ to RGB
993       {
994         double XYZ[3]  = {X,Y,Z};
995         double RGB[3];
996         babl_space_from_xyz (space, XYZ, RGB);
997         R = RGB[0];
998         G = RGB[1];
999         B = RGB[2];
1000       }
1001 
1002       ((double *) dst)[0] = R;
1003       ((double *) dst)[1] = G;
1004       ((double *) dst)[2] = B;
1005       ((double *) dst)[3] = 1.0;
1006 
1007       src += sizeof (double) * 3;
1008       dst += sizeof (double) * 4;
1009     }
1010 }
1011 
1012 
1013 static void
Yuva_to_rgba(const Babl * conversion,char * src,char * dst,long n)1014 Yuva_to_rgba (const Babl *conversion,char *src,
1015               char *dst,
1016               long  n)
1017 {
1018   const Babl *space = babl_conversion_get_destination_space (conversion);
1019   while (n--)
1020     {
1021       double Y     = ((double *) src)[0];
1022       double u     = ((double *) src)[1];
1023       double v     = ((double *) src)[2];
1024       double alpha = ((double *) src)[3];
1025 
1026       double X, Z;
1027 
1028       //convert Yuv to XYZ
1029       Yuv_to_XYZ (Y, u, v, &X, &Y, &Z);
1030 
1031       {
1032         //convert XYZ to RGB
1033         double XYZ[3]  = {X,Y,Z};
1034         babl_space_from_xyz (space, XYZ, (double*)dst);
1035       }
1036       ((double *) dst)[3] = alpha;
1037 
1038       src += sizeof (double) * 4;
1039       dst += sizeof (double) * 4;
1040     }
1041 }
1042 
1043 
1044 /* xyYf -> rgbf */
1045 
1046 static void
xyYf_to_rgbf(const Babl * conversion,float * src,float * dst,long samples)1047 xyYf_to_rgbf (const Babl *conversion,float *src,
1048                 float *dst,
1049                 long   samples)
1050 {
1051   const Babl *space = babl_conversion_get_source_space (conversion);
1052   float m_0_0 = space->space.XYZtoRGBf[0] * D50_WHITE_REF_X;
1053   float m_0_1 = space->space.XYZtoRGBf[1] * D50_WHITE_REF_Y;
1054   float m_0_2 = space->space.XYZtoRGBf[2] * D50_WHITE_REF_Z;
1055   float m_1_0 = space->space.XYZtoRGBf[3] * D50_WHITE_REF_X;
1056   float m_1_1 = space->space.XYZtoRGBf[4] * D50_WHITE_REF_Y;
1057   float m_1_2 = space->space.XYZtoRGBf[5] * D50_WHITE_REF_Z;
1058   float m_2_0 = space->space.XYZtoRGBf[6] * D50_WHITE_REF_X;
1059   float m_2_1 = space->space.XYZtoRGBf[7] * D50_WHITE_REF_Y;
1060   float m_2_2 = space->space.XYZtoRGBf[8] * D50_WHITE_REF_Z;
1061   long n = samples;
1062 
1063   while (n--)
1064     {
1065       float X, Z, r, g, b;
1066       float x = src[0];
1067       float y = src[1];
1068       float Y = src[2];
1069 
1070       if (near_zero (y))
1071         {
1072           X = 0.0f;
1073           Y = 0.0f;
1074           Z = 0.0f;
1075         }
1076       else
1077         {
1078           X = (x * Y) / y;
1079           //Y = Y;
1080           Z = ((1 - x - y) * Y) / y;
1081         }
1082 
1083       r = m_0_0 * X + m_0_1 * Y + m_0_2 * Z;
1084       g = m_1_0 * X + m_1_1 * Y + m_1_2 * Z;
1085       b = m_2_0 * X + m_2_1 * Y + m_2_2 * Z;
1086 
1087       dst[0] = r;
1088       dst[1] = g;
1089       dst[2] = b;
1090 
1091       src += 3;
1092       dst += 3;
1093     }
1094 }
1095 
1096 
1097 
1098 static void
xyYf_to_rgbaf(const Babl * conversion,float * src,float * dst,long samples)1099 xyYf_to_rgbaf (const Babl *conversion,
1100                float      *src,
1101                float      *dst,
1102                long        samples)
1103 {
1104   const Babl *space = babl_conversion_get_source_space (conversion);
1105   float m_0_0 = space->space.XYZtoRGBf[0] * D50_WHITE_REF_X;
1106   float m_0_1 = space->space.XYZtoRGBf[1] * D50_WHITE_REF_Y;
1107   float m_0_2 = space->space.XYZtoRGBf[2] * D50_WHITE_REF_Z;
1108   float m_1_0 = space->space.XYZtoRGBf[3] * D50_WHITE_REF_X;
1109   float m_1_1 = space->space.XYZtoRGBf[4] * D50_WHITE_REF_Y;
1110   float m_1_2 = space->space.XYZtoRGBf[5] * D50_WHITE_REF_Z;
1111   float m_2_0 = space->space.XYZtoRGBf[6] * D50_WHITE_REF_X;
1112   float m_2_1 = space->space.XYZtoRGBf[7] * D50_WHITE_REF_Y;
1113   float m_2_2 = space->space.XYZtoRGBf[8] * D50_WHITE_REF_Z;
1114   long n = samples;
1115 
1116   while (n--)
1117     {
1118       float X, Z, r, g, b;
1119       float x = src[0];
1120       float y = src[1];
1121       float Y = src[2];
1122 
1123 
1124       if (near_zero (Y))
1125         {
1126           X = 0.0f;
1127           Y = 0.0f;
1128           Z = 0.0f;
1129         }
1130       else
1131         {
1132           X = (x * Y) / y;
1133           //Y = Y;
1134           Z = ((1 - x - y) * Y) / y;
1135         }
1136 
1137       r = m_0_0 * X + m_0_1 * Y + m_0_2 * Z;
1138       g = m_1_0 * X + m_1_1 * Y + m_1_2 * Z;
1139       b = m_2_0 * X + m_2_1 * Y + m_2_2 * Z;
1140 
1141       dst[0] =    r;
1142       dst[1] =    g;
1143       dst[2] =    b;
1144       dst[3] = 1.0f;
1145 
1146       src += 3;
1147       dst += 4;
1148     }
1149 }
1150 
1151 static void
xyYaf_to_rgbaf(const Babl * conversion,float * src,float * dst,long samples)1152 xyYaf_to_rgbaf (const Babl *conversion,
1153                 float      *src,
1154                 float      *dst,
1155                 long        samples)
1156 {
1157   const Babl *space = babl_conversion_get_source_space (conversion);
1158   float m_0_0 = space->space.XYZtoRGBf[0] * D50_WHITE_REF_X;
1159   float m_0_1 = space->space.XYZtoRGBf[1] * D50_WHITE_REF_Y;
1160   float m_0_2 = space->space.XYZtoRGBf[2] * D50_WHITE_REF_Z;
1161   float m_1_0 = space->space.XYZtoRGBf[3] * D50_WHITE_REF_X;
1162   float m_1_1 = space->space.XYZtoRGBf[4] * D50_WHITE_REF_Y;
1163   float m_1_2 = space->space.XYZtoRGBf[5] * D50_WHITE_REF_Z;
1164   float m_2_0 = space->space.XYZtoRGBf[6] * D50_WHITE_REF_X;
1165   float m_2_1 = space->space.XYZtoRGBf[7] * D50_WHITE_REF_Y;
1166   float m_2_2 = space->space.XYZtoRGBf[8] * D50_WHITE_REF_Z;
1167   long n = samples;
1168 
1169   while (n--)
1170     {
1171       float X, Z, r, g, b;
1172       float x = src[0];
1173       float y = src[1];
1174       float Y = src[2];
1175       float a = src[3];
1176 
1177       if (near_zero (Y))
1178         {
1179           X = 0.0f;
1180           Y = 0.0f;
1181           Z = 0.0f;
1182         }
1183       else
1184         {
1185           X = (x * Y) / y;
1186           //Y = Y;
1187           Z = ((1 - x - y) * Y) / y;
1188         }
1189 
1190       r = m_0_0 * X + m_0_1 * Y + m_0_2 * Z;
1191       g = m_1_0 * X + m_1_1 * Y + m_1_2 * Z;
1192       b = m_2_0 * X + m_2_1 * Y + m_2_2 * Z;
1193 
1194       dst[0] = r;
1195       dst[1] = g;
1196       dst[2] = b;
1197       dst[3] = a;
1198 
1199       src += 4;
1200       dst += 4;
1201     }
1202 }
1203 
1204 
1205 
1206 /* Yuvf -> rgbf */
1207 
1208 static void
Yuvf_to_rgbf(const Babl * conversion,float * src,float * dst,long samples)1209 Yuvf_to_rgbf (const Babl *conversion,float *src,
1210               float *dst,
1211               long   samples)
1212 {
1213   const Babl *space = babl_conversion_get_source_space (conversion);
1214   float m_0_0 = space->space.XYZtoRGBf[0] * D50_WHITE_REF_X;
1215   float m_0_1 = space->space.XYZtoRGBf[1] * D50_WHITE_REF_Y;
1216   float m_0_2 = space->space.XYZtoRGBf[2] * D50_WHITE_REF_Z;
1217   float m_1_0 = space->space.XYZtoRGBf[3] * D50_WHITE_REF_X;
1218   float m_1_1 = space->space.XYZtoRGBf[4] * D50_WHITE_REF_Y;
1219   float m_1_2 = space->space.XYZtoRGBf[5] * D50_WHITE_REF_Z;
1220   float m_2_0 = space->space.XYZtoRGBf[6] * D50_WHITE_REF_X;
1221   float m_2_1 = space->space.XYZtoRGBf[7] * D50_WHITE_REF_Y;
1222   float m_2_2 = space->space.XYZtoRGBf[8] * D50_WHITE_REF_Z;
1223   long n = samples;
1224 
1225   while (n--)
1226     {
1227       float X, Z, r, g, b;
1228       float Y = src[0];
1229       float u = src[1];
1230       float v = src[2];
1231 
1232       if (near_zero (v))
1233         {
1234           X = 0.0f;
1235           Y = 0.0f;
1236           Z = 0.0f;
1237         }
1238       else
1239         {
1240           X = ((9.0 * u * Y)/(4.0 * v));
1241           //Y = Y;
1242           Z = -(((20.0 * v + 3.0 * u - 12.0) * Y)/(4.0 * v));
1243         }
1244 
1245       r = m_0_0 * X + m_0_1 * Y + m_0_2 * Z;
1246       g = m_1_0 * X + m_1_1 * Y + m_1_2 * Z;
1247       b = m_2_0 * X + m_2_1 * Y + m_2_2 * Z;
1248 
1249       dst[0] = r;
1250       dst[1] = g;
1251       dst[2] = b;
1252 
1253       src += 3;
1254       dst += 3;
1255     }
1256 }
1257 
1258 
1259 
1260 static void
Yuvf_to_rgbaf(const Babl * conversion,float * src,float * dst,long samples)1261 Yuvf_to_rgbaf (const Babl *conversion,
1262                float      *src,
1263                float      *dst,
1264                long        samples)
1265 {
1266   const Babl *space = babl_conversion_get_source_space (conversion);
1267   float m_0_0 = space->space.XYZtoRGBf[0] * D50_WHITE_REF_X;
1268   float m_0_1 = space->space.XYZtoRGBf[1] * D50_WHITE_REF_Y;
1269   float m_0_2 = space->space.XYZtoRGBf[2] * D50_WHITE_REF_Z;
1270   float m_1_0 = space->space.XYZtoRGBf[3] * D50_WHITE_REF_X;
1271   float m_1_1 = space->space.XYZtoRGBf[4] * D50_WHITE_REF_Y;
1272   float m_1_2 = space->space.XYZtoRGBf[5] * D50_WHITE_REF_Z;
1273   float m_2_0 = space->space.XYZtoRGBf[6] * D50_WHITE_REF_X;
1274   float m_2_1 = space->space.XYZtoRGBf[7] * D50_WHITE_REF_Y;
1275   float m_2_2 = space->space.XYZtoRGBf[8] * D50_WHITE_REF_Z;
1276   long n = samples;
1277 
1278   while (n--)
1279     {
1280       float X, Z, r, g, b;
1281       float Y = src[0];
1282       float u = src[1];
1283       float v = src[2];
1284 
1285       if (near_zero (v))
1286         {
1287           X = 0.0f;
1288           Y = 0.0f;
1289           Z = 0.0f;
1290         }
1291       else
1292         {
1293           X = ((9.0 * u * Y)/(4.0 * v));
1294           //Y = Y;
1295           Z = -(((20.0 * v + 3.0 * u - 12.0) * Y)/(4.0 * v));
1296         }
1297 
1298       r = m_0_0 * X + m_0_1 * Y + m_0_2 * Z;
1299       g = m_1_0 * X + m_1_1 * Y + m_1_2 * Z;
1300       b = m_2_0 * X + m_2_1 * Y + m_2_2 * Z;
1301 
1302       dst[0] =    r;
1303       dst[1] =    g;
1304       dst[2] =    b;
1305       dst[3] = 1.0f;
1306 
1307       src += 3;
1308       dst += 4;
1309     }
1310 }
1311 
1312 static void
Yuvaf_to_rgbaf(const Babl * conversion,float * src,float * dst,long samples)1313 Yuvaf_to_rgbaf (const Babl *conversion,
1314                 float      *src,
1315                 float      *dst,
1316                 long        samples)
1317 {
1318   const Babl *space = babl_conversion_get_source_space (conversion);
1319   float m_0_0 = space->space.XYZtoRGBf[0] * D50_WHITE_REF_X;
1320   float m_0_1 = space->space.XYZtoRGBf[1] * D50_WHITE_REF_Y;
1321   float m_0_2 = space->space.XYZtoRGBf[2] * D50_WHITE_REF_Z;
1322   float m_1_0 = space->space.XYZtoRGBf[3] * D50_WHITE_REF_X;
1323   float m_1_1 = space->space.XYZtoRGBf[4] * D50_WHITE_REF_Y;
1324   float m_1_2 = space->space.XYZtoRGBf[5] * D50_WHITE_REF_Z;
1325   float m_2_0 = space->space.XYZtoRGBf[6] * D50_WHITE_REF_X;
1326   float m_2_1 = space->space.XYZtoRGBf[7] * D50_WHITE_REF_Y;
1327   float m_2_2 = space->space.XYZtoRGBf[8] * D50_WHITE_REF_Z;
1328   long n = samples;
1329 
1330   while (n--)
1331     {
1332       float X, Z, r, g, b;
1333       float Y = src[0];
1334       float u = src[1];
1335       float v = src[2];
1336       float a = src[3];
1337 
1338       if (near_zero (v))
1339         {
1340           X = 0.0f;
1341           Y = 0.0f;
1342           Z = 0.0f;
1343         }
1344       else
1345         {
1346           X = ((9.0 * u * Y)/(4.0 * v));
1347           //Y = Y;
1348           Z = -(((20.0 * v + 3.0 * u - 12.0) * Y)/(4.0 * v));
1349         }
1350 
1351       r = m_0_0 * X + m_0_1 * Y + m_0_2 * Z;
1352       g = m_1_0 * X + m_1_1 * Y + m_1_2 * Z;
1353       b = m_2_0 * X + m_2_1 * Y + m_2_2 * Z;
1354 
1355       dst[0] = r;
1356       dst[1] = g;
1357       dst[2] = b;
1358       dst[3] = a;
1359 
1360       src += 4;
1361       dst += 4;
1362     }
1363 }
1364 
1365 
1366 /* rgb <-> LAB */
1367 
1368 static void
rgba_to_lab(const Babl * conversion,char * src,char * dst,long n)1369 rgba_to_lab (const Babl *conversion,
1370              char       *src,
1371              char       *dst,
1372              long        n)
1373 {
1374   const Babl *space = babl_conversion_get_source_space (conversion);
1375   while (n--)
1376     {
1377       double XYZ[3], L, a, b;
1378 
1379       babl_space_to_xyz (space, (double*)src, XYZ);
1380       XYZ_to_LAB (XYZ[0], XYZ[1], XYZ[2], &L, &a, &b);
1381 
1382       ((double *) dst)[0] = L;
1383       ((double *) dst)[1] = a;
1384       ((double *) dst)[2] = b;
1385 
1386       src += sizeof (double) * 4;
1387       dst += sizeof (double) * 3;
1388     }
1389 }
1390 
1391 
1392 static void
lab_to_rgba(const Babl * conversion,char * src,char * dst,long n)1393 lab_to_rgba (const Babl *conversion,
1394              char       *src,
1395              char       *dst,
1396              long        n)
1397 {
1398   const Babl *space = babl_conversion_get_destination_space (conversion);
1399   while (n--)
1400     {
1401       double L = ((double *) src)[0];
1402       double a = ((double *) src)[1];
1403       double b = ((double *) src)[2];
1404 
1405       double X, Y, Z, R, G, B;
1406 
1407       //convert Lab to XYZ
1408       LAB_to_XYZ (L, a, b, &X, &Y, &Z);
1409 
1410       //convert XYZ to RGB
1411       {
1412         double XYZ[3]  = {X,Y,Z};
1413         double RGB[3];
1414         babl_space_from_xyz (space, XYZ, RGB);
1415         R = RGB[0];
1416         G = RGB[1];
1417         B = RGB[2];
1418       }
1419 
1420       ((double *) dst)[0] = R;
1421       ((double *) dst)[1] = G;
1422       ((double *) dst)[2] = B;
1423       ((double *) dst)[3] = 1.0;
1424 
1425       src += sizeof (double) * 3;
1426       dst += sizeof (double) * 4;
1427     }
1428 }
1429 
1430 static void
rgba_to_laba(const Babl * conversion,char * src,char * dst,long n)1431 rgba_to_laba (const Babl *conversion,
1432               char       *src,
1433               char       *dst,
1434               long        n)
1435 {
1436   const Babl *space = babl_conversion_get_source_space (conversion);
1437   while (n--)
1438     {
1439       double alpha = ((double *) src)[3];
1440       double XYZ[3], L, a, b;
1441 
1442       //convert RGB to XYZ
1443       babl_space_to_xyz (space, (double*)src, XYZ);
1444 
1445       //convert XYZ to Lab
1446       XYZ_to_LAB (XYZ[0], XYZ[1], XYZ[2], &L, &a, &b);
1447 
1448       ((double *) dst)[0] = L;
1449       ((double *) dst)[1] = a;
1450       ((double *) dst)[2] = b;
1451       ((double *) dst)[3] = alpha;
1452 
1453       src += sizeof (double) * 4;
1454       dst += sizeof (double) * 4;
1455     }
1456 }
1457 
1458 static void
laba_to_rgba(const Babl * conversion,char * src,char * dst,long n)1459 laba_to_rgba (const Babl *conversion,
1460               char       *src,
1461               char       *dst,
1462               long        n)
1463 {
1464   const Babl *space = babl_conversion_get_destination_space (conversion);
1465   while (n--)
1466     {
1467       double L     = ((double *) src)[0];
1468       double a     = ((double *) src)[1];
1469       double b     = ((double *) src)[2];
1470       double alpha = ((double *) src)[3];
1471 
1472       double X, Y, Z;
1473 
1474       //convert Lab to XYZ
1475       LAB_to_XYZ (L, a, b, &X, &Y, &Z);
1476 
1477       {
1478         //convert XYZ to RGB
1479         double XYZ[3]  = {X,Y,Z};
1480         babl_space_from_xyz (space, XYZ, (double*)dst);
1481       }
1482       ((double *) dst)[3] = alpha;
1483 
1484       src += sizeof (double) * 4;
1485       dst += sizeof (double) * 4;
1486     }
1487 }
1488 
1489 
1490 /* rgb <-> LCh */
1491 
1492 static inline void
CHab_to_ab(double C,double H,double * to_a,double * to_b)1493 CHab_to_ab (double  C,
1494             double  H,
1495             double *to_a,
1496             double *to_b)
1497 {
1498   *to_a = cos (H * RADIANS_PER_DEGREE) * C;
1499   *to_b = sin (H * RADIANS_PER_DEGREE) * C;
1500 }
1501 
1502 static inline void
ab_to_CHab(double a,double b,double * to_C,double * to_H)1503 ab_to_CHab (double  a,
1504             double  b,
1505             double *to_C,
1506             double *to_H)
1507 {
1508   *to_C = sqrt ( (a * a) + (b * b) );
1509   *to_H = atan2 (b, a) * DEGREES_PER_RADIAN;
1510 
1511   // Keep H within the range 0-360
1512   if (*to_H < 0.0)
1513       *to_H += 360;
1514 }
1515 
1516 static void
rgba_to_lchab(const Babl * conversion,char * src,char * dst,long n)1517 rgba_to_lchab (const Babl *conversion,
1518                char       *src,
1519                char       *dst,
1520                long        n)
1521 {
1522   const Babl *space = babl_conversion_get_source_space (conversion);
1523 
1524   while (n--)
1525     {
1526       double XYZ[3], L, a, b, C, H;
1527 
1528       //convert RGB to XYZ
1529       babl_space_to_xyz (space, (double *)src, XYZ);
1530 
1531       //convert XYZ to Lab
1532       XYZ_to_LAB (XYZ[0], XYZ[1], XYZ[2], &L, &a, &b);
1533 
1534 
1535       //convert Lab to LCH(ab)
1536       ab_to_CHab (a, b, &C, &H);
1537 
1538       ((double *) dst)[0] = L;
1539       ((double *) dst)[1] = C;
1540       ((double *) dst)[2] = H;
1541 
1542       src += sizeof (double) * 4;
1543       dst += sizeof (double) * 3;
1544     }
1545 }
1546 
1547 static void
lchab_to_rgba(const Babl * conversion,char * src,char * dst,long n)1548 lchab_to_rgba (const Babl *conversion,
1549                char       *src,
1550                char       *dst,
1551                long        n)
1552 {
1553   const Babl *space = babl_conversion_get_source_space (conversion);
1554 
1555   while (n--)
1556     {
1557       double L = ((double *) src)[0];
1558       double C = ((double *) src)[1];
1559       double H = ((double *) src)[2];
1560       double a, b, X, Y, Z;
1561 
1562       //Convert LCH(ab) to Lab
1563       CHab_to_ab (C, H, &a, &b);
1564 
1565       //Convert LAB to XYZ
1566       LAB_to_XYZ (L, a, b, &X, &Y, &Z);
1567 
1568       //Convert XYZ to RGB
1569       {
1570         double XYZ[3]  = {X,Y,Z};
1571         babl_space_from_xyz (space, XYZ, (double*)dst);
1572       }
1573 
1574       ((double *) dst)[3] = 1.0;
1575 
1576       src += sizeof (double) * 3;
1577       dst += sizeof (double) * 4;
1578     }
1579 }
1580 
1581 static void
rgba_to_lchaba(const Babl * conversion,char * src,char * dst,long n)1582 rgba_to_lchaba (const Babl *conversion,
1583                 char       *src,
1584                 char       *dst,
1585                 long        n)
1586 {
1587   const Babl *space = babl_conversion_get_source_space (conversion);
1588 
1589   while (n--)
1590     {
1591       double alpha = ((double *) src)[3];
1592       double XYZ[3], L, a, b, C, H;
1593 
1594       //convert RGB to XYZ
1595       babl_space_to_xyz (space, (double*)src, XYZ);
1596 
1597       //convert XYZ to Lab
1598       XYZ_to_LAB (XYZ[0], XYZ[1], XYZ[2], &L, &a, &b);
1599 
1600       //convert Lab to LCH(ab)
1601       ab_to_CHab (a, b, &C, &H);
1602 
1603       ((double *) dst)[0] = L;
1604       ((double *) dst)[1] = C;
1605       ((double *) dst)[2] = H;
1606       ((double *) dst)[3] = alpha;
1607 
1608       src += sizeof (double) * 4;
1609       dst += sizeof (double) * 4;
1610     }
1611 }
1612 
1613 static void
lchaba_to_rgba(const Babl * conversion,char * src,char * dst,long n)1614 lchaba_to_rgba (const Babl *conversion,
1615                 char       *src,
1616                 char       *dst,
1617                 long        n)
1618 {
1619   const Babl *space = babl_conversion_get_destination_space (conversion);
1620   while (n--)
1621     {
1622       double L     = ((double *) src)[0];
1623       double C     = ((double *) src)[1];
1624       double H     = ((double *) src)[2];
1625       double alpha = ((double *) src)[3];
1626       double a, b, X, Y, Z;
1627 
1628       //Convert LCH(ab) to Lab
1629       CHab_to_ab (C, H, &a, &b);
1630 
1631       //Convert Lab to XYZ
1632       LAB_to_XYZ (L, a, b, &X, &Y, &Z);
1633 
1634       //Convert XYZ to RGB
1635       {
1636         double XYZ[3]  = {X,Y,Z};
1637         babl_space_from_xyz (space, XYZ, (double*)dst);
1638       }
1639       ((double *) dst)[3] = alpha;
1640 
1641       src += sizeof (double) * 4;
1642       dst += sizeof (double) * 4;
1643     }
1644 }
1645 
1646 
1647 /******** end double RGB/CIE color space conversions ******************/
1648 
1649 /******** begin floating point RGB/CIE color space conversions ********/
1650 
1651 /* origin: http://www.hackersdelight.org/hdcodetxt/acbrt.c.txt
1652  * permissions: http://www.hackersdelight.org/permissions.htm
1653  */
1654 /* _cbrtf(x)
1655  * Return cube root of x
1656  */
1657 
1658 static inline float
_cbrtf(float x)1659 _cbrtf (float x)
1660 {
1661   union { float f; uint32_t i; } u = { x };
1662 
1663   u.i = u.i / 4 + u.i / 16;
1664   u.i = u.i + u.i / 16;
1665   u.i = u.i + u.i / 256;
1666   u.i = 0x2a5137a0 + u.i;
1667   u.f = 0.33333333f * (2.0f * u.f + x / (u.f * u.f));
1668   u.f = 0.33333333f * (2.0f * u.f + x / (u.f * u.f));
1669 
1670   return u.f;
1671 }
1672 
1673 static inline float
cubef(float f)1674 cubef (float f)
1675 {
1676   return f * f * f;
1677 }
1678 
1679 static void
Yf_to_Lf(const Babl * conversion,float * src,float * dst,long samples)1680 Yf_to_Lf (const Babl *conversion,
1681           float      *src,
1682           float      *dst,
1683           long        samples)
1684 {
1685   long n = samples;
1686 
1687   while (n--)
1688     {
1689       float yr = src[0];
1690       float L  = yr > LAB_EPSILON ? 116.0f * _cbrtf (yr) - 16 : LAB_KAPPA * yr;
1691 
1692       dst[0] = L;
1693 
1694       src++;
1695       dst++;
1696     }
1697 }
1698 
1699 static void
Yaf_to_Lf(const Babl * conversion,float * src,float * dst,long samples)1700 Yaf_to_Lf (const Babl *conversion,
1701            float      *src,
1702            float      *dst,
1703            long        samples)
1704 {
1705   long n = samples;
1706 
1707   while (n--)
1708     {
1709       float yr = src[0];
1710       float L  = yr > LAB_EPSILON ? 116.0f * _cbrtf (yr) - 16 : LAB_KAPPA * yr;
1711 
1712       dst[0] = L;
1713 
1714       src += 2;
1715       dst += 1;
1716     }
1717 }
1718 
1719 static void
Yaf_to_Laf(const Babl * conversion,float * src,float * dst,long samples)1720 Yaf_to_Laf (const Babl *conversion,
1721             float      *src,
1722             float      *dst,
1723             long        samples)
1724 {
1725   long n = samples;
1726 
1727   while (n--)
1728     {
1729       float yr = src[0];
1730       float a  = src[1];
1731       float L  = yr > LAB_EPSILON ? 116.0f * _cbrtf (yr) - 16 : LAB_KAPPA * yr;
1732 
1733       dst[0] = L;
1734       dst[1] = a;
1735 
1736       src += 2;
1737       dst += 2;
1738     }
1739 }
1740 
1741 static void
rgbf_to_Labf(const Babl * conversion,float * src,float * dst,long samples)1742 rgbf_to_Labf (const Babl *conversion,
1743               float      *src,
1744               float      *dst,
1745               long        samples)
1746 {
1747   const Babl *space = babl_conversion_get_source_space (conversion);
1748   float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
1749   float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
1750   float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
1751   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
1752   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
1753   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
1754   float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
1755   float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
1756   float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
1757   long n = samples;
1758 
1759   while (n--)
1760     {
1761       float r = src[0];
1762       float g = src[1];
1763       float b = src[2];
1764 
1765       float xr = m_0_0 * r + m_0_1 * g + m_0_2 * b;
1766       float yr = m_1_0 * r + m_1_1 * g + m_1_2 * b;
1767       float zr = m_2_0 * r + m_2_1 * g + m_2_2 * b;
1768 
1769       float fx = xr > LAB_EPSILON ? _cbrtf (xr) : (LAB_KAPPA * xr + 16.0f) / 116.0f;
1770       float fy = yr > LAB_EPSILON ? _cbrtf (yr) : (LAB_KAPPA * yr + 16.0f) / 116.0f;
1771       float fz = zr > LAB_EPSILON ? _cbrtf (zr) : (LAB_KAPPA * zr + 16.0f) / 116.0f;
1772 
1773       float L = 116.0f * fy - 16.0f;
1774       float A = 500.0f * (fx - fy);
1775       float B = 200.0f * (fy - fz);
1776 
1777       dst[0] = L;
1778       dst[1] = A;
1779       dst[2] = B;
1780 
1781       src += 3;
1782       dst += 3;
1783     }
1784 }
1785 
1786 static void
rgbaf_to_Lf(const Babl * conversion,float * src,float * dst,long samples)1787 rgbaf_to_Lf (const Babl *conversion,
1788              float      *src,
1789              float      *dst,
1790              long        samples)
1791 {
1792   const Babl *space = babl_conversion_get_source_space (conversion);
1793   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
1794   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
1795   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
1796   long n = samples;
1797 
1798   while (n--)
1799     {
1800       float r = src[0];
1801       float g = src[1];
1802       float b = src[2];
1803 
1804       float yr = m_1_0 * r + m_1_1 * g + m_1_2 * b;
1805       float L = yr > LAB_EPSILON ? 116.0f * _cbrtf (yr) - 16 : LAB_KAPPA * yr;
1806 
1807       dst[0] = L;
1808 
1809       src += 4;
1810       dst += 1;
1811     }
1812 }
1813 
1814 static void
rgbaf_to_Labf(const Babl * conversion,float * src,float * dst,long samples)1815 rgbaf_to_Labf (const Babl *conversion,
1816                float      *src,
1817                float      *dst,
1818                long        samples)
1819 {
1820   const Babl *space = babl_conversion_get_source_space (conversion);
1821   float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
1822   float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
1823   float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
1824   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
1825   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
1826   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
1827   float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
1828   float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
1829   float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
1830   long n = samples;
1831 
1832   while (n--)
1833     {
1834       float r = src[0];
1835       float g = src[1];
1836       float b = src[2];
1837 
1838       float xr = m_0_0 * r + m_0_1 * g + m_0_2 * b;
1839       float yr = m_1_0 * r + m_1_1 * g + m_1_2 * b;
1840       float zr = m_2_0 * r + m_2_1 * g + m_2_2 * b;
1841 
1842       float fx = xr > LAB_EPSILON ? _cbrtf (xr) : (LAB_KAPPA * xr + 16.0f) / 116.0f;
1843       float fy = yr > LAB_EPSILON ? _cbrtf (yr) : (LAB_KAPPA * yr + 16.0f) / 116.0f;
1844       float fz = zr > LAB_EPSILON ? _cbrtf (zr) : (LAB_KAPPA * zr + 16.0f) / 116.0f;
1845 
1846       float L = 116.0f * fy - 16.0f;
1847       float A = 500.0f * (fx - fy);
1848       float B = 200.0f * (fy - fz);
1849 
1850       dst[0] = L;
1851       dst[1] = A;
1852       dst[2] = B;
1853 
1854       src += 4;
1855       dst += 3;
1856     }
1857 }
1858 
1859 static void
rgbaf_to_Labaf(const Babl * conversion,float * src,float * dst,long samples)1860 rgbaf_to_Labaf (const Babl *conversion,
1861                 float      *src,
1862                 float      *dst,
1863                 long        samples)
1864 {
1865   const Babl *space = babl_conversion_get_source_space (conversion);
1866   float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
1867   float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
1868   float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
1869   float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
1870   float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
1871   float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
1872   float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
1873   float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
1874   float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
1875   long n = samples;
1876 
1877   while (n--)
1878     {
1879       float r = src[0];
1880       float g = src[1];
1881       float b = src[2];
1882       float a = src[3];
1883 
1884       float xr = m_0_0 * r + m_0_1 * g + m_0_2 * b;
1885       float yr = m_1_0 * r + m_1_1 * g + m_1_2 * b;
1886       float zr = m_2_0 * r + m_2_1 * g + m_2_2 * b;
1887 
1888       float fx = xr > LAB_EPSILON ? _cbrtf (xr) : (LAB_KAPPA * xr + 16.0f) / 116.0f;
1889       float fy = yr > LAB_EPSILON ? _cbrtf (yr) : (LAB_KAPPA * yr + 16.0f) / 116.0f;
1890       float fz = zr > LAB_EPSILON ? _cbrtf (zr) : (LAB_KAPPA * zr + 16.0f) / 116.0f;
1891 
1892       float L = 116.0f * fy - 16.0f;
1893       float A = 500.0f * (fx - fy);
1894       float B = 200.0f * (fy - fz);
1895 
1896       dst[0] = L;
1897       dst[1] = A;
1898       dst[2] = B;
1899       dst[3] = a;
1900 
1901       src += 4;
1902       dst += 4;
1903     }
1904 }
1905 
1906 static void
Labf_to_Lf(const Babl * conversion,float * src,float * dst,long samples)1907 Labf_to_Lf (const Babl *conversion,
1908             float      *src,
1909             float      *dst,
1910             long        samples)
1911 {
1912   long n = samples;
1913 
1914   while (n--)
1915     {
1916       dst[0] = src[0];
1917 
1918       src += 3;
1919       dst += 1;
1920     }
1921 }
1922 
1923 static void
Labaf_to_Lf(const Babl * conversion,float * src,float * dst,long samples)1924 Labaf_to_Lf (const Babl *conversion,
1925              float      *src,
1926              float      *dst,
1927              long        samples)
1928 {
1929   long n = samples;
1930 
1931   while (n--)
1932     {
1933       dst[0] = src[0];
1934 
1935       src += 4;
1936       dst += 1;
1937     }
1938 }
1939 
1940 static void
Labf_to_rgbf(const Babl * conversion,float * src,float * dst,long samples)1941 Labf_to_rgbf (const Babl *conversion,
1942               float      *src,
1943               float      *dst,
1944               long        samples)
1945 {
1946   const Babl *space = babl_conversion_get_source_space (conversion);
1947   float m_0_0 = space->space.XYZtoRGBf[0] * D50_WHITE_REF_X;
1948   float m_0_1 = space->space.XYZtoRGBf[1] * D50_WHITE_REF_Y;
1949   float m_0_2 = space->space.XYZtoRGBf[2] * D50_WHITE_REF_Z;
1950   float m_1_0 = space->space.XYZtoRGBf[3] * D50_WHITE_REF_X;
1951   float m_1_1 = space->space.XYZtoRGBf[4] * D50_WHITE_REF_Y;
1952   float m_1_2 = space->space.XYZtoRGBf[5] * D50_WHITE_REF_Z;
1953   float m_2_0 = space->space.XYZtoRGBf[6] * D50_WHITE_REF_X;
1954   float m_2_1 = space->space.XYZtoRGBf[7] * D50_WHITE_REF_Y;
1955   float m_2_2 = space->space.XYZtoRGBf[8] * D50_WHITE_REF_Z;
1956   long n = samples;
1957 
1958   while (n--)
1959     {
1960       float L = src[0];
1961       float A = src[1];
1962       float B = src[2];
1963 
1964       float fy = (L + 16.0f) / 116.0f;
1965       float fx = fy + A / 500.0f;
1966       float fz = fy - B / 200.0f;
1967 
1968       float yr = L > LAB_KAPPA * LAB_EPSILON ? cubef (fy) : L / LAB_KAPPA;
1969       float xr = cubef (fx) > LAB_EPSILON ? cubef (fx) : (fx * 116.0f - 16.0f) / LAB_KAPPA;
1970       float zr = cubef (fz) > LAB_EPSILON ? cubef (fz) : (fz * 116.0f - 16.0f) / LAB_KAPPA;
1971 
1972       float r = m_0_0 * xr + m_0_1 * yr + m_0_2 * zr;
1973       float g = m_1_0 * xr + m_1_1 * yr + m_1_2 * zr;
1974       float b = m_2_0 * xr + m_2_1 * yr + m_2_2 * zr;
1975 
1976       dst[0] = r;
1977       dst[1] = g;
1978       dst[2] = b;
1979 
1980       src += 3;
1981       dst += 3;
1982     }
1983 }
1984 
1985 
1986 static void
Labf_to_rgbaf(const Babl * conversion,float * src,float * dst,long samples)1987 Labf_to_rgbaf (const Babl *conversion,float *src,
1988                float *dst,
1989                long   samples)
1990 {
1991   const Babl *space = babl_conversion_get_source_space (conversion);
1992   float m_0_0 = space->space.XYZtoRGBf[0] * D50_WHITE_REF_X;
1993   float m_0_1 = space->space.XYZtoRGBf[1] * D50_WHITE_REF_Y;
1994   float m_0_2 = space->space.XYZtoRGBf[2] * D50_WHITE_REF_Z;
1995   float m_1_0 = space->space.XYZtoRGBf[3] * D50_WHITE_REF_X;
1996   float m_1_1 = space->space.XYZtoRGBf[4] * D50_WHITE_REF_Y;
1997   float m_1_2 = space->space.XYZtoRGBf[5] * D50_WHITE_REF_Z;
1998   float m_2_0 = space->space.XYZtoRGBf[6] * D50_WHITE_REF_X;
1999   float m_2_1 = space->space.XYZtoRGBf[7] * D50_WHITE_REF_Y;
2000   float m_2_2 = space->space.XYZtoRGBf[8] * D50_WHITE_REF_Z;
2001   long n = samples;
2002 
2003   while (n--)
2004     {
2005       float L = src[0];
2006       float A = src[1];
2007       float B = src[2];
2008 
2009       float fy = (L + 16.0f) / 116.0f;
2010       float fx = fy + A / 500.0f;
2011       float fz = fy - B / 200.0f;
2012 
2013       float yr = L > LAB_KAPPA * LAB_EPSILON ? cubef (fy) : L / LAB_KAPPA;
2014       float xr = cubef (fx) > LAB_EPSILON ? cubef (fx) : (fx * 116.0f - 16.0f) / LAB_KAPPA;
2015       float zr = cubef (fz) > LAB_EPSILON ? cubef (fz) : (fz * 116.0f - 16.0f) / LAB_KAPPA;
2016 
2017       float r = m_0_0 * xr + m_0_1 * yr + m_0_2 * zr;
2018       float g = m_1_0 * xr + m_1_1 * yr + m_1_2 * zr;
2019       float b = m_2_0 * xr + m_2_1 * yr + m_2_2 * zr;
2020 
2021       dst[0] = r;
2022       dst[1] = g;
2023       dst[2] = b;
2024       dst[3] = 1.0f;
2025 
2026       src += 3;
2027       dst += 4;
2028     }
2029 }
2030 
2031 static void
Labaf_to_rgbaf(const Babl * conversion,float * src,float * dst,long samples)2032 Labaf_to_rgbaf (const Babl *conversion,
2033                 float      *src,
2034                 float      *dst,
2035                 long        samples)
2036 {
2037   const Babl *space = babl_conversion_get_source_space (conversion);
2038   float m_0_0 = space->space.XYZtoRGBf[0] * D50_WHITE_REF_X;
2039   float m_0_1 = space->space.XYZtoRGBf[1] * D50_WHITE_REF_Y;
2040   float m_0_2 = space->space.XYZtoRGBf[2] * D50_WHITE_REF_Z;
2041   float m_1_0 = space->space.XYZtoRGBf[3] * D50_WHITE_REF_X;
2042   float m_1_1 = space->space.XYZtoRGBf[4] * D50_WHITE_REF_Y;
2043   float m_1_2 = space->space.XYZtoRGBf[5] * D50_WHITE_REF_Z;
2044   float m_2_0 = space->space.XYZtoRGBf[6] * D50_WHITE_REF_X;
2045   float m_2_1 = space->space.XYZtoRGBf[7] * D50_WHITE_REF_Y;
2046   float m_2_2 = space->space.XYZtoRGBf[8] * D50_WHITE_REF_Z;
2047   long n = samples;
2048 
2049   while (n--)
2050     {
2051       float L = src[0];
2052       float A = src[1];
2053       float B = src[2];
2054       float a = src[3];
2055 
2056       float fy = (L + 16.0f) / 116.0f;
2057       float fx = fy + A / 500.0f;
2058       float fz = fy - B / 200.0f;
2059 
2060       float yr = L > LAB_KAPPA * LAB_EPSILON ? cubef (fy) : L / LAB_KAPPA;
2061       float xr = cubef (fx) > LAB_EPSILON ? cubef (fx) : (fx * 116.0f - 16.0f) / LAB_KAPPA;
2062       float zr = cubef (fz) > LAB_EPSILON ? cubef (fz) : (fz * 116.0f - 16.0f) / LAB_KAPPA;
2063 
2064       float r = m_0_0 * xr + m_0_1 * yr + m_0_2 * zr;
2065       float g = m_1_0 * xr + m_1_1 * yr + m_1_2 * zr;
2066       float b = m_2_0 * xr + m_2_1 * yr + m_2_2 * zr;
2067 
2068       dst[0] = r;
2069       dst[1] = g;
2070       dst[2] = b;
2071       dst[3] = a;
2072 
2073       src += 4;
2074       dst += 4;
2075     }
2076 }
2077 
2078 static void
Labf_to_Lchabf(const Babl * conversion,float * src,float * dst,long samples)2079 Labf_to_Lchabf (const Babl *conversion,
2080                 float      *src,
2081                 float      *dst,
2082                 long        samples)
2083 {
2084   long n = samples;
2085 
2086   while (n--)
2087     {
2088       float L = src[0];
2089       float A = src[1];
2090       float B = src[2];
2091 
2092       float C = sqrtf (A * A + B * B);
2093       float H = atan2f (B, A) * DEGREES_PER_RADIAN;
2094 
2095       // Keep H within the range 0-360
2096       if (H < 0.0f)
2097         H += 360.0f;
2098 
2099       dst[0] = L;
2100       dst[1] = C;
2101       dst[2] = H;
2102 
2103       src += 3;
2104       dst += 3;
2105     }
2106 }
2107 
2108 static void
Lchabf_to_Labf(const Babl * conversion,float * src,float * dst,long samples)2109 Lchabf_to_Labf (const Babl *conversion,
2110                 float      *src,
2111                 float      *dst,
2112                 long        samples)
2113 {
2114   long n = samples;
2115 
2116   while (n--)
2117     {
2118       float L = src[0];
2119       float C = src[1];
2120       float H = src[2];
2121 
2122       float A = C * cosf (H * RADIANS_PER_DEGREE);
2123       float B = C * sinf (H * RADIANS_PER_DEGREE);
2124 
2125       dst[0] = L;
2126       dst[1] = A;
2127       dst[2] = B;
2128 
2129       src += 3;
2130       dst += 3;
2131     }
2132 }
2133 
2134 static void
Labaf_to_Lchabaf(const Babl * conversion,float * src,float * dst,long samples)2135 Labaf_to_Lchabaf (const Babl *conversion,
2136                   float      *src,
2137                   float      *dst,
2138                   long        samples)
2139 {
2140   long n = samples;
2141 
2142   while (n--)
2143     {
2144       float L = src[0];
2145       float A = src[1];
2146       float B = src[2];
2147       float a = src[3];
2148 
2149       float C = sqrtf (A * A + B * B);
2150       float H = atan2f (B, A) * DEGREES_PER_RADIAN;
2151 
2152       // Keep H within the range 0-360
2153       if (H < 0.0f)
2154         H += 360.0f;
2155 
2156       dst[0] = L;
2157       dst[1] = C;
2158       dst[2] = H;
2159       dst[3] = a;
2160 
2161       src += 4;
2162       dst += 4;
2163     }
2164 }
2165 
2166 static void
Lchabaf_to_Labaf(const Babl * conversion,float * src,float * dst,long samples)2167 Lchabaf_to_Labaf (const Babl *conversion,
2168                   float      *src,
2169                   float      *dst,
2170                   long        samples)
2171 {
2172   long n = samples;
2173 
2174   while (n--)
2175     {
2176       float L = src[0];
2177       float C = src[1];
2178       float H = src[2];
2179       float a = src[3];
2180 
2181       float A = C * cosf (H * RADIANS_PER_DEGREE);
2182       float B = C * sinf (H * RADIANS_PER_DEGREE);
2183 
2184       dst[0] = L;
2185       dst[1] = A;
2186       dst[2] = B;
2187       dst[3] = a;
2188 
2189       src += 4;
2190       dst += 4;
2191     }
2192 }
2193 
2194 #if defined(USE_SSE2)
2195 
2196 /* This is an SSE2 version of Halley's method for approximating the
2197  * cube root of an IEEE float implementation.
2198  *
2199  * The scalar version is as follows:
2200  *
2201  * static inline float
2202  * _cbrt_5f (float x)
2203  * {
2204  *   union { float f; uint32_t i; } u = { x };
2205  *
2206  *   u.i = u.i / 3 + 709921077;
2207  *   return u.f;
2208  * }
2209  *
2210  * static inline float
2211  * _cbrta_halleyf (float a, float R)
2212  * {
2213  *   float a3 = a * a * a;
2214  *   float b = a * (a3 + R + R) / (a3 + a3 + R);
2215  *   return b;
2216  * }
2217  *
2218  * static inline float
2219  * _cbrtf (float x)
2220  * {
2221  *   float a;
2222  *
2223  *   a = _cbrt_5f (x);
2224  *   a = _cbrta_halleyf (a, x);
2225  *   a = _cbrta_halleyf (a, x);
2226  *   return a;
2227  * }
2228  *
2229  * The above scalar version seems to have originated from
2230  * http://metamerist.com/cbrt/cbrt.htm but that's not accessible
2231  * anymore. At present there's a copy in CubeRoot.cpp in the Skia
2232  * sources that's licensed under a BSD-style license. There's some
2233  * discussion on the implementation at
2234  * http://www.voidcn.com/article/p-gpwztojr-wt.html.
2235  *
2236  * Note that Darktable also has an SSE2 version of the same algorithm,
2237  * but uses only a single iteration of Halley's method, which is too
2238  * coarse.
2239  */
2240 /* Return cube roots of the four single-precision floating point
2241  * components of x.
2242  */
2243 static inline __m128
_cbrtf_ps_sse2(__m128 x)2244 _cbrtf_ps_sse2 (__m128 x)
2245 {
2246   const __m128i magic = _mm_set1_epi32 (709921077);
2247 
2248   __m128i xi = _mm_castps_si128 (x);
2249   __m128 xi_3 = _mm_div_ps (_mm_cvtepi32_ps (xi), _mm_set1_ps (3.0f));
2250   __m128i ai = _mm_add_epi32 (_mm_cvtps_epi32 (xi_3), magic);
2251   __m128 a = _mm_castsi128_ps (ai);
2252 
2253   __m128 a3 = _mm_mul_ps (_mm_mul_ps (a, a), a);
2254   __m128 divisor = _mm_add_ps (_mm_add_ps (a3, a3), x);
2255   a = _mm_div_ps (_mm_mul_ps (a, _mm_add_ps (a3, _mm_add_ps (x, x))), divisor);
2256 
2257   a3 = _mm_mul_ps (_mm_mul_ps (a, a), a);
2258   divisor = _mm_add_ps (_mm_add_ps (a3, a3), x);
2259   a = _mm_div_ps (_mm_mul_ps (a, _mm_add_ps (a3, _mm_add_ps (x, x))), divisor);
2260 
2261   return a;
2262 }
2263 
2264 static inline __m128
lab_r_to_f_sse2(__m128 r)2265 lab_r_to_f_sse2 (__m128 r)
2266 {
2267   const __m128 epsilon = _mm_set1_ps (LAB_EPSILON);
2268   const __m128 kappa = _mm_set1_ps (LAB_KAPPA);
2269 
2270   const __m128 f_big = _cbrtf_ps_sse2 (r);
2271 
2272   const __m128 f_small = _mm_div_ps (_mm_add_ps (_mm_mul_ps (kappa, r), _mm_set1_ps (16.0f)),
2273                                      _mm_set1_ps (116.0f));
2274 
2275   const __m128 mask = _mm_cmpgt_ps (r, epsilon);
2276   const __m128 f = _mm_or_ps (_mm_and_ps (mask, f_big), _mm_andnot_ps (mask, f_small));
2277   return f;
2278 }
2279 
2280 static void
Yf_to_Lf_sse2(const Babl * conversion,const float * src,float * dst,long samples)2281 Yf_to_Lf_sse2 (const Babl  *conversion,
2282                const float *src,
2283                float       *dst,
2284                long         samples)
2285 {
2286   long i = 0;
2287   long remainder;
2288 
2289   if (((uintptr_t) src % 16) + ((uintptr_t) dst % 16) == 0)
2290     {
2291       const long n = (samples / 4) * 4;
2292 
2293       for ( ; i < n; i += 4)
2294         {
2295           __m128 Y = _mm_load_ps (src);
2296 
2297           __m128 fy = lab_r_to_f_sse2 (Y);
2298 
2299           __m128 L = _mm_sub_ps (_mm_mul_ps (_mm_set1_ps (116.0f), fy), _mm_set1_ps (16.0f));
2300 
2301           _mm_store_ps (dst, L);
2302 
2303           src += 4;
2304           dst += 4;
2305         }
2306     }
2307 
2308   remainder = samples - i;
2309   while (remainder--)
2310     {
2311       float yr = src[0];
2312       float L  = yr > LAB_EPSILON ? 116.0f * _cbrtf (yr) - 16 : LAB_KAPPA * yr;
2313 
2314       dst[0] = L;
2315 
2316       src++;
2317       dst++;
2318     }
2319 }
2320 
2321 static void
Yaf_to_Lf_sse2(const Babl * conversion,const float * src,float * dst,long samples)2322 Yaf_to_Lf_sse2 (const Babl  *conversion,
2323                 const float *src,
2324                 float       *dst,
2325                 long         samples)
2326 {
2327   long i = 0;
2328   long remainder;
2329 
2330   if (((uintptr_t) src % 16) + ((uintptr_t) dst % 16) == 0)
2331     {
2332       const long n = (samples / 4) * 4;
2333 
2334       for ( ; i < n; i += 4)
2335         {
2336           __m128 YaYa0 = _mm_load_ps (src);
2337           __m128 YaYa1 = _mm_load_ps (src + 4);
2338 
2339           __m128 Y = _mm_shuffle_ps (YaYa0, YaYa1, _MM_SHUFFLE (2, 0, 2, 0));
2340 
2341           __m128 fy = lab_r_to_f_sse2 (Y);
2342 
2343           __m128 L = _mm_sub_ps (_mm_mul_ps (_mm_set1_ps (116.0f), fy), _mm_set1_ps (16.0f));
2344 
2345           _mm_store_ps (dst, L);
2346 
2347           src += 8;
2348           dst += 4;
2349         }
2350     }
2351 
2352   remainder = samples - i;
2353   while (remainder--)
2354     {
2355       float yr = src[0];
2356       float L  = yr > LAB_EPSILON ? 116.0f * _cbrtf (yr) - 16 : LAB_KAPPA * yr;
2357 
2358       dst[0] = L;
2359 
2360       src += 2;
2361       dst += 1;
2362     }
2363 }
2364 
2365 static void
rgbaf_to_Lf_sse2(const Babl * conversion,const float * src,float * dst,long samples)2366 rgbaf_to_Lf_sse2 (const Babl  *conversion,
2367                   const float *src,
2368                   float       *dst,
2369                   long         samples)
2370 {
2371   const Babl *space = babl_conversion_get_source_space (conversion);
2372   const float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
2373   const float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
2374   const float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
2375   long i = 0;
2376   long remainder;
2377 
2378   if (((uintptr_t) src % 16) + ((uintptr_t) dst % 16) == 0)
2379     {
2380       const long    n = (samples / 4) * 4;
2381       const __m128 m_1_0_v = _mm_set1_ps (m_1_0);
2382       const __m128 m_1_1_v = _mm_set1_ps (m_1_1);
2383       const __m128 m_1_2_v = _mm_set1_ps (m_1_2);
2384 
2385       for ( ; i < n; i += 4)
2386         {
2387           __m128 rgba0 = _mm_load_ps (src);
2388           __m128 rgba1 = _mm_load_ps (src + 4);
2389           __m128 rgba2 = _mm_load_ps (src + 8);
2390           __m128 rgba3 = _mm_load_ps (src + 12);
2391 
2392           __m128 r = rgba0;
2393           __m128 g = rgba1;
2394           __m128 b = rgba2;
2395           __m128 a = rgba3;
2396           _MM_TRANSPOSE4_PS (r, g, b, a);
2397 
2398           {
2399             __m128 yr = _mm_add_ps (_mm_add_ps (_mm_mul_ps (m_1_0_v, r), _mm_mul_ps (m_1_1_v, g)),
2400                                     _mm_mul_ps (m_1_2_v, b));
2401 
2402             __m128 fy = lab_r_to_f_sse2 (yr);
2403 
2404             __m128 L = _mm_sub_ps (_mm_mul_ps (_mm_set1_ps (116.0f), fy), _mm_set1_ps (16.0f));
2405 
2406             _mm_store_ps (dst, L);
2407           }
2408 
2409           src += 16;
2410           dst += 4;
2411         }
2412     }
2413 
2414   remainder = samples - i;
2415   while (remainder--)
2416     {
2417       float r = src[0];
2418       float g = src[1];
2419       float b = src[2];
2420 
2421       float yr = m_1_0 * r + m_1_1 * g + m_1_2 * b;
2422       float L = yr > LAB_EPSILON ? 116.0f * _cbrtf (yr) - 16 : LAB_KAPPA * yr;
2423 
2424       dst[0] = L;
2425 
2426       src += 4;
2427       dst += 1;
2428     }
2429 }
2430 
2431 static void
rgbaf_to_Labaf_sse2(const Babl * conversion,const float * src,float * dst,long samples)2432 rgbaf_to_Labaf_sse2 (const Babl  *conversion,
2433                      const float *src,
2434                      float       *dst,
2435                      long         samples)
2436 {
2437   const Babl *space = babl_conversion_get_source_space (conversion);
2438   const float m_0_0 = space->space.RGBtoXYZf[0] / D50_WHITE_REF_X;
2439   const float m_0_1 = space->space.RGBtoXYZf[1] / D50_WHITE_REF_X;
2440   const float m_0_2 = space->space.RGBtoXYZf[2] / D50_WHITE_REF_X;
2441   const float m_1_0 = space->space.RGBtoXYZf[3] / D50_WHITE_REF_Y;
2442   const float m_1_1 = space->space.RGBtoXYZf[4] / D50_WHITE_REF_Y;
2443   const float m_1_2 = space->space.RGBtoXYZf[5] / D50_WHITE_REF_Y;
2444   const float m_2_0 = space->space.RGBtoXYZf[6] / D50_WHITE_REF_Z;
2445   const float m_2_1 = space->space.RGBtoXYZf[7] / D50_WHITE_REF_Z;
2446   const float m_2_2 = space->space.RGBtoXYZf[8] / D50_WHITE_REF_Z;
2447   long i = 0;
2448   long remainder;
2449 
2450   if (((uintptr_t) src % 16) + ((uintptr_t) dst % 16) == 0)
2451     {
2452       const long    n = (samples / 4) * 4;
2453       const __m128 m_0_0_v = _mm_set1_ps (m_0_0);
2454       const __m128 m_0_1_v = _mm_set1_ps (m_0_1);
2455       const __m128 m_0_2_v = _mm_set1_ps (m_0_2);
2456       const __m128 m_1_0_v = _mm_set1_ps (m_1_0);
2457       const __m128 m_1_1_v = _mm_set1_ps (m_1_1);
2458       const __m128 m_1_2_v = _mm_set1_ps (m_1_2);
2459       const __m128 m_2_0_v = _mm_set1_ps (m_2_0);
2460       const __m128 m_2_1_v = _mm_set1_ps (m_2_1);
2461       const __m128 m_2_2_v = _mm_set1_ps (m_2_2);
2462 
2463       for ( ; i < n; i += 4)
2464         {
2465           __m128 Laba0;
2466           __m128 Laba1;
2467           __m128 Laba2;
2468           __m128 Laba3;
2469 
2470           __m128 rgba0 = _mm_load_ps (src);
2471           __m128 rgba1 = _mm_load_ps (src + 4);
2472           __m128 rgba2 = _mm_load_ps (src + 8);
2473           __m128 rgba3 = _mm_load_ps (src + 12);
2474 
2475           __m128 r = rgba0;
2476           __m128 g = rgba1;
2477           __m128 b = rgba2;
2478           __m128 a = rgba3;
2479           _MM_TRANSPOSE4_PS (r, g, b, a);
2480 
2481           {
2482             __m128 xr = _mm_add_ps (_mm_add_ps (_mm_mul_ps (m_0_0_v, r), _mm_mul_ps (m_0_1_v, g)),
2483                                     _mm_mul_ps (m_0_2_v, b));
2484             __m128 yr = _mm_add_ps (_mm_add_ps (_mm_mul_ps (m_1_0_v, r), _mm_mul_ps (m_1_1_v, g)),
2485                                     _mm_mul_ps (m_1_2_v, b));
2486             __m128 zr = _mm_add_ps (_mm_add_ps (_mm_mul_ps (m_2_0_v, r), _mm_mul_ps (m_2_1_v, g)),
2487                                     _mm_mul_ps (m_2_2_v, b));
2488 
2489             __m128 fx = lab_r_to_f_sse2 (xr);
2490             __m128 fy = lab_r_to_f_sse2 (yr);
2491             __m128 fz = lab_r_to_f_sse2 (zr);
2492 
2493             __m128 L = _mm_sub_ps (_mm_mul_ps (_mm_set1_ps (116.0f), fy), _mm_set1_ps (16.0f));
2494             __m128 A = _mm_mul_ps (_mm_set1_ps (500.0f), _mm_sub_ps (fx, fy));
2495             __m128 B = _mm_mul_ps (_mm_set1_ps (200.0f), _mm_sub_ps (fy, fz));
2496 
2497             Laba0 = L;
2498             Laba1 = A;
2499             Laba2 = B;
2500             Laba3 = a;
2501             _MM_TRANSPOSE4_PS (Laba0, Laba1, Laba2, Laba3);
2502           }
2503 
2504           _mm_store_ps (dst, Laba0);
2505           _mm_store_ps (dst + 4, Laba1);
2506           _mm_store_ps (dst + 8, Laba2);
2507           _mm_store_ps (dst + 12, Laba3);
2508 
2509           src += 16;
2510           dst += 16;
2511         }
2512     }
2513 
2514   remainder = samples - i;
2515   while (remainder--)
2516     {
2517       float r = src[0];
2518       float g = src[1];
2519       float b = src[2];
2520       float a = src[3];
2521 
2522       float xr = m_0_0 * r + m_0_1 * g + m_0_2 * b;
2523       float yr = m_1_0 * r + m_1_1 * g + m_1_2 * b;
2524       float zr = m_2_0 * r + m_2_1 * g + m_2_2 * b;
2525 
2526       float fx = xr > LAB_EPSILON ? _cbrtf (xr) : (LAB_KAPPA * xr + 16.0f) / 116.0f;
2527       float fy = yr > LAB_EPSILON ? _cbrtf (yr) : (LAB_KAPPA * yr + 16.0f) / 116.0f;
2528       float fz = zr > LAB_EPSILON ? _cbrtf (zr) : (LAB_KAPPA * zr + 16.0f) / 116.0f;
2529 
2530       float L = 116.0f * fy - 16.0f;
2531       float A = 500.0f * (fx - fy);
2532       float B = 200.0f * (fy - fz);
2533 
2534       dst[0] = L;
2535       dst[1] = A;
2536       dst[2] = B;
2537       dst[3] = a;
2538 
2539       src += 4;
2540       dst += 4;
2541     }
2542 }
2543 
2544 #endif /* defined(USE_SSE2) */
2545 
2546 static void
conversions(void)2547 conversions (void)
2548 {
2549   /* babl_model */
2550 
2551   babl_conversion_new (
2552     babl_model ("RGBA"),
2553     babl_model ("CIE Lab"),
2554     "linear", rgba_to_lab,
2555     NULL
2556   );
2557   babl_conversion_new (
2558     babl_model ("CIE Lab"),
2559     babl_model ("RGBA"),
2560     "linear", lab_to_rgba,
2561     NULL
2562   );
2563   babl_conversion_new (
2564     babl_model ("RGBA"),
2565     babl_model ("CIE Lab alpha"),
2566     "linear", rgba_to_laba,
2567     NULL
2568   );
2569     babl_conversion_new (
2570     babl_model ("CIE Lab alpha"),
2571     babl_model ("RGBA"),
2572     "linear", laba_to_rgba,
2573     NULL
2574   );
2575   babl_conversion_new (
2576     babl_model ("RGBA"),
2577     babl_model ("CIE LCH(ab)"),
2578     "linear", rgba_to_lchab,
2579     NULL
2580   );
2581   babl_conversion_new (
2582     babl_model ("CIE LCH(ab)"),
2583     babl_model ("RGBA"),
2584     "linear", lchab_to_rgba,
2585     NULL
2586   );
2587   babl_conversion_new (
2588     babl_model ("RGBA"),
2589     babl_model ("CIE LCH(ab) alpha"),
2590     "linear", rgba_to_lchaba,
2591     NULL
2592   );
2593   babl_conversion_new (
2594     babl_model ("CIE LCH(ab) alpha"),
2595     babl_model ("RGBA"),
2596     "linear", lchaba_to_rgba,
2597     NULL
2598   );
2599   babl_conversion_new (
2600     babl_model ("RGBA"),
2601     babl_model ("CIE XYZ"),
2602     "linear", rgba_to_xyz,
2603     NULL
2604   );
2605   babl_conversion_new (
2606     babl_model ("CIE XYZ"),
2607     babl_model ("RGBA"),
2608     "linear", xyz_to_rgba,
2609     NULL
2610   );
2611   babl_conversion_new (
2612     babl_model ("RGBA"),
2613     babl_model ("CIE XYZ alpha"),
2614     "linear", rgba_to_xyza,
2615     NULL
2616   );
2617   babl_conversion_new (
2618     babl_model ("CIE XYZ alpha"),
2619     babl_model ("RGBA"),
2620     "linear", xyza_to_rgba,
2621     NULL
2622   );
2623 
2624   /* CIE xyY */
2625   babl_conversion_new (
2626     babl_model ("RGBA"),
2627     babl_model ("CIE xyY"),
2628     "linear", rgba_to_xyY,
2629     NULL
2630   );
2631   babl_conversion_new (
2632     babl_model ("CIE xyY"),
2633     babl_model ("RGBA"),
2634     "linear", xyY_to_rgba,
2635     NULL
2636   );
2637   babl_conversion_new (
2638     babl_model ("RGBA"),
2639     babl_model ("CIE xyY alpha"),
2640     "linear", rgba_to_xyYa,
2641     NULL
2642   );
2643   babl_conversion_new (
2644     babl_model ("CIE xyY alpha"),
2645     babl_model ("RGBA"),
2646     "linear", xyYa_to_rgba,
2647     NULL
2648   );
2649 
2650   /* CIE 1976 UCS */
2651   babl_conversion_new (
2652     babl_model ("RGBA"),
2653     babl_model ("CIE Yuv"),
2654     "linear", rgba_to_Yuv,
2655     NULL
2656   );
2657   babl_conversion_new (
2658     babl_model ("CIE Yuv"),
2659     babl_model ("RGBA"),
2660     "linear", Yuv_to_rgba,
2661     NULL
2662   );
2663   babl_conversion_new (
2664     babl_model ("RGBA"),
2665     babl_model ("CIE Yuv alpha"),
2666     "linear", rgba_to_Yuva,
2667     NULL
2668   );
2669   babl_conversion_new (
2670     babl_model ("CIE Yuv alpha"),
2671     babl_model ("RGBA"),
2672     "linear", Yuva_to_rgba,
2673     NULL
2674   );
2675 
2676   /* babl_format */
2677 
2678   babl_conversion_new (
2679     babl_format ("RGB float"),
2680     babl_format ("CIE Lab float"),
2681     "linear", rgbf_to_Labf,
2682     NULL
2683   );
2684   babl_conversion_new (
2685     babl_format ("RGBA float"),
2686     babl_format ("CIE Lab float"),
2687     "linear", rgbaf_to_Labf,
2688     NULL
2689   );
2690   babl_conversion_new (
2691     babl_format ("RGBA float"),
2692     babl_format ("CIE Lab alpha float"),
2693     "linear", rgbaf_to_Labaf,
2694     NULL
2695   );
2696   babl_conversion_new (
2697     babl_format ("CIE Lab float"),
2698     babl_format ("RGB float"),
2699     "linear", Labf_to_rgbf,
2700     NULL
2701   );
2702   babl_conversion_new (
2703     babl_format ("CIE Lab float"),
2704     babl_format ("RGBA float"),
2705     "linear", Labf_to_rgbaf,
2706     NULL
2707   );
2708   babl_conversion_new (
2709     babl_format ("CIE Lab alpha float"),
2710     babl_format ("RGBA float"),
2711     "linear", Labaf_to_rgbaf,
2712     NULL
2713   );
2714   babl_conversion_new (
2715     babl_format ("Y float"),
2716     babl_format ("CIE L float"),
2717     "linear", Yf_to_Lf,
2718     NULL
2719   );
2720   babl_conversion_new (
2721     babl_format ("YA float"),
2722     babl_format ("CIE L float"),
2723     "linear", Yaf_to_Lf,
2724     NULL
2725   );
2726   babl_conversion_new (
2727     babl_format ("YA float"),
2728     babl_format ("CIE L alpha float"),
2729     "linear", Yaf_to_Laf,
2730     NULL
2731   );
2732   babl_conversion_new (
2733     babl_format ("RGBA float"),
2734     babl_format ("CIE L float"),
2735     "linear", rgbaf_to_Lf,
2736     NULL
2737   );
2738   babl_conversion_new (
2739     babl_format ("CIE Lab float"),
2740     babl_format ("CIE L float"),
2741     "linear", Labf_to_Lf,
2742     NULL
2743   );
2744   babl_conversion_new (
2745     babl_format ("CIE Lab alpha float"),
2746     babl_format ("CIE L float"),
2747     "linear", Labaf_to_Lf,
2748     NULL
2749   );
2750   babl_conversion_new (
2751     babl_format ("CIE Lab float"),
2752     babl_format ("CIE LCH(ab) float"),
2753     "linear", Labf_to_Lchabf,
2754     NULL
2755   );
2756   babl_conversion_new (
2757     babl_format ("CIE LCH(ab) float"),
2758     babl_format ("CIE Lab float"),
2759     "linear", Lchabf_to_Labf,
2760     NULL
2761   );
2762   babl_conversion_new (
2763     babl_format ("CIE Lab alpha float"),
2764     babl_format ("CIE LCH(ab) alpha float"),
2765     "linear", Labaf_to_Lchabaf,
2766     NULL
2767   );
2768   babl_conversion_new (
2769     babl_format ("CIE LCH(ab) alpha float"),
2770     babl_format ("CIE Lab alpha float"),
2771     "linear", Lchabaf_to_Labaf,
2772     NULL
2773   );
2774 
2775   /* CIE xyY */
2776   babl_conversion_new (
2777     babl_format ("RGB float"),
2778     babl_format ("CIE xyY float"),
2779     "linear", rgbf_to_xyYf,
2780     NULL
2781   );
2782   babl_conversion_new (
2783     babl_format ("RGBA float"),
2784     babl_format ("CIE xyY alpha float"),
2785     "linear", rgbaf_to_xyYaf,
2786     NULL
2787   );
2788   babl_conversion_new (
2789     babl_format ("RGBA float"),
2790     babl_format ("CIE xyY float"),
2791     "linear", rgbaf_to_xyYf,
2792     NULL
2793   );
2794   babl_conversion_new (
2795     babl_format ("CIE xyY float"),
2796     babl_format ("RGB float"),
2797     "linear", xyYf_to_rgbf,
2798     NULL
2799   );
2800   babl_conversion_new (
2801     babl_format ("CIE xyY float"),
2802     babl_format ("RGBA float"),
2803     "linear", xyYf_to_rgbaf,
2804     NULL
2805   );
2806   babl_conversion_new (
2807     babl_format ("CIE xyY alpha float"),
2808     babl_format ("RGBA float"),
2809     "linear", xyYaf_to_rgbaf,
2810     NULL
2811   );
2812   /* CIE 1976 UCS */
2813   babl_conversion_new (
2814     babl_format ("RGB float"),
2815     babl_format ("CIE Yuv float"),
2816     "linear", rgbf_to_Yuvf,
2817     NULL
2818   );
2819   babl_conversion_new (
2820     babl_format ("RGBA float"),
2821     babl_format ("CIE Yuv alpha float"),
2822     "linear", rgbaf_to_Yuvaf,
2823     NULL
2824   );
2825   babl_conversion_new (
2826     babl_format ("RGBA float"),
2827     babl_format ("CIE Yuv float"),
2828     "linear", rgbaf_to_Yuvf,
2829     NULL
2830   );
2831   babl_conversion_new (
2832     babl_format ("CIE Yuv float"),
2833     babl_format ("RGB float"),
2834     "linear", Yuvf_to_rgbf,
2835     NULL
2836   );
2837   babl_conversion_new (
2838     babl_format ("CIE Yuv float"),
2839     babl_format ("RGBA float"),
2840     "linear", Yuvf_to_rgbaf,
2841     NULL
2842   );
2843   babl_conversion_new (
2844     babl_format ("CIE Yuv alpha float"),
2845     babl_format ("RGBA float"),
2846     "linear", Yuvaf_to_rgbaf,
2847     NULL
2848   );
2849 #if defined(USE_SSE2)
2850 
2851   if (babl_cpu_accel_get_support () & BABL_CPU_ACCEL_X86_SSE2)
2852     {
2853       babl_conversion_new (
2854         babl_format ("RGBA float"),
2855         babl_format ("CIE Lab alpha float"),
2856         "linear", rgbaf_to_Labaf_sse2,
2857         NULL
2858       );
2859       babl_conversion_new (
2860         babl_format ("Y float"),
2861         babl_format ("CIE L float"),
2862         "linear", Yf_to_Lf_sse2,
2863         NULL
2864       );
2865       babl_conversion_new (
2866         babl_format ("YA float"),
2867         babl_format ("CIE L float"),
2868         "linear", Yaf_to_Lf_sse2,
2869         NULL
2870       );
2871       babl_conversion_new (
2872         babl_format ("RGBA float"),
2873         babl_format ("CIE L float"),
2874         "linear", rgbaf_to_Lf_sse2,
2875         NULL
2876       );
2877     }
2878 
2879 #endif /* defined(USE_SSE2) */
2880 
2881   rgbcie_init ();
2882 }
2883 
2884 static void
formats(void)2885 formats (void)
2886 {
2887   babl_format_new (
2888     "name", "CIE Lab float",
2889     babl_model ("CIE Lab"),
2890 
2891     babl_type ("float"),
2892     babl_component ("CIE L"),
2893     babl_component ("CIE a"),
2894     babl_component ("CIE b"),
2895     NULL);
2896 
2897   babl_format_new (
2898     "name", "CIE XYZ float",
2899     babl_model ("CIE XYZ"),
2900 
2901     babl_type ("float"),
2902     babl_component ("CIE X"),
2903     babl_component ("CIE Y"),
2904     babl_component ("CIE Z"),
2905     NULL);
2906 
2907   babl_format_new (
2908     "name", "CIE XYZ alpha float",
2909     babl_model ("CIE XYZ"),
2910 
2911     babl_type ("float"),
2912     babl_component ("CIE X"),
2913     babl_component ("CIE Y"),
2914     babl_component ("CIE Z"),
2915     babl_component ("A"),
2916     NULL);
2917 
2918   babl_format_new (
2919     "name", "CIE Lab alpha float",
2920     babl_model ("CIE Lab alpha"),
2921 
2922     babl_type ("float"),
2923     babl_component ("CIE L"),
2924     babl_component ("CIE a"),
2925     babl_component ("CIE b"),
2926     babl_component ("A"),
2927     NULL);
2928 
2929   babl_format_new (
2930     "name", "CIE LCH(ab) float",
2931     babl_model ("CIE LCH(ab)"),
2932 
2933     babl_type ("float"),
2934     babl_component ("CIE L"),
2935     babl_component ("CIE C(ab)"),
2936     babl_component ("CIE H(ab)"),
2937     NULL);
2938 
2939   babl_format_new (
2940     "name", "CIE LCH(ab) alpha float",
2941     babl_model ("CIE LCH(ab) alpha"),
2942 
2943     babl_type ("float"),
2944     babl_component ("CIE L"),
2945     babl_component ("CIE C(ab)"),
2946     babl_component ("CIE H(ab)"),
2947     babl_component ("A"),
2948     NULL);
2949 
2950   babl_format_new (
2951     "name", "CIE L float",
2952     babl_model ("CIE Lab"),
2953     babl_type ("float"),
2954     babl_component ("CIE L"),
2955     NULL);
2956 
2957   babl_format_new (
2958     "name", "CIE L alpha float",
2959     babl_model ("CIE Lab alpha"),
2960     babl_type ("float"),
2961     babl_component ("CIE L"),
2962     babl_component ("A"),
2963     NULL);
2964 
2965   babl_format_new (
2966     "name", "CIE Lab u8",
2967     babl_model ("CIE Lab"),
2968 
2969     babl_type ("CIE u8 L"),
2970     babl_component ("CIE L"),
2971     babl_type ("CIE u8 ab"),
2972     babl_component ("CIE a"),
2973     babl_type ("CIE u8 ab"),
2974     babl_component ("CIE b"),
2975     NULL);
2976 
2977   babl_format_new (
2978     "name", "CIE Lab u16",
2979     babl_model ("CIE Lab"),
2980 
2981     babl_type ("CIE u16 L"),
2982     babl_component ("CIE L"),
2983     babl_type ("CIE u16 ab"),
2984     babl_component ("CIE a"),
2985     babl_type ("CIE u16 ab"),
2986     babl_component ("CIE b"),
2987     NULL);
2988 
2989   babl_format_new (
2990     "name", "CIE xyY float",
2991     babl_model ("CIE xyY"),
2992 
2993     babl_type ("float"),
2994     babl_component ("CIE x"),
2995     babl_component ("CIE y"),
2996     babl_component ("CIE Y"),
2997     NULL);
2998 
2999   babl_format_new (
3000     "name", "CIE xyY alpha float",
3001     babl_model ("CIE xyY alpha"),
3002 
3003     babl_type ("float"),
3004     babl_component ("CIE x"),
3005     babl_component ("CIE y"),
3006     babl_component ("CIE Y"),
3007     babl_component ("A"),
3008     NULL);
3009 
3010   /* CIE 1976 UCS */
3011   babl_format_new (
3012     "name", "CIE Yuv float",
3013     babl_model ("CIE Yuv"),
3014 
3015     babl_type ("float"),
3016     babl_component ("CIE Y"),
3017     babl_component ("CIE u"),
3018     babl_component ("CIE v"),
3019     NULL);
3020 
3021   babl_format_new (
3022     "name", "CIE Yuv alpha float",
3023     babl_model ("CIE Yuv alpha"),
3024 
3025     babl_type ("float"),
3026     babl_component ("CIE Y"),
3027     babl_component ("CIE u"),
3028     babl_component ("CIE v"),
3029     babl_component ("A"),
3030     NULL);
3031 }
3032 
3033 
3034 /******** end floating point RGB/CIE color space conversions **********/
3035 
3036 /******** begin  integer RGB/CIE color space conversions **************/
3037 
3038 static inline void
convert_double_u8_scaled(const Babl * conversion,double min_val,double max_val,unsigned char min,unsigned char max,char * src,char * dst,int src_pitch,int dst_pitch,long n)3039 convert_double_u8_scaled (const Babl   *conversion,
3040                           double        min_val,
3041                           double        max_val,
3042                           unsigned char min,
3043                           unsigned char max,
3044                           char         *src,
3045                           char         *dst,
3046                           int           src_pitch,
3047                           int           dst_pitch,
3048                           long          n)
3049 {
3050   while (n--)
3051     {
3052       double        dval = *(double *) src;
3053       unsigned char u8val;
3054 
3055       if (dval < min_val)
3056         u8val = min;
3057       else if (dval > max_val)
3058         u8val = max;
3059       else
3060         u8val = rint ((dval - min_val) / (max_val - min_val) * (max - min) + min);
3061 
3062       *(unsigned char *) dst = u8val;
3063       src                   += src_pitch;
3064       dst                   += dst_pitch;
3065     }
3066 }
3067 
3068 static inline void
convert_u8_double_scaled(const Babl * conversion,double min_val,double max_val,unsigned char min,unsigned char max,char * src,char * dst,int src_pitch,int dst_pitch,long n)3069 convert_u8_double_scaled (const Babl   *conversion,
3070                           double        min_val,
3071                           double        max_val,
3072                           unsigned char min,
3073                           unsigned char max,
3074                           char         *src,
3075                           char         *dst,
3076                           int           src_pitch,
3077                           int           dst_pitch,
3078                           long          n)
3079 {
3080   while (n--)
3081     {
3082       int    u8val = *(unsigned char *) src;
3083       double dval;
3084 
3085       if (u8val < min)
3086         dval = min_val;
3087       else if (u8val > max)
3088         dval = max_val;
3089       else
3090         dval = (u8val - min) / (double) (max - min) * (max_val - min_val) + min_val;
3091 
3092       (*(double *) dst) = dval;
3093 
3094       dst += dst_pitch;
3095       src += src_pitch;
3096     }
3097 }
3098 
3099 #define MAKE_CONVERSIONS(name, min_val, max_val, min, max) \
3100   static void \
3101   convert_ ## name ## _double (const Babl *c, char *src, \
3102                                char *dst, \
3103                                int src_pitch, \
3104                                int dst_pitch, \
3105                                long n)        \
3106   { \
3107     convert_u8_double_scaled (c, min_val, max_val, min, max, \
3108                               src, dst, src_pitch, dst_pitch, n); \
3109   }                                                               \
3110   static void \
3111   convert_double_ ## name (const Babl *c, char *src, \
3112                            char *dst, \
3113                            int src_pitch, \
3114                            int dst_pitch, \
3115                            long n)        \
3116   { \
3117     convert_double_u8_scaled (c, min_val, max_val, min, max, \
3118                               src, dst, src_pitch, dst_pitch, n); \
3119   }
3120 
3121 /* source ICC.1:2004-10 */
3122 
3123 MAKE_CONVERSIONS (u8_l, 0.0, 100.0, 0x00, 0xff)
3124 MAKE_CONVERSIONS (u8_ab, -128.0, 127.0, 0x00, 0xff)
3125 
3126 #undef MAKE_CONVERSIONS
3127 
3128 static inline void
convert_float_u8_scaled(const Babl * conversion,float min_val,float max_val,unsigned char min,unsigned char max,char * src,char * dst,int src_pitch,int dst_pitch,long n)3129 convert_float_u8_scaled (const Babl     *conversion,
3130                           float          min_val,
3131                           float          max_val,
3132                           unsigned char  min,
3133                           unsigned char  max,
3134                           char          *src,
3135                           char          *dst,
3136                           int            src_pitch,
3137                           int            dst_pitch,
3138                           long           n)
3139 {
3140   while (n--)
3141     {
3142       float        dval = *(float *) src;
3143       unsigned char u8val;
3144 
3145       if (dval < min_val)
3146         u8val = min;
3147       else if (dval > max_val)
3148         u8val = max;
3149       else
3150         u8val = rint ((dval - min_val) / (max_val - min_val) * (max - min) + min);
3151 
3152       *(unsigned char *) dst = u8val;
3153       src                   += src_pitch;
3154       dst                   += dst_pitch;
3155     }
3156 }
3157 
3158 static inline void
convert_u8_float_scaled(const Babl * conversion,float min_val,float max_val,unsigned char min,unsigned char max,char * src,char * dst,int src_pitch,int dst_pitch,long n)3159 convert_u8_float_scaled (const Babl    *conversion,
3160                           float         min_val,
3161                           float         max_val,
3162                           unsigned char min,
3163                           unsigned char max,
3164                           char         *src,
3165                           char         *dst,
3166                           int           src_pitch,
3167                           int           dst_pitch,
3168                           long          n)
3169 {
3170   while (n--)
3171     {
3172       int    u8val = *(unsigned char *) src;
3173       float dval;
3174 
3175       if (u8val < min)
3176         dval = min_val;
3177       else if (u8val > max)
3178         dval = max_val;
3179       else
3180         dval = (u8val - min) / (float) (max - min) * (max_val - min_val) + min_val;
3181 
3182       (*(float *) dst) = dval;
3183 
3184       dst += dst_pitch;
3185       src += src_pitch;
3186     }
3187 }
3188 
3189 #define MAKE_CONVERSIONS(name, min_val, max_val, min, max) \
3190   static void \
3191   convert_ ## name ## _float (const Babl *c, char *src, \
3192                                char *dst, \
3193                                int src_pitch, \
3194                                int dst_pitch, \
3195                                long n)        \
3196   { \
3197     convert_u8_float_scaled (c, min_val, max_val, min, max, \
3198                               src, dst, src_pitch, dst_pitch, n); \
3199   }                                                               \
3200   static void \
3201   convert_float_ ## name (const Babl *c, char *src, \
3202                            char *dst, \
3203                            int src_pitch, \
3204                            int dst_pitch, \
3205                            long n)        \
3206   { \
3207     convert_float_u8_scaled (c, min_val, max_val, min, max, \
3208                               src, dst, src_pitch, dst_pitch, n); \
3209   }
3210 
3211 /* source ICC.1:2004-10 */
3212 
3213 MAKE_CONVERSIONS (u8_l, 0.0, 100.0, 0x00, 0xff)
3214 MAKE_CONVERSIONS (u8_ab, -128.0, 127.0, 0x00, 0xff)
3215 
3216 #undef MAKE_CONVERSIONS
3217 
3218 static void
types_u8(void)3219 types_u8 (void)
3220 {
3221   babl_type_new (
3222     "CIE u8 L",
3223     "integer",
3224     "unsigned",
3225     "bits", 8,
3226     "min_val", 0.0,
3227     "max_val", 100.0,
3228     NULL
3229   );
3230 
3231   babl_type_new (
3232     "CIE u8 ab",
3233     "integer",
3234     "unsigned",
3235     "bits", 8,
3236     "min_val", -128.0,
3237     "max_val", 127.0,
3238     NULL
3239   );
3240 
3241   babl_conversion_new (
3242     babl_type ("CIE u8 L"),
3243     babl_type ("double"),
3244     "plane", convert_u8_l_double,
3245     NULL
3246   );
3247   babl_conversion_new (
3248     babl_type ("double"),
3249     babl_type ("CIE u8 L"),
3250     "plane", convert_double_u8_l,
3251     NULL
3252   );
3253 
3254   babl_conversion_new (
3255     babl_type ("CIE u8 ab"),
3256     babl_type ("double"),
3257     "plane", convert_u8_ab_double,
3258     NULL
3259   );
3260   babl_conversion_new (
3261     babl_type ("double"),
3262     babl_type ("CIE u8 ab"),
3263     "plane", convert_double_u8_ab,
3264     NULL
3265   );
3266 
3267   babl_conversion_new (
3268     babl_type ("CIE u8 L"),
3269     babl_type ("float"),
3270     "plane", convert_u8_l_float,
3271     NULL
3272   );
3273   babl_conversion_new (
3274     babl_type ("float"),
3275     babl_type ("CIE u8 L"),
3276     "plane", convert_float_u8_l,
3277     NULL
3278   );
3279 
3280   babl_conversion_new (
3281     babl_type ("CIE u8 ab"),
3282     babl_type ("float"),
3283     "plane", convert_u8_ab_float,
3284     NULL
3285   );
3286   babl_conversion_new (
3287     babl_type ("float"),
3288     babl_type ("CIE u8 ab"),
3289     "plane", convert_float_u8_ab,
3290     NULL
3291   );
3292 }
3293 
3294 static inline void
convert_double_u16_scaled(const Babl * conversion,double min_val,double max_val,unsigned short min,unsigned short max,char * src,char * dst,int src_pitch,int dst_pitch,long n)3295 convert_double_u16_scaled (const Babl    *conversion,
3296                            double         min_val,
3297                            double         max_val,
3298                            unsigned short min,
3299                            unsigned short max,
3300                            char          *src,
3301                            char          *dst,
3302                            int            src_pitch,
3303                            int            dst_pitch,
3304                            long           n)
3305 {
3306   while (n--)
3307     {
3308       double         dval = *(double *) src;
3309       unsigned short u16val;
3310 
3311       if (dval < min_val)
3312         u16val = min;
3313       else if (dval > max_val)
3314         u16val = max;
3315       else
3316         u16val = rint ((dval - min_val) / (max_val - min_val) * (max - min) + min);
3317 
3318       *(unsigned short *) dst = u16val;
3319       dst                    += dst_pitch;
3320       src                    += src_pitch;
3321     }
3322 }
3323 
3324 static inline void
convert_u16_double_scaled(const Babl * conversion,double min_val,double max_val,unsigned short min,unsigned short max,char * src,char * dst,int src_pitch,int dst_pitch,long n)3325 convert_u16_double_scaled (const Babl    *conversion,
3326                            double         min_val,
3327                            double         max_val,
3328                            unsigned short min,
3329                            unsigned short max,
3330                            char          *src,
3331                            char          *dst,
3332                            int            src_pitch,
3333                            int            dst_pitch,
3334                            long           n)
3335 {
3336   while (n--)
3337     {
3338       int    u16val = *(unsigned short *) src;
3339       double dval;
3340 
3341       if (u16val < min)
3342         dval = min_val;
3343       else if (u16val > max)
3344         dval = max_val;
3345       else
3346         dval = (u16val - min) / (double) (max - min) * (max_val - min_val) + min_val;
3347 
3348       (*(double *) dst) = dval;
3349       dst              += dst_pitch;
3350       src              += src_pitch;
3351     }
3352 }
3353 
3354 #define MAKE_CONVERSIONS(name, min_val, max_val, min, max)      \
3355   static void \
3356   convert_ ## name ## _double (const Babl *c, char *src, \
3357                                char *dst, \
3358                                int src_pitch, \
3359                                int dst_pitch, \
3360                                long n)        \
3361   { \
3362     convert_u16_double_scaled (c, min_val, max_val, min, max, \
3363                                src, dst, src_pitch, dst_pitch, n); \
3364   }                                                               \
3365   static void \
3366   convert_double_ ## name (const Babl *c, char *src, \
3367                            char *dst, \
3368                            int src_pitch, \
3369                            int dst_pitch, \
3370                            long n)        \
3371   { \
3372     convert_double_u16_scaled (c, min_val, max_val, min, max, \
3373                                src, dst, src_pitch, dst_pitch, n); \
3374   }
3375 
3376 MAKE_CONVERSIONS (u16_l, 0.0, 100.0, 0x00, 0xffff)
3377 MAKE_CONVERSIONS (u16_ab, -128.0, 127.0, 0x00, 0xffff)
3378 
3379 #undef MAKE_CONVERSIONS
3380 
3381 
3382 static inline void
convert_float_u16_scaled(const Babl * conversion,float min_val,float max_val,unsigned short min,unsigned short max,char * src,char * dst,int src_pitch,int dst_pitch,long n)3383 convert_float_u16_scaled (const Babl     *conversion,
3384                            float          min_val,
3385                            float          max_val,
3386                            unsigned short min,
3387                            unsigned short max,
3388                            char          *src,
3389                            char          *dst,
3390                            int            src_pitch,
3391                            int            dst_pitch,
3392                            long           n)
3393 {
3394   while (n--)
3395     {
3396       float         dval = *(float *) src;
3397       unsigned short u16val;
3398 
3399       if (dval < min_val)
3400         u16val = min;
3401       else if (dval > max_val)
3402         u16val = max;
3403       else
3404         u16val = rint ((dval - min_val) / (max_val - min_val) * (max - min) + min);
3405 
3406       *(unsigned short *) dst = u16val;
3407       dst                    += dst_pitch;
3408       src                    += src_pitch;
3409     }
3410 }
3411 
3412 static inline void
convert_u16_float_scaled(const Babl * conversion,float min_val,float max_val,unsigned short min,unsigned short max,char * src,char * dst,int src_pitch,int dst_pitch,long n)3413 convert_u16_float_scaled (const Babl     *conversion,
3414                            float          min_val,
3415                            float          max_val,
3416                            unsigned short min,
3417                            unsigned short max,
3418                            char          *src,
3419                            char          *dst,
3420                            int            src_pitch,
3421                            int            dst_pitch,
3422                            long           n)
3423 {
3424   while (n--)
3425     {
3426       int    u16val = *(unsigned short *) src;
3427       float dval;
3428 
3429       if (u16val < min)
3430         dval = min_val;
3431       else if (u16val > max)
3432         dval = max_val;
3433       else
3434         dval = (u16val - min) / (float) (max - min) * (max_val - min_val) + min_val;
3435 
3436       (*(float *) dst) = dval;
3437       dst              += dst_pitch;
3438       src              += src_pitch;
3439     }
3440 }
3441 
3442 #define MAKE_CONVERSIONS(name, min_val, max_val, min, max)      \
3443   static void \
3444   convert_ ## name ## _float (const Babl *c, char *src, \
3445                                char *dst, \
3446                                int src_pitch, \
3447                                int dst_pitch, \
3448                                long n)        \
3449   { \
3450     convert_u16_float_scaled (c, min_val, max_val, min, max, \
3451                                src, dst, src_pitch, dst_pitch, n); \
3452   }                                                               \
3453   static void \
3454   convert_float_ ## name (const Babl *c, char *src, \
3455                            char *dst, \
3456                            int src_pitch, \
3457                            int dst_pitch, \
3458                            long n)        \
3459   { \
3460     convert_float_u16_scaled (c, min_val, max_val, min, max, \
3461                                src, dst, src_pitch, dst_pitch, n); \
3462   }
3463 
3464 MAKE_CONVERSIONS (u16_l, 0.0, 100.0, 0x00, 0xffff)
3465 MAKE_CONVERSIONS (u16_ab, -128.0, 127.0, 0x00, 0xffff)
3466 
3467 #undef MAKE_CONVERSIONS
3468 
3469 static void
types_u16(void)3470 types_u16 (void)
3471 {
3472   babl_type_new (
3473     "CIE u16 L",
3474     "integer",
3475     "unsigned",
3476     "bits", 16,
3477     "min_val", 0.0,
3478     "max_val", 100.0,
3479     NULL
3480   );
3481 
3482   babl_type_new (
3483     "CIE u16 ab",
3484     "integer",
3485     "unsigned",
3486     "bits", 16,
3487     "min_val", -128.0,
3488     "max_val", 127.0,
3489     NULL
3490   );
3491 
3492 
3493   babl_conversion_new (
3494     babl_type ("CIE u16 L"),
3495     babl_type ("double"),
3496     "plane", convert_u16_l_double,
3497     NULL
3498   );
3499   babl_conversion_new (
3500     babl_type ("double"),
3501     babl_type ("CIE u16 L"),
3502     "plane", convert_double_u16_l,
3503     NULL
3504   );
3505 
3506   babl_conversion_new (
3507     babl_type ("CIE u16 ab"),
3508     babl_type ("double"),
3509     "plane", convert_u16_ab_double,
3510     NULL
3511   );
3512   babl_conversion_new (
3513     babl_type ("double"),
3514     babl_type ("CIE u16 ab"),
3515     "plane", convert_double_u16_ab,
3516     NULL
3517   );
3518 
3519   babl_conversion_new (
3520     babl_type ("CIE u16 L"),
3521     babl_type ("float"),
3522     "plane", convert_u16_l_float,
3523     NULL
3524   );
3525   babl_conversion_new (
3526     babl_type ("float"),
3527     babl_type ("CIE u16 L"),
3528     "plane", convert_float_u16_l,
3529     NULL
3530   );
3531 
3532   babl_conversion_new (
3533     babl_type ("CIE u16 ab"),
3534     babl_type ("float"),
3535     "plane", convert_u16_ab_float,
3536     NULL
3537   );
3538   babl_conversion_new (
3539     babl_type ("float"),
3540     babl_type ("CIE u16 ab"),
3541     "plane", convert_float_u16_ab,
3542     NULL
3543   );
3544 }
3545 
3546 static void
types(void)3547 types (void)
3548 {
3549   types_u8 ();
3550   types_u16 ();
3551 }
3552 
3553 /******** end  integer RGB/CIE color space conversions ****************/
3554 
3555 static void
rgbxyzrgb_init(void)3556 rgbxyzrgb_init (void)
3557 {
3558 }
3559 
3560 static void
rgbcie_init(void)3561 rgbcie_init (void)
3562 {
3563   static int initialized = 0;
3564 
3565   if (!initialized)
3566     {
3567       rgbxyzrgb_init ();
3568       initialized = 1;
3569     }
3570 }
3571