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