Convert from our home-grown matrices to Eigen all over.

diff --git a/colorspace_conversion_effect.cpp b/colorspace_conversion_effect.cpp
index ca1661b7868b95c7a7765fe8712f164dade2db9e..3f075d2b721141ef6c13428c24b93c854d4879f3 100644 (file)
--- a/colorspace_conversion_effect.cpp
+++ b/colorspace_conversion_effect.cpp
@@ -1,8 +1,12 @@
 #include <assert.h>
 
 #include <assert.h>
 

+#include <Eigen/LU>
+

 #include "colorspace_conversion_effect.h"
 #include "util.h"
 
 #include "colorspace_conversion_effect.h"
 #include "util.h"
 

+using namespace Eigen;
+

 // Color coordinates from Rec. 709; sRGB uses the same primaries.
 double rec709_x_R = 0.640,  rec709_x_G = 0.300,  rec709_x_B = 0.150;
 double rec709_y_R = 0.330,  rec709_y_G = 0.600,  rec709_y_B = 0.060;
 // Color coordinates from Rec. 709; sRGB uses the same primaries.
 double rec709_x_R = 0.640,  rec709_x_G = 0.300,  rec709_x_B = 0.150;
 double rec709_y_R = 0.330,  rec709_y_G = 0.600,  rec709_y_B = 0.060;


@@ -24,13 +28,10 @@ ColorspaceConversionEffect::ColorspaceConversionEffect()
        register_int("destination_space", (int *)&destination_space);
 }
 
        register_int("destination_space", (int *)&destination_space);
 }
 

-void get_xyz_matrix(Colorspace space, Matrix3x3 m)
+Matrix3d get_xyz_matrix(Colorspace space)

 {
        if (space == COLORSPACE_XYZ) {
 {
        if (space == COLORSPACE_XYZ) {

-               m[0] = 1.0f; m[3] = 0.0f; m[6] = 0.0f;
-               m[1] = 0.0f; m[4] = 1.0f; m[7] = 0.0f;
-               m[2] = 0.0f; m[5] = 0.0f; m[8] = 1.0f;
-               return;
+               return Matrix3d::Identity();

        }
 
        double x_R, x_G, x_B;
        }
 
        double x_R, x_G, x_B;


@@ -60,9 +61,11 @@ void get_xyz_matrix(Colorspace space, Matrix3x3 m)
 
        // Find the XYZ coordinates of D65 (white point for both Rec. 601 and 709),
        // normalized so that Y=1.
 
        // Find the XYZ coordinates of D65 (white point for both Rec. 601 and 709),
        // normalized so that Y=1.

-       double d65_X = d65_x / d65_y;
-       double d65_Y = 1.0;
-       double d65_Z = (1.0 - d65_x - d65_y) / d65_y;
+       Vector3d d65_XYZ(
+               d65_x / d65_y,
+               1.0,
+               (1.0 - d65_x - d65_y) / d65_y
+       );

 
        // We have, for each primary (example is with red):
        //
 
        // We have, for each primary (example is with red):
        //


@@ -91,34 +94,37 @@ void get_xyz_matrix(Colorspace space, Matrix3x3 m)
        //
        // Which we can solve for (Y_R, Y_G, Y_B) by inverting a 3x3 matrix.
 
        //
        // Which we can solve for (Y_R, Y_G, Y_B) by inverting a 3x3 matrix.
 

-       Matrix3x3 temp, inverted;
-       temp[0] = x_R / y_R;
-       temp[3] = x_G / y_G;
-       temp[6] = x_B / y_B;
+       Matrix3d temp;
+       temp(0,0) = x_R / y_R;
+       temp(0,1) = x_G / y_G;
+       temp(0,2) = x_B / y_B;

 
 

-       temp[1] = 1.0;
-       temp[4] = 1.0;
-       temp[7] = 1.0;
+       temp(1,0) = 1.0;
+       temp(1,1) = 1.0;
+       temp(1,2) = 1.0;

 
 

-       temp[2] = z_R / y_R;
-       temp[5] = z_G / y_G;
-       temp[8] = z_B / y_B;
+       temp(2,0) = z_R / y_R;
+       temp(2,1) = z_G / y_G;
+       temp(2,2) = z_B / y_B;

 
 

-       invert_3x3_matrix(temp, inverted);
-       float Y_R, Y_G, Y_B;
-       multiply_3x3_matrix_float3(inverted, d65_X, d65_Y, d65_Z, &Y_R, &Y_G, &Y_B);
+       Vector3d Y_RGB = temp.inverse() * d65_XYZ;

 
        // Now convert xyY -> XYZ.
 
        // Now convert xyY -> XYZ.

-       double X_R = temp[0] * Y_R;
-       double Z_R = temp[2] * Y_R;
-       double X_G = temp[3] * Y_G;
-       double Z_G = temp[5] * Y_G;
-       double X_B = temp[6] * Y_B;
-       double Z_B = temp[8] * Y_B;
-
-       m[0] = X_R; m[3] = X_G; m[6] = X_B;
-       m[1] = Y_R; m[4] = Y_G; m[7] = Y_B;
-       m[2] = Z_R; m[5] = Z_G; m[8] = Z_B;
+       double X_R = temp(0,0) * Y_RGB[0];
+       double Z_R = temp(2,0) * Y_RGB[0];
+
+       double X_G = temp(0,1) * Y_RGB[1];
+       double Z_G = temp(2,1) * Y_RGB[1];
+
+       double X_B = temp(0,2) * Y_RGB[2];
+       double Z_B = temp(2,2) * Y_RGB[2];
+
+       Matrix3d m;
+       m(0,0) = X_R;      m(0,1) = X_G;      m(0,2) = X_B;
+       m(1,0) = Y_RGB[0]; m(1,1) = Y_RGB[1]; m(1,2) = Y_RGB[2];
+       m(2,0) = Z_R;      m(2,1) = Z_G;      m(2,2) = Z_B;
+
+       return m;

 }
 
 std::string ColorspaceConversionEffect::output_fragment_shader()
 }
 
 std::string ColorspaceConversionEffect::output_fragment_shader()


@@ -129,26 +135,10 @@ std::string ColorspaceConversionEffect::output_fragment_shader()
        //
        // Since we right-multiply the RGB column vector, the matrix
        // concatenation order needs to be the opposite of the operation order.
        //
        // Since we right-multiply the RGB column vector, the matrix
        // concatenation order needs to be the opposite of the operation order.

-       Matrix3x3 m;
-
-       Matrix3x3 source_space_to_xyz;
-       Matrix3x3 destination_space_to_xyz;
-       Matrix3x3 xyz_to_destination_space;
-
-       get_xyz_matrix(source_space, source_space_to_xyz);
-       get_xyz_matrix(destination_space, destination_space_to_xyz);
-       invert_3x3_matrix(destination_space_to_xyz, xyz_to_destination_space);
-       
-       multiply_3x3_matrices(xyz_to_destination_space, source_space_to_xyz, m);
-
-       char buf[1024];
-       sprintf(buf,
-               "const mat3 PREFIX(conversion_matrix) = mat3(\n"
-               "    %.8f, %.8f, %.8f,\n"
-               "    %.8f, %.8f, %.8f,\n"
-               "    %.8f, %.8f, %.8f);\n\n",
-               m[0], m[1], m[2],
-               m[3], m[4], m[5],
-               m[6], m[7], m[8]);
-       return buf + read_file("colorspace_conversion_effect.frag");
+       Matrix3d source_space_to_xyz = get_xyz_matrix(source_space);
+       Matrix3d xyz_to_destination_space = get_xyz_matrix(destination_space).inverse();
+       Matrix3d m = xyz_to_destination_space * source_space_to_xyz;
+
+       return output_glsl_mat3("PREFIX(conversion_matrix)", m) +
+               read_file("colorspace_conversion_effect.frag");

 }
 }

diff --git a/effect.cpp b/effect.cpp
index 1847f4831f4792b73edc62fcf6f4df07c60406c3..9ff200bacde5d01258ffa2b7040181f266824efa 100644 (file)
--- a/effect.cpp
+++ b/effect.cpp
@@ -68,7 +68,7 @@ void set_uniform_vec4_array(GLuint glsl_program_num, const std::string &prefix,
        check_error();
 }
 
        check_error();
 }
 

-void set_uniform_mat3(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const Matrix3x3 matrix)
+void set_uniform_mat3(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const Eigen::Matrix3d& matrix)

 {
        GLint location = get_uniform_location(glsl_program_num, prefix, key);
        if (location == -1) {
 {
        GLint location = get_uniform_location(glsl_program_num, prefix, key);
        if (location == -1) {


@@ -78,8 +78,10 @@ void set_uniform_mat3(GLuint glsl_program_num, const std::string &prefix, const
 
        // Convert to float (GLSL has no double matrices).
        float matrixf[9];
 
        // Convert to float (GLSL has no double matrices).
        float matrixf[9];

-       for (unsigned i = 0; i < 9; ++i) {
-               matrixf[i] = matrix[i];
+       for (unsigned y = 0; y < 3; ++y) {
+               for (unsigned x = 0; x < 3; ++x) {
+                       matrixf[y + x * 3] = matrix(y, x);
+               }

        }
 
        glUniformMatrix3fv(location, 1, GL_FALSE, matrixf);
        }
 
        glUniformMatrix3fv(location, 1, GL_FALSE, matrixf);

diff --git a/effect.h b/effect.h
index 3a1fc9d9b4bd548defad4e96a0b9ace03dee08b6..db0fa94d43397804194021c910419363e826d14e 100644 (file)
--- a/effect.h
+++ b/effect.h
@@ -16,6 +16,8 @@
 
 #include <assert.h>
 
 
 #include <assert.h>
 

+#include <Eigen/Core>
+

 #include "opengl.h"
 #include "util.h"
 
 #include "opengl.h"
 #include "util.h"
 


@@ -45,7 +47,7 @@ void set_uniform_float(GLuint glsl_program_num, const std::string &prefix, const
 void set_uniform_vec2(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const float *values);
 void set_uniform_vec3(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const float *values);
 void set_uniform_vec4_array(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const float *values, size_t num_values);
 void set_uniform_vec2(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const float *values);
 void set_uniform_vec3(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const float *values);
 void set_uniform_vec4_array(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const float *values, size_t num_values);

-void set_uniform_mat3(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const Matrix3x3 matrix);
+void set_uniform_mat3(GLuint glsl_program_num, const std::string &prefix, const std::string &key, const Eigen::Matrix3d &matrix);

 
 class Effect {
 public:
 
 class Effect {
 public:

diff --git a/util.cpp b/util.cpp
index 9eaa262734f284b53aebba1d1c2672538dcc3aee..0d852e725ad65cb989ce4872d54635596d754c06 100644 (file)
--- a/util.cpp
+++ b/util.cpp
@@ -102,52 +102,25 @@ GLuint compile_shader(const std::string &shader_src, GLenum type)
        return obj;
 }
 
        return obj;
 }
 

-void multiply_3x3_matrices(const Matrix3x3 a, const Matrix3x3 b, Matrix3x3 result)
+void print_3x3_matrix(const Eigen::Matrix3d& m)

 {
 {

-        result[0] = a[0] * b[0] + a[3] * b[1] + a[6] * b[2];
-        result[1] = a[1] * b[0] + a[4] * b[1] + a[7] * b[2];
-        result[2] = a[2] * b[0] + a[5] * b[1] + a[8] * b[2];
-
-        result[3] = a[0] * b[3] + a[3] * b[4] + a[6] * b[5];
-        result[4] = a[1] * b[3] + a[4] * b[4] + a[7] * b[5];
-        result[5] = a[2] * b[3] + a[5] * b[4] + a[8] * b[5];
-
-        result[6] = a[0] * b[6] + a[3] * b[7] + a[6] * b[8];
-        result[7] = a[1] * b[6] + a[4] * b[7] + a[7] * b[8];
-        result[8] = a[2] * b[6] + a[5] * b[7] + a[8] * b[8];
-}
-
-void multiply_3x3_matrix_float3(const Matrix3x3 M, float x0, float x1, float x2, float *y0, float *y1, float *y2)
-{
-       *y0 = M[0] * x0 + M[3] * x1 + M[6] * x2;
-       *y1 = M[1] * x0 + M[4] * x1 + M[7] * x2;
-       *y2 = M[2] * x0 + M[5] * x1 + M[8] * x2;
-}
-
-void invert_3x3_matrix(const Matrix3x3 m, Matrix3x3 result)
-{
-       double inv_det = 1.0 / (
-               m[6] * m[1] * m[5] - m[6] * m[2] * m[4] -
-               m[3] * m[1] * m[8] + m[3] * m[2] * m[7] +
-               m[0] * m[4] * m[8] - m[0] * m[5] * m[7]);
-
-       result[0] = inv_det * (m[4] * m[8] - m[5] * m[7]);
-       result[1] = inv_det * (m[2] * m[7] - m[1] * m[8]);
-       result[2] = inv_det * (m[1] * m[5] - m[2] * m[4]);
-
-       result[3] = inv_det * (m[6] * m[5] - m[3] * m[8]);
-       result[4] = inv_det * (m[0] * m[8] - m[6] * m[2]);
-       result[5] = inv_det * (m[3] * m[2] - m[0] * m[5]);
-
-       result[6] = inv_det * (m[3] * m[7] - m[6] * m[4]);
-       result[7] = inv_det * (m[6] * m[1] - m[0] * m[7]);
-       result[8] = inv_det * (m[0] * m[4] - m[3] * m[1]);
+       printf("%6.4f %6.4f %6.4f\n", m(0,0), m(0,1), m(0,2));
+       printf("%6.4f %6.4f %6.4f\n", m(1,0), m(1,1), m(1,2));
+       printf("%6.4f %6.4f %6.4f\n", m(2,0), m(2,1), m(2,2));
+       printf("\n");

 }
 
 }
 

-void print_3x3_matrix(const Matrix3x3 m)
+std::string output_glsl_mat3(const std::string &name, const Eigen::Matrix3d &m)

 {
 {

-       printf("%6.4f %6.4f %6.4f\n", m[0], m[3], m[6]);
-       printf("%6.4f %6.4f %6.4f\n", m[1], m[4], m[7]);
-       printf("%6.4f %6.4f %6.4f\n", m[2], m[5], m[8]);
-       printf("\n");
+       char buf[1024];
+       sprintf(buf,
+               "const mat3 %s = mat3(\n"
+               "    %.8f, %.8f, %.8f,\n"
+               "    %.8f, %.8f, %.8f,\n"
+               "    %.8f, %.8f, %.8f);\n\n",
+               name.c_str(),
+               m(0,0), m(1,0), m(2,0),
+               m(0,1), m(1,1), m(2,1),
+               m(0,2), m(1,2), m(2,2));
+       return buf;

 }
 }

diff --git a/util.h b/util.h
index 4834bbdf4a86f4dc6048183f6f74d8b1259e3f10..b53bac76ce742367c618e47d1a4578d926bb8cd8 100644 (file)
--- a/util.h
+++ b/util.h
@@ -7,6 +7,7 @@
 #include <stdlib.h>
 
 #include <string>
 #include <stdlib.h>
 
 #include <string>

+#include <Eigen/Core>

 
 #include "opengl.h"
 
 
 #include "opengl.h"
 


@@ -19,9 +20,6 @@ void hsv2rgb(float h, float s, float v, float *r, float *g, float *b);
 // (ie. color luminance is as if S=0).
 void hsv2rgb_normalized(float h, float s, float v, float *r, float *g, float *b);
 
 // (ie. color luminance is as if S=0).
 void hsv2rgb_normalized(float h, float s, float v, float *r, float *g, float *b);
 

-// Column major (same as OpenGL).
-typedef double Matrix3x3[9];
-

 // Read a file from disk and return its contents.
 // Dies if the file does not exist.
 std::string read_file(const std::string &filename);
 // Read a file from disk and return its contents.
 // Dies if the file does not exist.
 std::string read_file(const std::string &filename);


@@ -30,17 +28,11 @@ std::string read_file(const std::string &filename);
 // and return the object number.
 GLuint compile_shader(const std::string &shader_src, GLenum type);
 
 // and return the object number.
 GLuint compile_shader(const std::string &shader_src, GLenum type);
 

-// Compute a * b.
-void multiply_3x3_matrices(const Matrix3x3 a, const Matrix3x3 b, Matrix3x3 result);
-
-// Compute M * [x0 x1 x2]'.
-void multiply_3x3_matrix_float3(const Matrix3x3 M, float x0, float x1, float x2, float *y0, float *y1, float *y2);
-
-// Compute m^-1. Result is undefined if the matrix is singular or near-singular.
-void invert_3x3_matrix(const Matrix3x3 m, Matrix3x3 result);
-

 // Print a 3x3 matrix to standard output. Useful for debugging.
 // Print a 3x3 matrix to standard output. Useful for debugging.

-void print_3x3_matrix(const Matrix3x3 m);
+void print_3x3_matrix(const Eigen::Matrix3d &m);
+
+// Output a GLSL 3x3 matrix declaration.
+std::string output_glsl_mat3(const std::string &name, const Eigen::Matrix3d &m);

 
 #ifdef NDEBUG
 #define check_error()
 
 #ifdef NDEBUG
 #define check_error()

diff --git a/white_balance_effect.cpp b/white_balance_effect.cpp
index 661212053b866ec03eb5c056ad7111d778b01530..d022984b71676a2abb5d38a2c1a1014a6860216f 100644 (file)
--- a/white_balance_effect.cpp
+++ b/white_balance_effect.cpp
@@ -1,42 +1,44 @@
 #include <math.h>
 #include <assert.h>
 
 #include <math.h>
 #include <assert.h>
 

+#include <Eigen/LU>
+

 #include "white_balance_effect.h"
 #include "util.h"
 #include "opengl.h"
 
 #include "white_balance_effect.h"
 #include "util.h"
 #include "opengl.h"
 

+using namespace Eigen;
+

 namespace {
 
 // Temperature is in Kelvin. Formula from http://en.wikipedia.org/wiki/Planckian_locus#Approximation .
 namespace {
 
 // Temperature is in Kelvin. Formula from http://en.wikipedia.org/wiki/Planckian_locus#Approximation .

-void convert_color_temperature_to_xyz(float T, float *x, float *y, float *z)
+Vector3d convert_color_temperature_to_xyz(float T)

 {
        double invT = 1.0 / T;
 {
        double invT = 1.0 / T;

-       double xc, yc;
+       double x, y;

 
        assert(T >= 1000.0f);
        assert(T <= 15000.0f);
 
        if (T <= 4000.0f) {
 
        assert(T >= 1000.0f);
        assert(T <= 15000.0f);
 
        if (T <= 4000.0f) {

-               xc = ((-0.2661239e9 * invT - 0.2343580e6) * invT + 0.8776956e3) * invT + 0.179910;
+               x = ((-0.2661239e9 * invT - 0.2343580e6) * invT + 0.8776956e3) * invT + 0.179910;

        } else {
        } else {

-               xc = ((-3.0258469e9 * invT + 2.1070379e6) * invT + 0.2226347e3) * invT + 0.240390;
+               x = ((-3.0258469e9 * invT + 2.1070379e6) * invT + 0.2226347e3) * invT + 0.240390;

        }
 
        if (T <= 2222.0f) {
        }
 
        if (T <= 2222.0f) {

-               yc = ((-1.1063814 * xc - 1.34811020) * xc + 2.18555832) * xc - 0.20219683;
+               y = ((-1.1063814 * x - 1.34811020) * x + 2.18555832) * x - 0.20219683;

        } else if (T <= 4000.0f) {
        } else if (T <= 4000.0f) {

-               yc = ((-0.9549476 * xc - 1.37418593) * xc + 2.09137015) * xc - 0.16748867;
+               y = ((-0.9549476 * x - 1.37418593) * x + 2.09137015) * x - 0.16748867;

        } else {
        } else {

-               yc = (( 3.0817580 * xc - 5.87338670) * xc + 3.75112997) * xc - 0.37001483;
+               y = (( 3.0817580 * x - 5.87338670) * x + 3.75112997) * x - 0.37001483;

        }
 
        }
 

-       *x = xc;
-       *y = yc;
-       *z = 1.0 - xc - yc;
+       return Vector3d(x, y, 1.0 - x - y);

 }
 
 // Assuming sRGB primaries, from Wikipedia.
 }
 
 // Assuming sRGB primaries, from Wikipedia.

-static const Matrix3x3 rgb_to_xyz_matrix = {
+double rgb_to_xyz_matrix[9] = {

        0.4124, 0.2126, 0.0193, 
        0.3576, 0.7152, 0.1192,
        0.1805, 0.0722, 0.9505,
        0.4124, 0.2126, 0.0193, 
        0.3576, 0.7152, 0.1192,
        0.1805, 0.0722, 0.9505,


@@ -54,7 +56,7 @@ static const Matrix3x3 rgb_to_xyz_matrix = {
  * in the first place) assumes the D65 illuminant, and so the D65 illuminant
  * also gives R=G=B in sRGB.
  */
  * in the first place) assumes the D65 illuminant, and so the D65 illuminant
  * also gives R=G=B in sRGB.
  */

-static const Matrix3x3 xyz_to_lms_matrix = {
+const double xyz_to_lms_matrix[9] = {

         0.4002, -0.2263,    0.0,
         0.7076,  1.1653,    0.0,
        -0.0808,  0.0457, 0.9182,
         0.4002, -0.2263,    0.0,
         0.7076,  1.1653,    0.0,
        -0.0808,  0.0457, 0.9182,


@@ -70,14 +72,18 @@ static const Matrix3x3 xyz_to_lms_matrix = {
  * for the reference color. Since L'=M'=S' and the Y row of the LMS-to-XYZ matrix
  * sums to unity, we know that Y'=L', and it's easy to find the fL that sets Y'=Y.
  */
  * for the reference color. Since L'=M'=S' and the Y row of the LMS-to-XYZ matrix
  * sums to unity, we know that Y'=L', and it's easy to find the fL that sets Y'=Y.
  */

-static void compute_lms_scaling_factors(float x, float y, float z, float *scale_l, float *scale_m, float *scale_s)
+Vector3d compute_lms_scaling_factors(const Vector3d &xyz)

 {
 {

-       float l, m, s;
-       multiply_3x3_matrix_float3(xyz_to_lms_matrix, x, y, z, &l, &m, &s);
+       Vector3d lms = Map<const Matrix3d>(xyz_to_lms_matrix) * xyz;
+       double l = lms[0];
+       double m = lms[1];
+       double s = lms[2];
+
+       double scale_l = xyz[1] / l;
+       double scale_m = scale_l * (l / m);
+       double scale_s = scale_l * (l / s);

 
 

-       *scale_l = y / l;
-       *scale_m = *scale_l * (l / m);
-       *scale_s = *scale_l * (l / s);
+       return Vector3d(scale_l, scale_m, scale_s);

 }
 
 }  // namespace
 }
 
 }  // namespace


@@ -97,14 +103,9 @@ std::string WhiteBalanceEffect::output_fragment_shader()
 
 void WhiteBalanceEffect::set_gl_state(GLuint glsl_program_num, const std::string &prefix, unsigned *sampler_num)
 {
 
 void WhiteBalanceEffect::set_gl_state(GLuint glsl_program_num, const std::string &prefix, unsigned *sampler_num)
 {

-       float x, y, z;
-       multiply_3x3_matrix_float3(rgb_to_xyz_matrix, neutral_color.r, neutral_color.g, neutral_color.b, &x, &y, &z);
-
-       float l, m, s;
-       multiply_3x3_matrix_float3(xyz_to_lms_matrix, x, y, z, &l, &m, &s);
-
-       float l_scale, m_scale, s_scale;
-       compute_lms_scaling_factors(x, y, z, &l_scale, &m_scale, &s_scale);
+       Vector3d rgb(neutral_color.r, neutral_color.g, neutral_color.b);
+       Vector3d xyz = Map<const Matrix3d>(rgb_to_xyz_matrix) * rgb;
+       Vector3d lms_scale = compute_lms_scaling_factors(xyz);

 
        /*
         * Now apply the color balance. Simply put, we find the chromacity point
 
        /*
         * Now apply the color balance. Simply put, we find the chromacity point


@@ -114,37 +115,26 @@ void WhiteBalanceEffect::set_gl_state(GLuint glsl_program_num, const std::string
         * since the D65 illuminant does not exactly match the results of T=6500K);
         * we normalize so that T=6500K really is a no-op.
         */
         * since the D65 illuminant does not exactly match the results of T=6500K);
         * we normalize so that T=6500K really is a no-op.
         */

-       float white_x, white_y, white_z, l_scale_white, m_scale_white, s_scale_white;
-       convert_color_temperature_to_xyz(output_color_temperature, &white_x, &white_y, &white_z);
-       compute_lms_scaling_factors(white_x, white_y, white_z, &l_scale_white, &m_scale_white, &s_scale_white);
-       
-       float ref_x, ref_y, ref_z, l_scale_ref, m_scale_ref, s_scale_ref;
-       convert_color_temperature_to_xyz(6500.0f, &ref_x, &ref_y, &ref_z);
-       compute_lms_scaling_factors(ref_x, ref_y, ref_z, &l_scale_ref, &m_scale_ref, &s_scale_ref);
-       
-       l_scale *= l_scale_ref / l_scale_white;
-       m_scale *= m_scale_ref / m_scale_white;
-       s_scale *= s_scale_ref / s_scale_white;
+       Vector3d white_xyz = convert_color_temperature_to_xyz(output_color_temperature);
+       Vector3d lms_scale_white = compute_lms_scaling_factors(white_xyz);
+
+       Vector3d ref_xyz = convert_color_temperature_to_xyz(6500.0f);
+       Vector3d lms_scale_ref = compute_lms_scaling_factors(ref_xyz);
+
+       lms_scale[0] *= lms_scale_ref[0] / lms_scale_white[0];
+       lms_scale[1] *= lms_scale_ref[1] / lms_scale_white[1];
+       lms_scale[2] *= lms_scale_ref[2] / lms_scale_white[2];

        
        /*
         * Concatenate all the different linear operations into a single 3x3 matrix.
         * Note that since we postmultiply our vectors, the order of the matrices
         * has to be the opposite of the execution order.
         */
        
        /*
         * Concatenate all the different linear operations into a single 3x3 matrix.
         * Note that since we postmultiply our vectors, the order of the matrices
         * has to be the opposite of the execution order.
         */

-       Matrix3x3 lms_to_xyz_matrix, xyz_to_rgb_matrix;
-       invert_3x3_matrix(xyz_to_lms_matrix, lms_to_xyz_matrix);
-       invert_3x3_matrix(rgb_to_xyz_matrix, xyz_to_rgb_matrix);
-
-       Matrix3x3 temp, temp2, corr_matrix;
-       Matrix3x3 lms_scale_matrix = {
-               l_scale,    0.0f,    0.0f,
-                  0.0f, m_scale,    0.0f,
-                  0.0f,    0.0f, s_scale,
-       };
-       multiply_3x3_matrices(xyz_to_rgb_matrix, lms_to_xyz_matrix, temp);
-       multiply_3x3_matrices(temp, lms_scale_matrix, temp2);
-       multiply_3x3_matrices(temp2, xyz_to_lms_matrix, temp);
-       multiply_3x3_matrices(temp, rgb_to_xyz_matrix, corr_matrix);
-
+       Matrix3d corr_matrix =
+               Map<const Matrix3d>(rgb_to_xyz_matrix).inverse() *
+               Map<const Matrix3d>(xyz_to_lms_matrix).inverse() *
+               lms_scale.asDiagonal() *
+               Map<const Matrix3d>(xyz_to_lms_matrix) *
+               Map<const Matrix3d>(rgb_to_xyz_matrix);

        set_uniform_mat3(glsl_program_num, prefix, "correction_matrix", corr_matrix);
 }
        set_uniform_mat3(glsl_program_num, prefix, "correction_matrix", corr_matrix);
 }

diff --git a/ycbcr_input.cpp b/ycbcr_input.cpp
index 8d58c7de8ab10c6d09adf68fa0bc8ef467f8be51..b1ad4ac7ad3249627fc5c7d309c0de8fc6bbe081 100644 (file)
--- a/ycbcr_input.cpp
+++ b/ycbcr_input.cpp
@@ -1,10 +1,14 @@
 #include <string.h>
 #include <assert.h>
 
 #include <string.h>
 #include <assert.h>
 

+#include <Eigen/LU>
+

 #include "ycbcr_input.h"
 #include "util.h"
 #include "opengl.h"
 
 #include "ycbcr_input.h"
 #include "util.h"
 #include "opengl.h"
 

+using namespace Eigen;
+

 YCbCrInput::YCbCrInput(const ImageFormat &image_format,
                        const YCbCrFormat &ycbcr_format,
                        unsigned width, unsigned height)
 YCbCrInput::YCbCrInput(const ImageFormat &image_format,
                        const YCbCrFormat &ycbcr_format,
                        unsigned width, unsigned height)


@@ -175,38 +179,29 @@ std::string YCbCrInput::output_fragment_shader()
        }
 
        // Matrix to convert RGB to YCbCr. See e.g. Rec. 601.
        }
 
        // Matrix to convert RGB to YCbCr. See e.g. Rec. 601.

-       Matrix3x3 rgb_to_ycbcr;
-       rgb_to_ycbcr[0] = coeff[0];
-       rgb_to_ycbcr[3] = coeff[1];
-       rgb_to_ycbcr[6] = coeff[2];
+       Matrix3d rgb_to_ycbcr;
+       rgb_to_ycbcr(0,0) = coeff[0];
+       rgb_to_ycbcr(0,1) = coeff[1];
+       rgb_to_ycbcr(0,2) = coeff[2];

 
        float cb_fac = (224.0 / 219.0) / (coeff[0] + coeff[1] + 1.0f - coeff[2]);
 
        float cb_fac = (224.0 / 219.0) / (coeff[0] + coeff[1] + 1.0f - coeff[2]);

-       rgb_to_ycbcr[1] = -coeff[0] * cb_fac;
-       rgb_to_ycbcr[4] = -coeff[1] * cb_fac;
-       rgb_to_ycbcr[7] = (1.0f - coeff[2]) * cb_fac;
+       rgb_to_ycbcr(1,0) = -coeff[0] * cb_fac;
+       rgb_to_ycbcr(1,1) = -coeff[1] * cb_fac;
+       rgb_to_ycbcr(1,2) = (1.0f - coeff[2]) * cb_fac;

 
        float cr_fac = (224.0 / 219.0) / (1.0f - coeff[0] + coeff[1] + coeff[2]);
 
        float cr_fac = (224.0 / 219.0) / (1.0f - coeff[0] + coeff[1] + coeff[2]);

-       rgb_to_ycbcr[2] = (1.0f - coeff[0]) * cr_fac;
-       rgb_to_ycbcr[5] = -coeff[1] * cr_fac;
-       rgb_to_ycbcr[8] = -coeff[2] * cr_fac;
+       rgb_to_ycbcr(2,0) = (1.0f - coeff[0]) * cr_fac;
+       rgb_to_ycbcr(2,1) = -coeff[1] * cr_fac;
+       rgb_to_ycbcr(2,2) = -coeff[2] * cr_fac;

 
        // Inverting the matrix gives us what we need to go from YCbCr back to RGB.
 
        // Inverting the matrix gives us what we need to go from YCbCr back to RGB.

-       Matrix3x3 ycbcr_to_rgb;
-       invert_3x3_matrix(rgb_to_ycbcr, ycbcr_to_rgb);
+       Matrix3d ycbcr_to_rgb = rgb_to_ycbcr.inverse();

 
        std::string frag_shader;
 
 
        std::string frag_shader;
 

-       char buf[1024];
-       sprintf(buf,
-               "const mat3 PREFIX(inv_ycbcr_matrix) = mat3(\n"
-               "    %.8f, %.8f, %.8f,\n"
-               "    %.8f, %.8f, %.8f,\n"
-               "    %.8f, %.8f, %.8f);\n",
-               ycbcr_to_rgb[0], ycbcr_to_rgb[1], ycbcr_to_rgb[2],
-               ycbcr_to_rgb[3], ycbcr_to_rgb[4], ycbcr_to_rgb[5],
-               ycbcr_to_rgb[6], ycbcr_to_rgb[7], ycbcr_to_rgb[8]);
-       frag_shader = buf;
+       frag_shader = output_glsl_mat3("PREFIX(inv_ycbcr_matrix)", ycbcr_to_rgb);

 
 

+       char buf[256];

        sprintf(buf, "const vec3 PREFIX(offset) = vec3(%.8f, %.8f, %.8f);\n",
                offset[0], offset[1], offset[2]);
        frag_shader += buf;
        sprintf(buf, "const vec3 PREFIX(offset) = vec3(%.8f, %.8f, %.8f);\n",
                offset[0], offset[1], offset[2]);
        frag_shader += buf;