1 files changed, 49 insertions, 0 deletions
diff --git a/libglitter_per_channel_desaturate_double.c b/libglitter_per_channel_desaturate_double.c
new file mode 100644
index 0000000..8561708
--- /dev/null
+++ b/libglitter_per_channel_desaturate_double.c
@@ -0,0 +1,49 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic ignored "-Wunsuffixed-float-constants"
+#endif
+
+
+void
+libglitter_per_channel_desaturate_double(double **rasters, size_t nrasters, size_t rowsize,
+                                         size_t cellsize, size_t width, size_t height,
+                                         const double *saturations, const double *primary_ys)
+{
+	size_t y, x, ch, raster_y, raster_i;
+	double intensity;
+
+	rowsize *= cellsize;
+
+	if (nrasters == 3) {
+		double *raster1 = rasters[0], saturation1 = saturations[0];
+		double *raster2 = rasters[1], saturation2 = saturations[1];
+		double *raster3 = rasters[2], saturation3 = saturations[2];
+		double primary1_y = primary_ys ? primary_ys[0] : (double)0.212673370378408277403536885686;
+		double primary2_y = primary_ys ? primary_ys[1] : (double)0.715151730491031756287156895269;
+		double primary3_y = primary_ys ? primary_ys[2] : (double)0.072174899130559869164791564344;
+		for (y = 0, raster_y = 0; y < height; y++, raster_y += rowsize) {
+			for (x = 0, raster_i = raster_y; x < width; x++, raster_i += cellsize) {
+				intensity  = raster1[raster_i] * primary1_y;
+				intensity += raster2[raster_i] * primary2_y;
+				intensity += raster3[raster_i] * primary3_y;
+				raster1[raster_i] = fma(raster1[raster_i] - intensity, saturation1, intensity);
+				raster2[raster_i] = fma(raster2[raster_i] - intensity, saturation2, intensity);
+				raster3[raster_i] = fma(raster3[raster_i] - intensity, saturation3, intensity);
+			}
+		}
+
+	} else {
+		for (y = 0, raster_y = 0; y < height; y++, raster_y += rowsize) {
+			for (x = 0, raster_i = raster_y; x < width; x++, raster_i += cellsize) {
+				intensity = 0;
+				for (ch = 0; ch < nrasters; ch++)
+					intensity += rasters[ch][raster_i] * primary_ys[ch];
+				for (ch = 0; ch < nrasters; ch++)
+					rasters[ch][raster_i] = fma(rasters[ch][raster_i] - intensity, saturations[ch], intensity);
+			}
+		}
+	}
+}