1583922461053dd1c793eaa51f93b28275b1cab5
1 #include <Eigen/Core>
2 #include <Eigen/LU>
3 #include <GL/glew.h>
4 #include <assert.h>
6 #include "d65.h"
7 #include "util.h"
8 #include "white_balance_effect.h"
10 using namespace Eigen;
12 namespace {
14 // Temperature is in Kelvin. Formula from http://en.wikipedia.org/wiki/Planckian_locus#Approximation .
15 Vector3d convert_color_temperature_to_xyz(float T)
16 {
17         double invT = 1.0 / T;
18         double x, y;
20         assert(T >= 1000.0f);
21         assert(T <= 15000.0f);
23         if (T <= 4000.0f) {
24                 x = ((-0.2661239e9 * invT - 0.2343580e6) * invT + 0.8776956e3) * invT + 0.179910;
25         } else {
26                 x = ((-3.0258469e9 * invT + 2.1070379e6) * invT + 0.2226347e3) * invT + 0.240390;
27         }
29         if (T <= 2222.0f) {
30                 y = ((-1.1063814 * x - 1.34811020) * x + 2.18555832) * x - 0.20219683;
31         } else if (T <= 4000.0f) {
32                 y = ((-0.9549476 * x - 1.37418593) * x + 2.09137015) * x - 0.16748867;
33         } else {
34                 y = (( 3.0817580 * x - 5.87338670) * x + 3.75112997) * x - 0.37001483;
35         }
37         return Vector3d(x, y, 1.0 - x - y);
38 }
40 // Assuming sRGB primaries, from Wikipedia.
41 const double rgb_to_xyz_matrix = {
42         0.4124, 0.2126, 0.0193,
43         0.3576, 0.7152, 0.1192,
44         0.1805, 0.0722, 0.9505,
45 };
47 /*
48  * There are several different perceptual color spaces with different intended
49  * uses; for instance, CIECAM02 uses one space (CAT02) for purposes of computing
50  * chromatic adaptation (the effect that the human eye perceives an object as
51  * the same color even under differing illuminants), but a different space
52  * (Hunt-Pointer-Estevez, or HPE) for the actual perception post-adaptation.
53  *
54  * CIECAM02 chromatic adaptation, while related to the transformation we want,
55  * is a more complex phenomenon that depends on factors like the viewing conditions
56  * (e.g. amount of surrounding light), and can no longer be implemented by just scaling
57  * each component in LMS space. The simpler way out is to use the HPE matrix,
58  * which is intended to be close to the actual cone response; this results in
59  * the “von Kries transformation” when we couple it with normalization in LMS space.
60  *
62  * von Kries transformation with using another matrix, the Bradford matrix,
63  * and generally finds that the Bradford method gives a better result,
64  * as in giving better matches with the true result (as calculated using
65  * spectral matching) when converting between various CIE illuminants.
66  * The actual perceptual differences were found to be minor, though.
67  * We use the Bradford tranformation matrix from that page, and compute the
68  * inverse ourselves. (The Bradford matrix is also used in CMCCAT97.)
69  */
70 const double xyz_to_lms_matrix = {
71          0.7328, -0.7036,  0.0030,
72          0.4296,  1.6975,  0.0136,
73         -0.1624,  0.0061,  0.9834,
74 };
76 /*
77  * For a given reference color (given in XYZ space), compute scaling factors
78  * for L, M and S. What we want at the output is turning the reference color
79  * into a scaled version of the D65 illuminant (giving it R=G=B in sRGB), or
80  *
81  *   (sL ref_L, sM ref_M, sS ref_S) = (s d65_L, s d65_M, s d65_S)
82  *
83  * This removes two degrees of freedom from our system, and we only need to find s.
84  * A reasonable last constraint would be to preserve Y, approximately the brightness,
85  * for the reference color. Thus, we choose our D65 illuminant's Y such that it is
86  * equal to the reference color's Y, and the rest is easy.
87  */
88 Vector3d compute_lms_scaling_factors(const Vector3d &ref_xyz)
89 {
90         Vector3d ref_lms = Map<const Matrix3d>(xyz_to_lms_matrix) * ref_xyz;
91         Vector3d d65_lms = Map<const Matrix3d>(xyz_to_lms_matrix) *
92                 (ref_xyz * Vector3d(d65_X, d65_Y, d65_Z));  // d65_Y = 1.0.
94         double scale_l = d65_lms / ref_lms;
95         double scale_m = d65_lms / ref_lms;
96         double scale_s = d65_lms / ref_lms;
98         return Vector3d(scale_l, scale_m, scale_s);
99 }
101 }  // namespace
103 WhiteBalanceEffect::WhiteBalanceEffect()
104         : neutral_color(0.5f, 0.5f, 0.5f),
105           output_color_temperature(6500.0f)
106 {
107         register_vec3("neutral_color", (float *)&neutral_color);
108         register_float("output_color_temperature", &output_color_temperature);
109 }
112 {
114 }
116 void WhiteBalanceEffect::set_gl_state(GLuint glsl_program_num, const std::string &prefix, unsigned *sampler_num)
117 {
118         Vector3d rgb(neutral_color.r, neutral_color.g, neutral_color.b);
119         Vector3d xyz = Map<const Matrix3d>(rgb_to_xyz_matrix) * rgb;
120         Vector3d lms_scale = compute_lms_scaling_factors(xyz);
122         /*
123          * Now apply the color balance. Simply put, we find the chromacity point
124          * for the desired white temperature, see what LMS scaling factors they
125          * would have given us, and then reverse that transform. For T=6500K,
126          * the default, this gives us nearly an identity transform (but only nearly,
127          * since the D65 illuminant does not exactly match the results of T=6500K);
128          * we normalize so that T=6500K really is a no-op.
129          */
130         Vector3d white_xyz = convert_color_temperature_to_xyz(output_color_temperature);
131         Vector3d lms_scale_white = compute_lms_scaling_factors(white_xyz);
133         Vector3d ref_xyz = convert_color_temperature_to_xyz(6500.0f);
134         Vector3d lms_scale_ref = compute_lms_scaling_factors(ref_xyz);
136         lms_scale *= lms_scale_ref / lms_scale_white;
137         lms_scale *= lms_scale_ref / lms_scale_white;
138         lms_scale *= lms_scale_ref / lms_scale_white;
140         /*
141          * Concatenate all the different linear operations into a single 3x3 matrix.
142          * Note that since we postmultiply our vectors, the order of the matrices
143          * has to be the opposite of the execution order.
144          */
145         Matrix3d corr_matrix =
146                 Map<const Matrix3d>(rgb_to_xyz_matrix).inverse() *
147                 Map<const Matrix3d>(xyz_to_lms_matrix).inverse() *
148                 lms_scale.asDiagonal() *
149                 Map<const Matrix3d>(xyz_to_lms_matrix) *
150                 Map<const Matrix3d>(rgb_to_xyz_matrix);
151         set_uniform_mat3(glsl_program_num, prefix, "correction_matrix", corr_matrix);
152 }