Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

bugfix : Fix filmlight RGB / color balance RGB #8976

Merged
merged 6 commits into from May 15, 2021
Merged

bugfix : Fix filmlight RGB / color balance RGB #8976

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
1fe5fcf
colorbalancergb: gamut-map only at input and in JzAzBz
aurelienpierre May 13, 2021
7def6e2
remove redundant y == 0 checks for atan2(x, y)
aurelienpierre May 13, 2021
9a4bea5
chromatic_adaptation.h : refactor functions
aurelienpierre May 15, 2021
74e2766
libs refactor : move dot product (vector and matrix) to math.h
aurelienpierre May 15, 2021
933351e
colorbalancergb: fix color spaces conversions
aurelienpierre May 15, 2021
0149e6f
chromatic_adaptation: change variable name to please Win
aurelienpierre May 15, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
161 changes: 117 additions & 44 deletions data/kernels/colorspace.h
View file
Open in desktop
Expand Up @@ -17,6 +17,17 @@
along with darktable. If not, see <http://www.gnu.org/licenses/>.
*/

inline float4 matrix_dot(const float4 vector, const float4 matrix[3])
{
float4 output;
const float4 vector_copy = { vector.x, vector.y, vector.z, 0.f };
output.x = dot(vector_copy, matrix[0]);
output.y = dot(vector_copy, matrix[1]);
output.z = dot(vector_copy, matrix[2]);
output.w = vector.w;
return output;
}

inline float4 Lab_2_LCH(float4 Lab)
{
float H = atan2(Lab.z, Lab.y);
Expand Down Expand Up @@ -416,65 +427,127 @@ inline float4 JzAzBz_to_JzCzhz(float4 JzAzBz)
}


inline float4 gradingRGB_to_Ych(float4 RGB)
// Convert CIE 1931 2° XYZ D65 to CIE 2006 LMS D65 (cone space)
/*
* The CIE 1931 XYZ 2° observer D65 is converted to CIE 2006 LMS D65 using the approximation by
* Richard A. Kirk, Chromaticity coordinates for graphic arts based on CIE 2006 LMS
* with even spacing of Munsell colours
* https://doi.org/10.2352/issn.2169-2629.2019.27.38
*/

inline float4 XYZ_to_LMS(const float4 XYZ)
{
const float4 D65 = { 0.18600766f, 0.5908061f, 0.22318624f, 0.f };
float4 Ych;
Ych.x = fmax(0.67282368f * RGB.x + 0.47812261f * RGB.y + 0.01044966f * RGB.z, 0.f);
const float a = RGB.x + RGB.y + RGB.z;
RGB = (a == 0.f) ? (float4)0.f : RGB / a;
const float4 XYZ_D65_to_LMS_2006_D65[3]
= { { 0.257085f, 0.859943f, -0.031061f, 0.f },
{ -0.394427f, 1.175800f, 0.106423f, 0.f },
{ 0.064856f, -0.076250f, 0.559067f, 0.f } };

const float *pD65 = &D65;
return matrix_dot(XYZ, XYZ_D65_to_LMS_2006_D65);
}

RGB.x -= pD65[0];
RGB.y -= pD65[1];

Ych.y = hypot(RGB.y, RGB.x);
Ych.z = (Ych.x == 0.f) ? 0.f : atan2(RGB.y, RGB.x);
Ych.w = RGB.w;
return Ych;
inline float4 LMS_to_XYZ(const float4 LMS)
{
const float4 LMS_2006_D65_to_XYZ_D65[3]
= { { 1.80794659f, -1.29971660f, 0.34785879f, 0.f },
{ 0.61783960f, 0.39595453f, -0.04104687f, 0.f },
{ -0.12546960f, 0.20478038f, 1.74274183f, 0.f } };

return matrix_dot(LMS, LMS_2006_D65_to_XYZ_D65);
}


inline float4 Ych_to_gradingRGB(const float4 Ych)
/*
* Convert from CIE 2006 LMS D65 to Filmlight RGB defined in
* Richard A. Kirk, Chromaticity coordinates for graphic arts based on CIE 2006 LMS
* with even spacing of Munsell colours
* https://doi.org/10.2352/issn.2169-2629.2019.27.38
*/

inline float4 gradingRGB_to_LMS(const float4 RGB)
{
const float4 D65 = { 0.18600766f, 0.5908061f, 0.22318624f, 0.f };
const float4 filmlightRGB_D65_to_LMS_D65[3]
= { { 0.95f, 0.38f, 0.00f, 0.f },
{ 0.05f, 0.62f, 0.03f, 0.f },
{ 0.00f, 0.00f, 0.97f, 0.f } };

const float *pD65 = &D65;
float4 RGB;
RGB.x = Ych.y * native_cos(Ych.z) + pD65[0];
RGB.y = Ych.y * native_sin(Ych.z) + pD65[1];
RGB.z = 1.f - RGB.x - RGB.y;
return matrix_dot(RGB, filmlightRGB_D65_to_LMS_D65);
}

const float a = (0.67282368f * RGB.x + 0.47812261f * RGB.y + 0.01044966f * RGB.z);
RGB = (a == 0.f) ? (float4)0.f : RGB * Ych.x / a;
RGB.w = Ych.w;
return RGB;
inline float4 LMS_to_gradingRGB(const float4 LMS)
{
const float4 LMS_D65_to_filmlightRGB_D65[3]
= { { 1.0877193f, -0.66666667f, 0.02061856f, 0.f },
{ -0.0877193f, 1.66666667f, -0.05154639f, 0.f },
{ 0.f, 0.f, 1.03092784f, 0.f } };

return matrix_dot(LMS, LMS_D65_to_filmlightRGB_D65);
}


/*
* Re-express the Filmlight RGB triplet as Yrg luminance/chromacity coordinates
*/

inline float4 LMS_to_Yrg(const float4 LMS)
{
// compute luminance
const float Y = 0.68990272f * LMS.x + 0.34832189f * LMS.y;

// normalize LMS
const float a = LMS.x + LMS.y + LMS.z;
const float4 lms = (a == 0.f) ? 0.f : LMS / a;

// convert to Filmlight rgb (normalized)
const float4 rgb = LMS_to_gradingRGB(lms);

return (float4)(Y, rgb.x, rgb.y, LMS.w);
}

/* Same as above but compute only Yrg */
inline float4 gradingRGB_to_Yrg(float4 RGB)

inline float4 Yrg_to_LMS(const float4 Yrg)
{
float4 Yrg;
Yrg.x = fmax(0.67282368f * RGB.x + 0.47812261f * RGB.y + 0.01044966f * RGB.z, 0.f);
const float a = RGB.x + RGB.y + RGB.z;
RGB = (a == 0.f) ? (float4)0.f : RGB / a;

Yrg.y = RGB.x;
Yrg.z = RGB.y;
Yrg.w = RGB.w;
return Yrg;
const float Y = Yrg.x;

// reform rgb (normalized) from chroma
const float r = Yrg.y;
const float g = Yrg.z;
const float b = 1.f - r - g;
const float4 rgb = { r, g, b, 0.f };

// convert to lms (normalized)
const float4 lms = gradingRGB_to_LMS(rgb);

// denormalize to LMS
const float denom = (0.68990272f * lms.x + 0.34832189f * lms.y);
const float a = (denom == 0.f) ? 0.f : Y / denom;
return lms * a;
}

inline float4 Yrg_to_gradingRGB(const float4 Yrg)

/*
* Re-express Filmlight Yrg in polar coordinates Ych
*/

inline float4 Yrg_to_Ych(const float4 Yrg)
{
float4 RGB;
RGB.x = Yrg.y;
RGB.y = Yrg.z;
RGB.z = 1.f - Yrg.y - Yrg.z;
const float D65[4] = { 0.21962576f, 0.54487092f, 0.23550333f, 0.f };
const float Y = Yrg.x;
const float r = Yrg.y - D65[0];
const float g = Yrg.z - D65[1];
const float c = hypot(g, r);
const float h = atan2(g, r);
return (float4)(Y, c, h, Yrg.w);
}

const float a = (0.67282368f * RGB.x + 0.47812261f * RGB.y + 0.01044966f * RGB.z);
RGB = (a == 0.f) ? (float4)0.f : RGB * Yrg.x / a;
RGB.w = Yrg.w;
return RGB;

inline float4 Ych_to_Yrg(const float4 Ych)
{
const float D65[4] = { 0.21962576f, 0.54487092f, 0.23550333f, 0.f };
const float Y = Ych.x;
const float c = Ych.y;
const float h = Ych.z;
const float r = c * native_cos(h) + D65[0];
const float g = c * native_sin(h) + D65[1];
return (float4)(Y, r, g, Ych.w);
}
100 changes: 47 additions & 53 deletions data/kernels/extended.cl
View file
Open in desktop
Expand Up @@ -766,15 +766,6 @@ inline float soft_clip(const float x, const float soft_threshold, const float ha
}


inline float4 matrix_dot(const float4 vector, const float4 matrix[3])
{
float4 output;
output.x = dot(vector, matrix[0]);
output.y = dot(vector, matrix[1]);
output.z = dot(vector, matrix[2]);
return output;
}

kernel void
colorbalancergb (read_only image2d_t in, write_only image2d_t out,
const int width, const int height,
Expand All @@ -796,15 +787,27 @@ colorbalancergb (read_only image2d_t in, write_only image2d_t out,
if(x >= width || y >= height) return;
const float4 pix_in = read_imagef(in, sampleri, (int2)(x, y));

float4 Ych, RGB;
float4 XYZ_D65 = 0.f;
float4 LMS = 0.f;
float4 RGB = 0.f;
float4 Yrg = 0.f;
float4 Ych = 0.f;

Ych = fmax(pix_in, 0.f);
RGB = matrix_product_float4(Ych, matrix_in);
Ych = gradingRGB_to_Ych(RGB);
// clip pipeline RGB
RGB = fmax(pix_in, 0.f);

// go to CIE 2006 LMS D65
LMS = matrix_product_float4(RGB, matrix_in);

// go to Filmlight Yrg
Yrg = LMS_to_Yrg(LMS);

// go to Ych
Ych = Yrg_to_Ych(Yrg);

// Sanitize input : no negative luminance
float Y = fmax(Ych.x, 0.f);
const float4 opacities = opacity_masks(native_powr(Y, 0.4101205819200422f), // center middle grey in 50 %
Ych.x = fmax(Ych.x, 0.f);
const float4 opacities = opacity_masks(native_powr(Ych.x, 0.4101205819200422f), // center middle grey in 50 %
shadows_weight, highlights_weight, midtones_weight, mask_grey_fulcrum);
const float4 opacities_comp = (float4)1.f - opacities;

Expand All @@ -815,15 +818,20 @@ colorbalancergb (read_only image2d_t in, write_only image2d_t out,
if(Ych.z > M_PI_F) Ych.z -= 2.f * M_PI_F;
else if(Ych.z < -M_PI_F) Ych.z += 2.f * M_PI_F;

// Get max allowed chroma in working RGB gamut at current output hue
const float max_chroma_h = Y * lookup_gamut(gamut_lut, Ych.z);

// Linear chroma : distance to achromatic at constant luminance in scene-referred
const float chroma_boost = chroma_global + dot(opacities, chroma);
const float vib = vibrance * (1.0f - native_powr(Ych.y, fabs(vibrance)));
const float chroma_factor = fmax(1.f + chroma_boost + vib, 0.f);
Ych.y = soft_clip(Ych.y * chroma_factor, max_chroma_h, max_chroma_h * 4.f);
RGB = Ych_to_gradingRGB(Ych);
Ych.y *= chroma_factor;

// Go to Yrg
Yrg = Ych_to_Yrg(Ych);

// Go to LMS
LMS = Yrg_to_LMS(Yrg);

// Go to Filmlight RGB
RGB = LMS_to_gradingRGB(LMS);

// Color balance

Expand All @@ -837,36 +845,27 @@ colorbalancergb (read_only image2d_t in, write_only image2d_t out,
// midtones : power with sign preservation
RGB = sign(RGB) * native_powr(fabs(RGB) / white_fulcrum, midtones) * white_fulcrum;

// for the Y midtones power (gamma), we need to go in Ych again because RGB doesn't preserve color
Ych = gradingRGB_to_Yrg(RGB);
Y = Ych.x = native_powr(fmax(Ych.x / white_fulcrum, 0.f), midtones_Y) * white_fulcrum;
// for the non-linear ops we need to go in Yrg again because RGB doesn't preserve color
LMS = gradingRGB_to_LMS(RGB);
Yrg = LMS_to_Yrg(LMS);

// then the contrast
Y = Ych.x = grey_fulcrum * native_powr(Ych.x / grey_fulcrum, contrast);
RGB = Yrg_to_gradingRGB(Ych);
// Y midtones power (gamma)
Yrg.x = native_powr(fmax(Yrg.x / white_fulcrum, 0.f), midtones_Y) * white_fulcrum;

// Perceptual color adjustments
// Y fulcrumed contrast
Yrg.x = grey_fulcrum * native_powr(Yrg.x / grey_fulcrum, contrast);

// grading RGB to CIE 1931 XYZ 2° D65
const float4 RGB_to_XYZ_D65[3] = { { 1.64004888f, -0.10969806f, 0.49329934f, 0.f },
{ 0.61055787f, 0.47749658f, -0.08730269f, 0.f },
{ -0.10698534f, 0.07785058f, 1.66590006f, 0.f } };
LMS = Yrg_to_LMS(Yrg);
XYZ_D65 = LMS_to_XYZ(LMS);

const float4 XYZ_to_RGB_D65[3] = { { 0.54392489f, 0.14993776f, -0.15320716f, 0.f },
{ -0.68327274f, 1.88816348f, 0.30127843f, 0.f },
{ 0.06686186f, -0.07860825f, 0.57635773f, 0.f } };
// Perceptual color adjustments

// Go to JzAzBz for perceptual saturation
// We can't use gradingRGB_to_XYZ() since it also does chromatic adaptation to D50
// and JzAzBz uses D65, same as grading RGB. So we use the matrices above instead
float4 Jab;
RGB.w = Ych.w = 0.f; // init the 4th element for the dot product
Ych.xyz = matrix_dot(RGB, RGB_to_XYZ_D65).xyz;
Jab = XYZ_to_JzAzBz(Ych);
float4 Jab = XYZ_to_JzAzBz(XYZ_D65);

// Convert to JCh
float JC[2] = { Jab.x, hypot(Jab.y, Jab.z) }; // brightness/chroma vector
const float h = (JC[1] == 0.f) ? 0.f : atan2(Jab.z, Jab.y); // hue : (a, b) angle
const float h = atan2(Jab.z, Jab.y); // hue : (a, b) angle

// Project JC onto S, the saturation eigenvector, with orthogonal vector O.
// Note : O should be = (C * cosf(T) - J * sinf(T)) = 0 since S is the eigenvector,
Expand All @@ -892,24 +891,19 @@ colorbalancergb (read_only image2d_t in, write_only image2d_t out,
JC[0] = fmax(SO[0] * M_rot_inv[0][0] + SO[1] * M_rot_inv[0][1], 0.f);
JC[1] = fmax(SO[0] * M_rot_inv[1][0] + SO[1] * M_rot_inv[1][1], 0.f);

// Gamut mapping
const float out_max_chroma_h = JC[0] * lookup_gamut(gamut_lut, h);
JC[1] = soft_clip(JC[1], 0.8f * out_max_chroma_h, out_max_chroma_h);

// Project back to JzAzBz
Jab.x = JC[0];
Jab.y = JC[1] * native_cos(h);
Jab.z = JC[1] * native_sin(h);

Ych = JzAzBz_2_XYZ(Jab);
RGB.xyz = matrix_dot(Ych, XYZ_to_RGB_D65).xyz;
Ych = gradingRGB_to_Ych(RGB);
Y = fmax(Ych.x, 0.f);

// Gamut mapping
// Note : no need to check hue is in [-PI; PI], gradingRGB_to_Ych uses atan2f()
// which always returns angles in [-PI; PI]
const float out_max_chroma_h = Y * lookup_gamut(gamut_lut, Ych.z);
Ych.y = soft_clip(Ych.y, out_max_chroma_h, out_max_chroma_h * 4.f);
XYZ_D65 = JzAzBz_2_XYZ(Jab);

RGB = Ych_to_gradingRGB(Ych);
RGB = matrix_product_float4(RGB, matrix_out);
// Project back to D50 pipeline RGB
RGB = matrix_product_float4(XYZ_D65, matrix_out);

if(mask_display)
{
Expand Down