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

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