darktable-org / darktable

darktable is an open source photography workflow application and raw developer
https://www.darktable.org
GNU General Public License v3.0
9.76k stars 1.14k forks source link

color calibration to report original illuminant CCT #14630

Open kofa73 opened 1 year ago

kofa73 commented 1 year ago

Is your feature request related to a problem? Please describe. Users are frequently confused and ask about the CCT dispalyed by color calibration. For example, for bright daylight shots, the displayed CCT is way below 5000 K, while white balance would give a more realistic reading of over 5000 K. The difference is about 1000 -1500 K.

If I understand correctly, this happens because the CCT displayed by color calibration is not the original, but a mapped value. If my understanding is right, the 'modern' chromatic adaptation process first maps the D65 white point to D50; then color calibration is used to map the mapped scene white point to D50; the CCT displayed is not the CCT of the original illuminant, but rather where it ended up after the D65 -> D50 conversion. That value, I think, may be useful for developers and for the purposes of the algorithm, but is not useful, is even confusing, for photographers who are used to thinking in terms of the scene illuminant.

Quoting colour picker values from https://discuss.pixls.us/t/confused-about-d50-d65-and-cct-in-white-balance-and-color-calibration-modules/37293/10, which uses a photo taken under controlled studio conditions, with calibrated (5000 - 5100 K) lighting, it seems the following mapping is carried out: 6500 K → 5000 K 5341 K → 4321 K 5216 K → 4216 K

With the 'modern' method: image With the 'legacy' method (white balance only): image

The documentation even says: To the left of the color patch is the CCT (correlated color temperature) approximation. This is the closest temperature, in kelvin, to the illuminant currently in use. That is something most would interpret as the 'unadapted', original scene illuminant.

I do not know if the later evaluations done by the module, to determine if the illuminant is close to the Planckian or Daylight locus, take this into account or not. If not, then they are seriously off.

Describe the solution you'd like Have the module display the 'unadapted'/'unmapped' CCT.

Alternatives Keep as-is, document that the CCT has been mapped and is not what a colorimeter would show.

fxthomas commented 1 year ago

This got me looking into it, because that actually explains some unintuitive results I had in color calibration a while back for scenes where I mostly knew the illuminant.

The temperature for the "color calibration" module is computed in src/common/illuminants.h:

static inline void WB_coeffs_to_illuminant_xy(const float CAM_to_XYZ[4][3], const dt_aligned_pixel_t WB,
                                              float *x, float *y)
{
  // Find the illuminant chromaticity x y from RAW WB coeffs and camera input matrice
  dt_aligned_pixel_t XYZ, LMS;
  // Simulate white point, aka convert (1, 1, 1) in camera space to XYZ
  // warning : we multiply the transpose of CAM_to_XYZ  since the pseudoinverse transposes it
  XYZ[0] = CAM_to_XYZ[0][0] / WB[0] + CAM_to_XYZ[1][0] / WB[1] + CAM_to_XYZ[2][0] / WB[2];
  XYZ[1] = CAM_to_XYZ[0][1] / WB[0] + CAM_to_XYZ[1][1] / WB[1] + CAM_to_XYZ[2][1] / WB[2];
  XYZ[2] = CAM_to_XYZ[0][2] / WB[0] + CAM_to_XYZ[1][2] / WB[1] + CAM_to_XYZ[2][2] / WB[2];

  // Matrices white point is D65. We need to convert it for darktable's pipe (D50)
  static const dt_aligned_pixel_t D65 = { 0.941238f, 1.040633f, 1.088932f, 0.f };
  const float p = powf(1.088932f / 0.818155f, 0.0834f);

  convert_XYZ_to_bradford_LMS(XYZ, LMS);
  bradford_adapt_D50(LMS, D65, p, FALSE, LMS);
  convert_bradford_LMS_to_XYZ(LMS, XYZ);

  // Get white point chromaticity
  XYZ[0] /= XYZ[1];
  XYZ[2] /= XYZ[1];
  XYZ[1] /= XYZ[1];

  *x = XYZ[0] / (XYZ[0] + XYZ[1] + XYZ[2]);
  *y = XYZ[1] / (XYZ[0] + XYZ[1] + XYZ[2]);
}

For the "white balance" module it's located in src/iop/temperature.c:

static cmsCIEXYZ mul2xyz(dt_iop_module_t *self,
                         const dt_iop_temperature_params_t *p)
{
  dt_iop_temperature_gui_data_t *g = (dt_iop_temperature_gui_data_t *)self->gui_data;

  double CAM[4];
  _temp_array_from_params(CAM, p);

  for(int k = 0; k < 4; k++)
    CAM[k] = CAM[k] > 0.0f ? 1.0 / CAM[k] : 0.0f;

  double XYZ[3];
  for(int k = 0; k < 3; k++)
  {
    XYZ[k] = 0.0;
    for(int i = 0; i < 4; i++)
    {
      XYZ[k] += g->CAM_to_XYZ[k][i] * CAM[i];
    }
  }

  return (cmsCIEXYZ){ XYZ[0], XYZ[1], XYZ[2] };
}

Both follow a fairly simple algorithm : the wb coefficients stored by the camera convert the raw sensor RGB values into equal RGB values for the illuminant ; the wb matrix is supposed to handle XYZ conversions for the entire colorspace.

For a pixel corresponding to the scene illuminant we have both $M \cdot (R G B) = (X Y Z)$ but also $k \cdot (R G B) = 1$. Therefore, $(X Y Z) = M \cdot (1 / k)$, which can then be converted into a CCT.

That's exactly what the "white balance module" is doing, because it's converting from RAW domain to D65 XYZ.

The "color calibration" module however, is after the "input profile" module, which puts the pipeline into the D50 illuminant by default (if I'm understanding this correctly). The xy coordinates are therefore adapted using a CAT (Bradford) which is directly visible in the code... and you're right that from a developer PoV it makes sense, because that's how the pipeline works.

fxthomas commented 1 year ago

Now on to the CCT conversion there's also this (used in "white balance") which uses an ad-hoc binary search:

// binary search inversion
static void XYZ_to_temperature(cmsCIEXYZ XYZ, float *TempK, float *tint)
{
  double maxtemp = DT_IOP_HIGHEST_TEMPERATURE, mintemp = DT_IOP_LOWEST_TEMPERATURE;
  cmsCIEXYZ _xyz;

  for(*TempK = (maxtemp + mintemp) / 2.0;
      (maxtemp - mintemp) > 1.0;
      *TempK = (maxtemp + mintemp) / 2.0)
  {
    _xyz = temperature_to_XYZ(*TempK);
    if(_xyz.Z / _xyz.X > XYZ.Z / XYZ.X)
      maxtemp = *TempK;
    else
      mintemp = *TempK;
  }

  // TODO: Fix this to move orthogonally to planckian locus
  *tint = (_xyz.Y / _xyz.X) / (XYZ.Y / XYZ.X);

  if(*TempK < DT_IOP_LOWEST_TEMPERATURE) *TempK = DT_IOP_LOWEST_TEMPERATURE;
  if(*TempK > DT_IOP_HIGHEST_TEMPERATURE) *TempK = DT_IOP_HIGHEST_TEMPERATURE;
  if(*tint < DT_IOP_LOWEST_TINT) *tint = DT_IOP_LOWEST_TINT;
  if(*tint > DT_IOP_HIGHEST_TINT) *tint = DT_IOP_HIGHEST_TINT;
}

...and this, used in "color calibration", which uses method from Hernandez-Andres, Lee & Romero (1999) (the link is dead, but I could hunt for the same model in the colour-science library):

static inline float xy_to_CCT(const float x, const float y)
{
  // Try to find correlated color temperature from chromaticity
  // Valid for 3000 K to 50000 K
  // Reference : https://www.usna.edu/Users/oceano/raylee/papers/RLee_AO_CCTpaper.pdf
  // Warning : we throw a number ever if it's grossly off. You need to check the error later.
  if(x < FLT_MAX)
  {
    const float n = (x - 0.3366f)/(y - 0.1735f);
    return (-949.86315f + 6253.80338f * expf(-n / 0.92159f)
            + 28.70599f * expf(-n / 0.20039f)
            + 0.00004f * expf(-n / 0.07125f));
  }
  else // we were called with coordinates flagged as invalid
    return 0.0f; // invalid chromaticity
}

They "should" be similar, but I wouldn't be surprised to also find differences between them given that they use very different methods with some remaining "fix me" comments ;)

github-actions[bot] commented 1 year ago

This issue has been marked as stale due to inactivity for the last 60 days. It will be automatically closed in 300 days if no update occurs. Please check if the master branch has fixed it and report again or close the issue.

kanyck commented 10 months ago

Is there any progress on this? It's really confusing. I don't look at this CCT numbers anymore, because they show nothing practically usable. I expect the numbers to reflect the scene illuminant with reasonable accuracy, or they shouldn't appeared in the UI at all. Wrong info is worse then no info.

github-actions[bot] commented 8 months ago

This issue has been marked as stale due to inactivity for the last 60 days. It will be automatically closed in 300 days if no update occurs. Please check if the master branch has fixed it and report again or close the issue.