git.sesse.net Git - movit/blob - white_balance_effect.cpp

   1 #include <math.h>
   2 #include <assert.h>
   3
   4 #include <Eigen/LU>
   5
   6 #include "white_balance_effect.h"
   7 #include "util.h"
   8 #include "opengl.h"
   9
  10 using namespace Eigen;
  11
  12 namespace {
  13
  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;
  19
  20         assert(T >= 1000.0f);
  21         assert(T <= 15000.0f);
  22
  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         }
  28
  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         }
  36
  37         return Vector3d(x, y, 1.0 - x - y);
  38 }
  39
  40 // Assuming sRGB primaries, from Wikipedia.
  41 double rgb_to_xyz_matrix[9] = {
  42         0.4124, 0.2126, 0.0193,
  43         0.3576, 0.7152, 0.1192,
  44         0.1805, 0.0722, 0.9505,
  45 };
  46
  47 /*
  48  * There are several different LMS spaces, at least according to Wikipedia.
  49  * Through practical testing, I've found most of them (like the CIECAM02 model)
  50  * to yield a result that is too reddish in practice, possibly because they
  51  * are intended for different illuminants than what sRGB assumes.
  52  *
  53  * This is the RLAB space, normalized to D65, which means that the standard
  54  * D65 illuminant (x=0.31271, y=0.32902, z=1-y-x) gives L=M=S under this transformation.
  55  * This makes sense because sRGB (which is used to derive those XYZ values
  56  * in the first place) assumes the D65 illuminant, and so the D65 illuminant
  57  * also gives R=G=B in sRGB.
  58  */
  59 const double xyz_to_lms_matrix[9] = {
  60          0.4002, -0.2263,    0.0,
  61          0.7076,  1.1653,    0.0,
  62         -0.0808,  0.0457, 0.9182,
  63 };
  64
  65 /*
  66  * For a given reference color (given in XYZ space),
  67  * compute scaling factors for L, M and S. What we want at the output is equal L, M and S
  68  * for the reference color (making it a neutral illuminant), or sL ref_L = sM ref_M = sS ref_S.
  69  * This removes two degrees of freedom for our system, and we only need to find fL.
  70  *
  71  * A reasonable last constraint would be to preserve Y, approximately the brightness,
  72  * for the reference color. Since L'=M'=S' and the Y row of the LMS-to-XYZ matrix
  73  * sums to unity, we know that Y'=L', and it's easy to find the fL that sets Y'=Y.
  74  */
  75 Vector3d compute_lms_scaling_factors(const Vector3d &xyz)
  76 {
  77         Vector3d lms = Map<const Matrix3d>(xyz_to_lms_matrix) * xyz;
  78         double l = lms[0];
  79         double m = lms[1];
  80         double s = lms[2];
  81
  82         double scale_l = xyz[1] / l;
  83         double scale_m = scale_l * (l / m);
  84         double scale_s = scale_l * (l / s);
  85
  86         return Vector3d(scale_l, scale_m, scale_s);
  87 }
  88
  89 }  // namespace
  90
  91 WhiteBalanceEffect::WhiteBalanceEffect()
  92         : neutral_color(0.5f, 0.5f, 0.5f),
  93           output_color_temperature(6500.0f)
  94 {
  95         register_vec3("neutral_color", (float *)&neutral_color);
  96         register_float("output_color_temperature", &output_color_temperature);
  97 }
  98
  99 std::string WhiteBalanceEffect::output_fragment_shader()
 100 {
 101         return read_file("white_balance_effect.frag");
 102 }
 103
 104 void WhiteBalanceEffect::set_gl_state(GLuint glsl_program_num, const std::string &prefix, unsigned *sampler_num)
 105 {
 106         Vector3d rgb(neutral_color.r, neutral_color.g, neutral_color.b);
 107         Vector3d xyz = Map<const Matrix3d>(rgb_to_xyz_matrix) * rgb;
 108         Vector3d lms_scale = compute_lms_scaling_factors(xyz);
 109
 110         /*
 111          * Now apply the color balance. Simply put, we find the chromacity point
 112          * for the desired white temperature, see what LMS scaling factors they
 113          * would have given us, and then reverse that transform. For T=6500K,
 114          * the default, this gives us nearly an identity transform (but only nearly,
 115          * since the D65 illuminant does not exactly match the results of T=6500K);
 116          * we normalize so that T=6500K really is a no-op.
 117          */
 118         Vector3d white_xyz = convert_color_temperature_to_xyz(output_color_temperature);
 119         Vector3d lms_scale_white = compute_lms_scaling_factors(white_xyz);
 120
 121         Vector3d ref_xyz = convert_color_temperature_to_xyz(6500.0f);
 122         Vector3d lms_scale_ref = compute_lms_scaling_factors(ref_xyz);
 123
 124         lms_scale[0] *= lms_scale_ref[0] / lms_scale_white[0];
 125         lms_scale[1] *= lms_scale_ref[1] / lms_scale_white[1];
 126         lms_scale[2] *= lms_scale_ref[2] / lms_scale_white[2];
 127
 128         /*
 129          * Concatenate all the different linear operations into a single 3x3 matrix.
 130          * Note that since we postmultiply our vectors, the order of the matrices
 131          * has to be the opposite of the execution order.
 132          */
 133         Matrix3d corr_matrix =
 134                 Map<const Matrix3d>(rgb_to_xyz_matrix).inverse() *
 135                 Map<const Matrix3d>(xyz_to_lms_matrix).inverse() *
 136                 lms_scale.asDiagonal() *
 137                 Map<const Matrix3d>(xyz_to_lms_matrix) *
 138                 Map<const Matrix3d>(rgb_to_xyz_matrix);
 139         set_uniform_mat3(glsl_program_num, prefix, "correction_matrix", corr_matrix);
 140 }