1/*
2    This file is part of darktable,
3    copyright (c) 2021 darktable developers.
4
5    darktable is free software: you can redistribute it and/or modify
6    it under the terms of the GNU General Public License as published by
7    the Free Software Foundation, either version 3 of the License, or
8    (at your option) any later version.
9
10    darktable is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU General Public License for more details.
14
15    You should have received a copy of the GNU General Public License
16    along with darktable.  If not, see <http://www.gnu.org/licenses/>.
17*/
18
19#include "colorspace.h"
20#include "common.h"
21
22typedef enum dt_iop_channelmixer_rgb_version_t
23{
24  CHANNELMIXERRGB_V_1 = 0, // $DESCRIPTION: "version 1 (2020)"
25  CHANNELMIXERRGB_V_2 = 1, // $DESCRIPTION: "version 2 (2021)"
26  CHANNELMIXERRGB_V_3 = 2, // $DESCRIPTION: "version 3 (Apr 2021)"
27} dt_iop_channelmixer_rgb_version_t;
28
29
30typedef enum dt_adaptation_t
31{
32  DT_ADAPTATION_LINEAR_BRADFORD = 0, // $DESCRIPTION: "linear Bradford (ICC v4)"
33  DT_ADAPTATION_CAT16           = 1, // $DESCRIPTION: "CAT16 (CIECAM16)"
34  DT_ADAPTATION_FULL_BRADFORD   = 2, // $DESCRIPTION: "non-linear Bradford"
35  DT_ADAPTATION_XYZ             = 3, // $DESCRIPTION: "XYZ "
36  DT_ADAPTATION_RGB             = 4, // $DESCRIPTION: "none (bypass)"
37  DT_ADAPTATION_LAST
38} dt_adaptation_t;
39
40
41#define INVERSE_SQRT_3 0.5773502691896258f
42#define TRUE 1
43#define FALSE 0
44#define NORM_MIN 1e-6f
45
46static inline float sqf(const float x)
47{
48  return x * x;
49}
50
51static inline float euclidean_norm(const float4 input)
52{
53  return fmax(native_sqrt(sqf(input.x) + sqf(input.y) + sqf(input.z)), NORM_MIN);
54}
55
56static inline float4 gamut_mapping(const float4 input, const float compression, const int clip)
57{
58  // Get the sum XYZ
59  const float sum = input.x + input.y + input.z;
60  const float Y = input.y;
61
62  float4 output;
63
64  if(sum > 0.f && Y > 0.f)
65  {
66    // Convert to xyY
67    float4 xyY = { input.x / sum, input.y / sum , Y, 0.0f };
68
69    // Convert to uvY
70    float4 uvY = dt_xyY_to_uvY(xyY);
71
72    // Get the chromaticity difference with white point uv
73    const float2 D50 = { 0.20915914598542354f, 0.488075320769787f };
74    const float2 delta = D50 - uvY.xy;
75    const float Delta = Y * (sqf(delta.x) + sqf(delta.y));
76
77    // Compress chromaticity (move toward white point)
78    const float correction = (compression == 0.0f) ? 0.f : native_powr(Delta, compression);
79
80    // Ensure the correction does not bring our uyY vector the other side of D50
81    // that would switch to the opposite color, so we clip at D50
82    const float2 tmp =  correction * delta + uvY.xy;
83    uvY.xy = (uvY.xy > D50) ? fmax(tmp, D50)
84                            : fmin(tmp, D50);
85
86    // Convert back to xyY
87    xyY = dt_uvY_to_xyY(uvY);
88
89    // Clip upon request
90    if(clip) xyY.xy = fmax(xyY.xy, 0.0f);
91
92    // Check sanity of y
93    // since we later divide by y, it can't be zero
94    xyY.y = fmax(xyY.y, NORM_MIN);
95
96    // Check sanity of x and y :
97    // since Z = Y (1 - x - y) / y, if x + y >= 1, Z will be negative
98    const float scale = xyY.x + xyY.y;
99    const int sanitize = (scale >= 1.f);
100    xyY.xy = (sanitize) ? xyY.xy / scale : xyY.xy;
101
102    // Convert back to XYZ
103    output = dt_xyY_to_XYZ(xyY);
104  }
105  else
106  {
107    // sum of channels == 0, and/or Y == 0 so we have black
108    output = (float4)0.f;
109  }
110  return output;
111}
112
113static inline float4 luma_chroma(const float4 input, const float4 saturation, const float4 lightness,
114                                 const dt_iop_channelmixer_rgb_version_t version)
115{
116  float4 output;
117
118  // Compute euclidean norm
119  float norm = euclidean_norm(input);
120  const float avg = fmax((input.x + input.y + input.z) / 3.0f, NORM_MIN);
121
122  if(norm > 0.f && avg > 0.f)
123  {
124    // Compute flat lightness adjustment
125    const float mix = dot(input, lightness);
126
127    // Compensate the norm to get color ratios (R, G, B) = (1, 1, 1) for grey (colorless) pixels.
128    if(version == CHANNELMIXERRGB_V_3) norm *= INVERSE_SQRT_3;
129
130    // Ratios
131    // WARNING : dot product below uses all 4 channels, you need to make sure
132    // input.w != NaN since saturation.w = 0.f
133    output = input / norm;
134
135    // Compute ratios and a flat colorfulness adjustment for the whole pixel
136    float coeff_ratio = 0.f;
137
138    if(version == CHANNELMIXERRGB_V_1)
139      coeff_ratio = dot((1.f - output), saturation);
140    else
141      coeff_ratio = dot(output, saturation) / 3.f;
142
143    // Adjust the RGB ratios with the pixel correction
144
145    // if the ratio was already invalid (negative), we accept the result to be invalid too
146    // otherwise bright saturated blues end up solid black
147    const float4 min_ratio = (output < 0.0f) ? output : 0.0f;
148    const float4 output_inverse = 1.0f - output;
149    output = fmax(output_inverse * coeff_ratio + output, min_ratio);
150
151    // The above interpolation between original pixel ratios and (1, 1, 1) might change the norm of the
152    // ratios. Compensate for that.
153    if(version == CHANNELMIXERRGB_V_3) norm /= euclidean_norm(output) * INVERSE_SQRT_3;
154
155    // Apply colorfulness adjustment channel-wise and repack with lightness to get LMS back
156    norm *= fmax(1.f + mix / avg, 0.f);
157    output *= norm;
158  }
159  else
160  {
161    // we have black, 0 stays 0, no luminance = no color
162    output = input;
163  }
164
165  return output;
166}
167
168#define unswitch_convert_any_LMS_to_XYZ(kind) \
169  ({ switch(kind) \
170    { \
171      case DT_ADAPTATION_FULL_BRADFORD: \
172      case DT_ADAPTATION_LINEAR_BRADFORD: \
173      { \
174        XYZ = convert_bradford_LMS_to_XYZ(LMS); \
175        break; \
176      }\
177      case DT_ADAPTATION_CAT16:\
178      { \
179        XYZ = convert_CAT16_LMS_to_XYZ(LMS); \
180        break; \
181      } \
182      case DT_ADAPTATION_XYZ: \
183      { \
184        XYZ = LMS; \
185        break; \
186      } \
187      case DT_ADAPTATION_RGB: \
188      case DT_ADAPTATION_LAST: \
189      { \
190        XYZ = matrix_product_float4(LMS, RGB_to_XYZ); \
191        break; \
192      } \
193    }})
194
195
196#define unswitch_convert_XYZ_to_any_LMS(kind) \
197  ({ switch(kind) \
198  { \
199    case DT_ADAPTATION_FULL_BRADFORD: \
200    case DT_ADAPTATION_LINEAR_BRADFORD: \
201    { \
202      LMS = convert_XYZ_to_bradford_LMS(XYZ); \
203      break; \
204    } \
205    case DT_ADAPTATION_CAT16: \
206    { \
207      LMS = convert_XYZ_to_CAT16_LMS(XYZ); \
208      break; \
209    } \
210    case DT_ADAPTATION_XYZ: \
211    { \
212      LMS = XYZ; \
213      break; \
214    } \
215    case DT_ADAPTATION_RGB: \
216    case DT_ADAPTATION_LAST: \
217    { \
218      LMS = matrix_product_float4(XYZ, XYZ_to_RGB); \
219      break; \
220    } \
221  }})
222
223
224static inline void downscale_vector(float4 *const vector, const float scaling)
225{
226  const int valid = (scaling > NORM_MIN) && !isnan(scaling);
227  *vector /= (valid) ? (scaling + NORM_MIN) : NORM_MIN;
228}
229
230static inline void upscale_vector(float4 *const vector, const float scaling)
231{
232  const int valid = (scaling > NORM_MIN) && !isnan(scaling);
233  *vector *= (valid) ? (scaling + NORM_MIN) : NORM_MIN;
234}
235
236static inline float4 chroma_adapt_bradford(const float4 RGB,
237                                           constant const float *const RGB_to_XYZ,
238                                           constant const float *const MIX,
239                                           const float4 illuminant, const float p,
240                                           const int full)
241{
242  // Convert from RGB to XYZ
243  const float4 XYZ = matrix_product_float4(RGB, RGB_to_XYZ);
244  const float Y = XYZ.y;
245
246  // Convert to LMS
247  float4 LMS = convert_XYZ_to_bradford_LMS(XYZ);
248
249  // Do white balance
250  downscale_vector(&LMS, Y);
251    bradford_adapt_D50(&LMS, illuminant, p, full);
252  upscale_vector(&LMS, Y);
253
254  // Compute the 3D mix - this is a rotation + homothety of the vector base
255  const float4 LMS_mixed = matrix_product_float4(LMS, MIX);
256
257  return convert_bradford_LMS_to_XYZ(LMS_mixed);
258};
259
260static inline float4 chroma_adapt_CAT16(const float4 RGB,
261                                        constant const float *const RGB_to_XYZ,
262                                        constant const float *const MIX,
263                                        const float4 illuminant, const float p,
264                                        const int full)
265{
266  // Convert from RGB to XYZ
267  const float4 XYZ = matrix_product_float4(RGB, RGB_to_XYZ);
268  const float Y = XYZ.y;
269
270  // Convert to LMS
271  float4 LMS = convert_XYZ_to_CAT16_LMS(XYZ);
272
273  // Do white balance
274  downscale_vector(&LMS, Y);
275    CAT16_adapt_D50(&LMS, illuminant, p, full); // force full-adaptation
276  upscale_vector(&LMS, Y);
277
278  // Compute the 3D mix - this is a rotation + homothety of the vector base
279  const float4 LMS_mixed = matrix_product_float4(LMS, MIX);
280
281  return convert_CAT16_LMS_to_XYZ(LMS_mixed);
282}
283
284static inline float4 chroma_adapt_XYZ(const float4 RGB,
285                                      constant const float *const RGB_to_XYZ,
286                                      constant const float *const MIX, const float4 illuminant)
287{
288  // Convert from RGB to XYZ
289  float4 XYZ_mixed = matrix_product_float4(RGB, RGB_to_XYZ);
290  const float Y = XYZ_mixed.y;
291
292  // Do white balance in XYZ
293  downscale_vector(&XYZ_mixed, Y);
294    XYZ_adapt_D50(&XYZ_mixed, illuminant);
295  upscale_vector(&XYZ_mixed, Y);
296
297  // Compute the 3D mix in XYZ - this is a rotation + homothety of the vector base
298  return matrix_product_float4(XYZ_mixed, MIX);
299}
300
301static inline float4 chroma_adapt_RGB(const float4 RGB,
302                               constant const float *const RGB_to_XYZ,
303                               constant const float *const MIX)
304{
305  // No white balance.
306
307  // Compute the 3D mix in RGB - this is a rotation + homothety of the vector base
308  float4 RGB_mixed = matrix_product_float4(RGB, MIX);
309
310  // Convert from RGB to XYZ
311  return matrix_product_float4(RGB_mixed, RGB_to_XYZ);
312}
313
314#define unswitch_chroma_adapt(kind) \
315 ({ switch(kind) \
316  { \
317    case DT_ADAPTATION_FULL_BRADFORD: \
318    { \
319      XYZ = chroma_adapt_bradford(RGB, RGB_to_XYZ, MIX, illuminant, p, TRUE); \
320      break; \
321    } \
322    case DT_ADAPTATION_LINEAR_BRADFORD: \
323    { \
324      XYZ = chroma_adapt_bradford(RGB, RGB_to_XYZ, MIX, illuminant, p, FALSE); \
325      break; \
326    } \
327    case DT_ADAPTATION_CAT16: \
328    { \
329      XYZ = chroma_adapt_CAT16(RGB, RGB_to_XYZ, MIX, illuminant, 1.f, FALSE); \
330      break; \
331    } \
332    case DT_ADAPTATION_XYZ: \
333    { \
334      XYZ = chroma_adapt_XYZ(RGB, RGB_to_XYZ, MIX, illuminant); \
335      break; \
336    } \
337    case DT_ADAPTATION_RGB: \
338    case DT_ADAPTATION_LAST: \
339    { \
340      XYZ = chroma_adapt_RGB(RGB, RGB_to_XYZ, MIX); \
341      break; \
342    } \
343  }})
344
345
346/*
347* The following kernels are 100% copy-pasted with the exception of
348* the first line : const dt_adaptation_t kind = ...
349* This ensures to unswitch the color space conversions branches for performance
350* while keeping the same overall code structure for maintenance.
351*
352* The reference C version in src/iop/channelmixerrgb.c does it differently
353* since C has an explicite -funswitchloop option, but OpenCL doesn't and
354* we have to do it manually using macros and duplicating kernels.
355*/
356
357kernel void
358channelmixerrgb_CAT16(read_only image2d_t in, write_only image2d_t out,
359                      const int width, const int height,
360                      constant const float *const RGB_to_XYZ,
361                      constant const float *const XYZ_to_RGB,
362                      constant const float *const MIX,
363                      const float4 illuminant,
364                      const float4 saturation,
365                      const float4 lightness,
366                      const float4 grey,
367                      const float p, const float gamut, const int clip, const int apply_grey,
368                      const dt_iop_channelmixer_rgb_version_t version)
369{
370  const dt_adaptation_t kind = DT_ADAPTATION_CAT16;
371
372  const int x = get_global_id(0);
373  const int y = get_global_id(1);
374
375  if(x >= width || y >= height) return;
376
377  float4 pix_in = read_imagef(in, sampleri, (int2)(x, y));
378
379  float4 XYZ, LMS;
380  float4 RGB = pix_in;
381  RGB.w = 0.f;
382
383  if(clip) RGB = fmax(RGB, 0.f);
384
385  /* WE START IN PIPELINE RGB */
386  unswitch_chroma_adapt(kind);
387
388  /* FROM HERE WE ARE MANDATORILY IN XYZ - DATA IS IN temp_one */
389
390  // Gamut mapping happens in XYZ space no matter what
391  XYZ = gamut_mapping(XYZ, gamut, clip);
392
393  // convert to LMS, XYZ or pipeline RGB
394  unswitch_convert_XYZ_to_any_LMS(kind);
395
396  /* FROM HERE WE ARE IN LMS, XYZ OR PIPELINE RGB depending on user param */
397
398  // Clip in LMS
399  if(clip) LMS = fmax(LMS, 0.0f);
400
401  // Apply lightness / saturation adjustment
402  LMS = luma_chroma(LMS, saturation, lightness, version);
403
404  // Clip in LMS
405  if(clip) LMS = fmax(LMS, 0.0f);
406
407  // Save
408  if(apply_grey)
409  {
410    // Turn LMS, XYZ or pipeline RGB into monochrome
411    const float grey_mix = fmax(dot(LMS, grey), 0.0f);
412
413    RGB.xyz = grey_mix;
414    RGB.w = pix_in.w; // alpha mask
415  }
416  else
417  {
418    // Convert back to XYZ
419    unswitch_convert_any_LMS_to_XYZ(kind);
420
421    /* FROM HERE WE ARE MANDATORILY IN XYZ */
422
423    // Clip in XYZ
424    if(clip) XYZ = fmax(XYZ, 0.0f);
425
426    // Convert back to RGB
427    RGB = matrix_product_float4(XYZ, XYZ_to_RGB);
428
429    if(clip) RGB = fmax(RGB, 0.f);
430    RGB.w = pix_in.w;
431  }
432
433  write_imagef(out, (int2)(x, y), RGB);
434}
435
436
437kernel void
438channelmixerrgb_bradford_linear(read_only image2d_t in, write_only image2d_t out,
439                                const int width, const int height,
440                                constant const float *const RGB_to_XYZ,
441                                constant const float *const XYZ_to_RGB,
442                                constant const float *const MIX,
443                                const float4 illuminant,
444                                const float4 saturation,
445                                const float4 lightness,
446                                const float4 grey,
447                                const float p, const float gamut, const int clip, const int apply_grey,
448                                const dt_iop_channelmixer_rgb_version_t version)
449{
450  const dt_adaptation_t kind = DT_ADAPTATION_LINEAR_BRADFORD;
451
452  const int x = get_global_id(0);
453  const int y = get_global_id(1);
454
455  if(x >= width || y >= height) return;
456
457  float4 pix_in = read_imagef(in, sampleri, (int2)(x, y));
458
459  float4 XYZ, LMS;
460  float4 RGB = pix_in;
461  RGB.w = 0.f;
462
463  if(clip) RGB = fmax(RGB, 0.f);
464
465  /* WE START IN PIPELINE RGB */
466  unswitch_chroma_adapt(kind);
467
468  /* FROM HERE WE ARE MANDATORILY IN XYZ - DATA IS IN temp_one */
469
470  // Gamut mapping happens in XYZ space no matter what
471  XYZ = gamut_mapping(XYZ, gamut, clip);
472
473  // convert to LMS, XYZ or pipeline RGB
474  unswitch_convert_XYZ_to_any_LMS(kind);
475
476  /* FROM HERE WE ARE IN LMS, XYZ OR PIPELINE RGB depending on user param */
477
478  // Clip in LMS
479  if(clip) LMS = fmax(LMS, 0.0f);
480
481  // Apply lightness / saturation adjustment
482  LMS = luma_chroma(LMS, saturation, lightness, version);
483
484  // Clip in LMS
485  if(clip) LMS = fmax(LMS, 0.0f);
486
487  // Save
488  if(apply_grey)
489  {
490    // Turn LMS, XYZ or pipeline RGB into monochrome
491    const float grey_mix = fmax(dot(LMS, grey), 0.0f);
492
493    RGB.xyz = grey_mix;
494    RGB.w = pix_in.w; // alpha mask
495  }
496  else
497  {
498    // Convert back to XYZ
499    unswitch_convert_any_LMS_to_XYZ(kind);
500
501    /* FROM HERE WE ARE MANDATORILY IN XYZ */
502
503    // Clip in XYZ
504    if(clip) XYZ = fmax(XYZ, 0.0f);
505
506    // Convert back to RGB
507    RGB = matrix_product_float4(XYZ, XYZ_to_RGB);
508
509    if(clip) RGB = fmax(RGB, 0.f);
510    RGB.w = pix_in.w;
511  }
512
513  write_imagef(out, (int2)(x, y), RGB);
514}
515
516kernel void
517channelmixerrgb_bradford_full(read_only image2d_t in, write_only image2d_t out,
518                              const int width, const int height,
519                              constant const float *const RGB_to_XYZ,
520                              constant const float *const XYZ_to_RGB,
521                              constant const float *const MIX,
522                              const float4 illuminant,
523                              const float4 saturation,
524                              const float4 lightness,
525                              const float4 grey,
526                              const float p, const float gamut, const int clip, const int apply_grey,
527                              const dt_iop_channelmixer_rgb_version_t version)
528{
529  const dt_adaptation_t kind = DT_ADAPTATION_FULL_BRADFORD;
530
531  const int x = get_global_id(0);
532  const int y = get_global_id(1);
533
534  if(x >= width || y >= height) return;
535
536  float4 pix_in = read_imagef(in, sampleri, (int2)(x, y));
537
538  float4 XYZ, LMS;
539  float4 RGB = pix_in;
540  RGB.w = 0.f;
541
542  if(clip) RGB = fmax(RGB, 0.f);
543
544  /* WE START IN PIPELINE RGB */
545  unswitch_chroma_adapt(kind);
546
547  /* FROM HERE WE ARE MANDATORILY IN XYZ - DATA IS IN temp_one */
548
549  // Gamut mapping happens in XYZ space no matter what
550  XYZ = gamut_mapping(XYZ, gamut, clip);
551
552  // convert to LMS, XYZ or pipeline RGB
553  unswitch_convert_XYZ_to_any_LMS(kind);
554
555  /* FROM HERE WE ARE IN LMS, XYZ OR PIPELINE RGB depending on user param */
556
557  // Clip in LMS
558  if(clip) LMS = fmax(LMS, 0.0f);
559
560  // Apply lightness / saturation adjustment
561  LMS = luma_chroma(LMS, saturation, lightness, version);
562
563  // Clip in LMS
564  if(clip) LMS = fmax(LMS, 0.0f);
565
566  // Save
567  if(apply_grey)
568  {
569    // Turn LMS, XYZ or pipeline RGB into monochrome
570    const float grey_mix = fmax(dot(LMS, grey), 0.0f);
571
572    RGB.xyz = grey_mix;
573    RGB.w = pix_in.w; // alpha mask
574  }
575  else
576  {
577    // Convert back to XYZ
578    unswitch_convert_any_LMS_to_XYZ(kind);
579
580    /* FROM HERE WE ARE MANDATORILY IN XYZ */
581
582    // Clip in XYZ
583    if(clip) XYZ = fmax(XYZ, 0.0f);
584
585    // Convert back to RGB
586    RGB = matrix_product_float4(XYZ, XYZ_to_RGB);
587
588    if(clip) RGB = fmax(RGB, 0.f);
589    RGB.w = pix_in.w;
590  }
591
592  write_imagef(out, (int2)(x, y), RGB);
593}
594
595
596kernel void
597channelmixerrgb_XYZ(read_only image2d_t in, write_only image2d_t out,
598                    const int width, const int height,
599                    constant const float *const RGB_to_XYZ,
600                    constant const float *const XYZ_to_RGB,
601                    constant const float *const MIX,
602                    const float4 illuminant,
603                    const float4 saturation,
604                    const float4 lightness,
605                    const float4 grey,
606                    const float p, const float gamut, const int clip, const int apply_grey,
607                    const dt_iop_channelmixer_rgb_version_t version)
608{
609  const dt_adaptation_t kind = DT_ADAPTATION_XYZ;
610
611  const int x = get_global_id(0);
612  const int y = get_global_id(1);
613
614  if(x >= width || y >= height) return;
615
616  float4 pix_in = read_imagef(in, sampleri, (int2)(x, y));
617
618  float4 XYZ, LMS;
619  float4 RGB = pix_in;
620  RGB.w = 0.f;
621
622  if(clip) RGB = fmax(RGB, 0.f);
623
624  /* WE START IN PIPELINE RGB */
625  unswitch_chroma_adapt(kind);
626
627  /* FROM HERE WE ARE MANDATORILY IN XYZ - DATA IS IN temp_one */
628
629  // Gamut mapping happens in XYZ space no matter what
630  XYZ = gamut_mapping(XYZ, gamut, clip);
631
632  // convert to LMS, XYZ or pipeline RGB
633  unswitch_convert_XYZ_to_any_LMS(kind);
634
635  /* FROM HERE WE ARE IN LMS, XYZ OR PIPELINE RGB depending on user param */
636
637  // Clip in LMS
638  if(clip) LMS = fmax(LMS, 0.0f);
639
640  // Apply lightness / saturation adjustment
641  LMS = luma_chroma(LMS, saturation, lightness, version);
642
643  // Clip in LMS
644  if(clip) LMS = fmax(LMS, 0.0f);
645
646  // Save
647  if(apply_grey)
648  {
649    // Turn LMS, XYZ or pipeline RGB into monochrome
650    const float grey_mix = fmax(dot(LMS, grey), 0.0f);
651
652    RGB.xyz = grey_mix;
653    RGB.w = pix_in.w; // alpha mask
654  }
655  else
656  {
657    // Convert back to XYZ
658    unswitch_convert_any_LMS_to_XYZ(kind);
659
660    /* FROM HERE WE ARE MANDATORILY IN XYZ */
661
662    // Clip in XYZ
663    if(clip) XYZ = fmax(XYZ, 0.0f);
664
665    // Convert back to RGB
666    RGB = matrix_product_float4(XYZ, XYZ_to_RGB);
667
668    if(clip) RGB = fmax(RGB, 0.f);
669    RGB.w = pix_in.w;
670  }
671
672  write_imagef(out, (int2)(x, y), RGB);
673}
674
675kernel void
676channelmixerrgb_RGB(read_only image2d_t in, write_only image2d_t out,
677                    const int width, const int height,
678                    constant const float *const RGB_to_XYZ,
679                    constant const float *const XYZ_to_RGB,
680                    constant const float *const MIX,
681                    const float4 illuminant,
682                    const float4 saturation,
683                    const float4 lightness,
684                    const float4 grey,
685                    const float p, const float gamut, const int clip, const int apply_grey,
686                    const dt_iop_channelmixer_rgb_version_t version)
687{
688  const dt_adaptation_t kind = DT_ADAPTATION_RGB;
689
690  const int x = get_global_id(0);
691  const int y = get_global_id(1);
692
693  if(x >= width || y >= height) return;
694
695  float4 pix_in = read_imagef(in, sampleri, (int2)(x, y));
696
697  float4 XYZ, LMS;
698  float4 RGB = pix_in;
699  RGB.w = 0.f;
700
701  if(clip) RGB = fmax(RGB, 0.f);
702
703  /* WE START IN PIPELINE RGB */
704  unswitch_chroma_adapt(kind);
705
706  /* FROM HERE WE ARE MANDATORILY IN XYZ - DATA IS IN temp_one */
707
708  // Gamut mapping happens in XYZ space no matter what
709  XYZ = gamut_mapping(XYZ, gamut, clip);
710
711  // convert to LMS, XYZ or pipeline RGB
712  unswitch_convert_XYZ_to_any_LMS(kind);
713
714  /* FROM HERE WE ARE IN LMS, XYZ OR PIPELINE RGB depending on user param */
715
716  // Clip in LMS
717  if(clip) LMS = fmax(LMS, 0.0f);
718
719  // Apply lightness / saturation adjustment
720  LMS = luma_chroma(LMS, saturation, lightness, version);
721
722  // Clip in LMS
723  if(clip) LMS = fmax(LMS, 0.0f);
724
725  // Save
726  if(apply_grey)
727  {
728    // Turn LMS, XYZ or pipeline RGB into monochrome
729    const float grey_mix = fmax(dot(LMS, grey), 0.0f);
730
731    RGB.xyz = grey_mix;
732    RGB.w = pix_in.w; // alpha mask
733  }
734  else
735  {
736    // Convert back to XYZ
737    unswitch_convert_any_LMS_to_XYZ(kind);
738
739    /* FROM HERE WE ARE MANDATORILY IN XYZ */
740
741    // Clip in XYZ
742    if(clip) XYZ = fmax(XYZ, 0.0f);
743
744    // Convert back to RGB
745    RGB = matrix_product_float4(XYZ, XYZ_to_RGB);
746
747    if(clip) RGB = fmax(RGB, 0.f);
748    RGB.w = pix_in.w;
749  }
750
751  write_imagef(out, (int2)(x, y), RGB);
752}
753