From 7941becb1722913c82716bad98f64e2e87b49adb Mon Sep 17 00:00:00 2001
From: Martijn Braam <martijn@brixit.nl>
Date: Tue, 25 Jul 2023 13:28:42 +0200
Subject: [PATCH] Use more colorscience from millipixels

---
 src/gles2_debayer.c | 15 +++++++++------
 src/matrix.h        | 29 +++++++++++++++++++++++++++++
 2 files changed, 38 insertions(+), 6 deletions(-)

diff --git a/src/gles2_debayer.c b/src/gles2_debayer.c
index e02bfe5..55e5436 100644
--- a/src/gles2_debayer.c
+++ b/src/gles2_debayer.c
@@ -3,6 +3,7 @@
 #include "camera.h"
 #include "dcp.h"
 #include "gl_util.h"
+#include "matrix.h"
 #include <stdlib.h>
 
 #define VERTEX_ATTRIBUTE 0
@@ -121,6 +122,7 @@ gles2_debayer_configure(GLES2Debayer *self,
         glViewport(0, 0, (int)dst_width, (int)dst_height);
         check_gl();
 
+        /* Rotation matrix for orientation correction */
         GLfloat rotation_list[4] = { 0, -1, 0, 1 };
         int rotation_index = 4 - (int)rotation / 90;
 
@@ -142,8 +144,9 @@ gles2_debayer_configure(GLES2Debayer *self,
         glUniform2f(self->uniform_pixel_size, pixel_size_x, pixel_size_y);
         check_gl();
 
+        /* Color calibration curves and matrices */
         float gamma = 1.0f;
-        for (int i = 0; i < calibration.tone_curve_length * 2; i += 2) {
+        for (int i = 2; i < calibration.tone_curve_length; i += 2) {
                 float g = calibration.tone_curve[i + 1] / calibration.tone_curve[i];
                 if (g > gamma) {
                         gamma = g;
@@ -153,11 +156,11 @@ gles2_debayer_configure(GLES2Debayer *self,
 
         if (calibration.color_matrix_1[0]) {
                 GLfloat transposed[9];
-                for (int i = 0; i < 3; ++i)
-                        for (int j = 0; j < 3; ++j)
-                                transposed[i + j * 3] =
-                                        calibration.color_matrix_1[j + i * 3];
-
+                float colormat_inv[9];
+                float colormat[9];
+                invert_matrix(calibration.color_matrix_1, colormat_inv);
+                multiply_matrices(xyz_to_srgb, colormat_inv, colormat);
+                transpose_matrix(colormat, transposed);
                 glUniformMatrix3fv(
                         self->uniform_color_matrix, 1, GL_FALSE, transposed);
         } else {
diff --git a/src/matrix.h b/src/matrix.h
index 19a705d..3ab68c7 100644
--- a/src/matrix.h
+++ b/src/matrix.h
@@ -1 +1,30 @@
+static float xyz_to_srgb[] = { 3.2404542f,  -1.5371385f, -0.4985314f,
+                               -0.9692660f, 1.8760108f,  0.0415560f,
+                               0.0556434f,  -0.2040259f, 1.0572252f };
+
 void multiply_matrices(float a[9], float b[9], float out[9]);
+
+void
+invert_matrix(const float in[9], float out[9])
+{
+        float det = in[0] * (in[4] * in[8] - in[5] * in[7]) -
+                    in[1] * (in[3] * in[8] - in[5] * in[6]) +
+                    in[2] * (in[3] * in[7] - in[4] * in[7]);
+        out[0] = (in[4] * in[8] - in[7] * in[5]) / det;
+        out[1] = (in[7] * in[2] - in[1] * in[8]) / det;
+        out[2] = (in[1] * in[5] - in[4] * in[2]) / det;
+        out[3] = (in[6] * in[5] - in[3] * in[8]) / det;
+        out[4] = (in[0] * in[8] - in[6] * in[5]) / det;
+        out[5] = (in[3] * in[2] - in[0] * in[5]) / det;
+        out[6] = (in[3] * in[7] - in[6] * in[4]) / det;
+        out[7] = (in[6] * in[1] - in[0] * in[7]) / det;
+        out[8] = (in[0] * in[4] - in[3] * in[1]) / det;
+}
+
+void
+transpose_matrix(const float in[9], float out[9])
+{
+        for (int i = 0; i < 3; ++i)
+                for (int j = 0; j < 3; ++j)
+                        out[i + j * 3] = in[j + i * 3];
+}
\ No newline at end of file