X-Git-Url: https://git.sesse.net/?p=movit;a=blobdiff_plain;f=white_balance_effect.cpp;h=69759adfc6148379698998674109a5ae17213e4b;hp=661212053b866ec03eb5c056ad7111d778b01530;hb=95edbfccb0843da3cc105dadc5bc6d8e102f6071;hpb=62c039a63a8afdea4ec2da081d0173d7dd1d4578 diff --git a/white_balance_effect.cpp b/white_balance_effect.cpp index 6612120..69759ad 100644 --- a/white_balance_effect.cpp +++ b/white_balance_effect.cpp @@ -1,83 +1,100 @@ -#include +#include +#include +#include #include -#include "white_balance_effect.h" +#include "colorspace_conversion_effect.h" +#include "d65.h" +#include "effect_util.h" +#include "image_format.h" #include "util.h" -#include "opengl.h" +#include "white_balance_effect.h" + +using namespace Eigen; +using namespace std; + +namespace movit { 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 xc, yc; + double x, y; 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.2343589e6) * invT + 0.8776956e3) * invT + 0.179910; } 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) { - 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) { - 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 { - 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. -static const Matrix3x3 rgb_to_xyz_matrix = { - 0.4124, 0.2126, 0.0193, - 0.3576, 0.7152, 0.1192, - 0.1805, 0.0722, 0.9505, -}; - /* - * There are several different LMS spaces, at least according to Wikipedia. - * Through practical testing, I've found most of them (like the CIECAM02 model) - * to yield a result that is too reddish in practice, possibly because they - * are intended for different illuminants than what sRGB assumes. + * There are several different perceptual color spaces with different intended + * uses; for instance, CIECAM02 uses one space (CAT02) for purposes of computing + * chromatic adaptation (the effect that the human eye perceives an object as + * the same color even under differing illuminants), but a different space + * (Hunt-Pointer-Estevez, or HPE) for the actual perception post-adaptation. * - * This is the RLAB space, normalized to D65, which means that the standard - * D65 illuminant (x=0.31271, y=0.32902, z=1-y-x) gives L=M=S under this transformation. - * This makes sense because sRGB (which is used to derive those XYZ values - * in the first place) assumes the D65 illuminant, and so the D65 illuminant - * also gives R=G=B in sRGB. + * CIECAM02 chromatic adaptation, while related to the transformation we want, + * is a more complex phenomenon that depends on factors like the viewing conditions + * (e.g. amount of surrounding light), and can no longer be implemented by just scaling + * each component in LMS space. The simpler way out is to use the HPE matrix, + * which is intended to be close to the actual cone response; this results in + * the “von Kries transformation” when we couple it with normalization in LMS space. + * + * http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html compares + * von Kries transformation with using another matrix, the Bradford matrix, + * and generally finds that the Bradford method gives a better result, + * as in giving better matches with the true result (as calculated using + * spectral matching) when converting between various CIE illuminants. + * The actual perceptual differences were found to be minor, though. + * We use the Bradford tranformation matrix from that page, and compute the + * inverse ourselves. (The Bradford matrix is also used in CMCCAT97.) */ -static const Matrix3x3 xyz_to_lms_matrix = { - 0.4002, -0.2263, 0.0, - 0.7076, 1.1653, 0.0, - -0.0808, 0.0457, 0.9182, +const double xyz_to_lms_matrix[9] = { + 0.7328, -0.7036, 0.0030, + 0.4296, 1.6975, 0.0136, + -0.1624, 0.0061, 0.9834, }; /* - * For a given reference color (given in XYZ space), - * compute scaling factors for L, M and S. What we want at the output is equal L, M and S - * for the reference color (making it a neutral illuminant), or sL ref_L = sM ref_M = sS ref_S. - * This removes two degrees of freedom for our system, and we only need to find fL. + * For a given reference color (given in XYZ space), compute scaling factors + * for L, M and S. What we want at the output is turning the reference color + * into a scaled version of the D65 illuminant (giving it R=G=B in sRGB), or + * + * (sL ref_L, sM ref_M, sS ref_S) = (s d65_L, s d65_M, s d65_S) * + * This removes two degrees of freedom from our system, and we only need to find s. * A reasonable last constraint would be to preserve Y, approximately the brightness, - * 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. Thus, we choose our D65 illuminant's Y such that it is + * equal to the reference color's Y, and the rest is easy. */ -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 &ref_xyz) { - float l, m, s; - multiply_3x3_matrix_float3(xyz_to_lms_matrix, x, y, z, &l, &m, &s); + Vector3d ref_lms = Map(xyz_to_lms_matrix) * ref_xyz; + Vector3d d65_lms = Map(xyz_to_lms_matrix) * + (ref_xyz[1] * Vector3d(d65_X, d65_Y, d65_Z)); // d65_Y = 1.0. + + double scale_l = d65_lms[0] / ref_lms[0]; + double scale_m = d65_lms[1] / ref_lms[1]; + double scale_s = d65_lms[2] / ref_lms[2]; - *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 @@ -90,21 +107,17 @@ WhiteBalanceEffect::WhiteBalanceEffect() register_float("output_color_temperature", &output_color_temperature); } -std::string WhiteBalanceEffect::output_fragment_shader() +string WhiteBalanceEffect::output_fragment_shader() { return read_file("white_balance_effect.frag"); } -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 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); + Matrix3d rgb_to_xyz_matrix = ColorspaceConversionEffect::get_xyz_matrix(COLORSPACE_sRGB); + Vector3d rgb(neutral_color.r, neutral_color.g, neutral_color.b); + Vector3d xyz = 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 @@ -114,37 +127,28 @@ 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. */ - 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. */ - 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 = + rgb_to_xyz_matrix.inverse() * + Map(xyz_to_lms_matrix).inverse() * + lms_scale.asDiagonal() * + Map(xyz_to_lms_matrix) * + rgb_to_xyz_matrix; set_uniform_mat3(glsl_program_num, prefix, "correction_matrix", corr_matrix); } + +} // namespace movit