Color management of sources of illumination (stars) in an image

Question

I am working with images of star fields. For a given star I want to obtain a representative RGB color in a given image working colorspace, which may not be sRGB. (If I only had to worry about sRGB the problem is already solved.) I can look up the details of a star in stellar catalogues and obtain the star's B-V (this is a photometric measure, the important thing about it is that there is a formula for calculating the effective temperature of a star from its B-V). I can convert the effective temperature into a set of CIE1931 xy coordinates which then get converted into XYZ (again, how to do this is not in question). I'm using lcms2 to perform the color space transformations.

To provide a reasonably natural representation of star color the Sun is taken as white. This is effectively a white point of D58, whereas the working colorspace of the image may be sRGB or Rec2020 (white point D65) or Adobe RGB (white point D50) or potentially others.

I'm not entirely sure how to get the colors looking correct.

If I convert the XYZ to the working ICC profile using the Relative Colorimetric intent the white point looks wrong, I think because it's making the white point of the working profile act as the white point of the data whereas I want them to have their own white point, different to that of the image I'm placing them in.

Instead, should I first be converting XYZ to an intermediate RGB ICC profile with a D58 white point (effectively using that as a source profile) and then converting that RGB value to the image working profile using an Absolute Colorimetric intent?

Answer 1

OK, I'm going to answer this myself as I got it working.

Step 1: look up the B and V magnitudes of each star from a stellar catalogue (e.g. Simbad). I'll skip over this because it's not central to the question and the catalogues have reasonably documented query APIs.

Step 2: Calculate B-V for the star:

results[n].BV = star[i].bmag - star[i].vmag;

Step 3: Convert B-V to temperature. The simplest way is using the formula implemented by the following function (reference in the comment):

// Reference: https://en.wikipedia.org/wiki/Color_index and https://arxiv.org/abs/1201.1809 (Ballesteros, F. J., 2012)
// Uses Ballesteros' formula based on considering stars as black bodies
double bvToT(double bv) {
    double t = 4600. * ((1. / ((0.92 * bv) + 1.7)) + (1. / ((0.92 * bv) + 0.62)));
    return t;
}

You can try to look up the spectral type for each star and use separate bvToT calculations for each spectral type, but the above is generally good enough especially if you're dealing statistically with a large number of stars.

Step 4: Convert temperature to CIE xy: There are several approaches you can use here. The first is to use Mitchell Charity's Black body temperature to CIExy conversion tables. There are versions of these using the 1931 2 degree CMF with Judd-Vos corrections and the 1964 10 degree CMF and they are available for D65 and D58 white references. Implementation of a look-up table and interpolation is trivial. The tables claim validity from 1,000K to 40,000K. Secondly you can use a cubic spline method by Kim et al, as per the function below (reference in the comment):

// Returns a valid xyY for 1667K <= t <= 25000K, otherwise xyY = { 0.0 }
// Uses Kim et al's cubic spline Planckian locus (https://en.wikipedia.org/wiki/Planckian_locus)
static void temp_to_xyY(cmsCIExyY *xyY, cmsFloat64Number t) {
    // Calculate x
    if (t < 1667.0)
        xyY->x = 0.0;
    else if (t < 4000.0)
        xyY->x = (-0.2661239e9 / (t * t * t)) - (0.2343589e6 / (t * t)) +( 0.8776956e3 / t) + 0.179910;
    else if (t < 25000.0)
        xyY->x = (-3.0258469e9 / (t * t * t)) + (2.1070379e6 / (t * t)) + (0.2226347e3 / t) + 0.240390;
    else
        xyY->x = 0.0;

    cmsFloat64Number x = xyY->x;
    // Calculate y
    if (t < 1667)
        xyY->y = 0.0;
    else if (t < 2222.0)
        xyY->y = (-1.1063814 * x * x * x) - (1.34811020 * x * x) + (2.18555832 * x) - 0.20219683;
    else if (t < 4000.0)
        xyY->y = (-0.9549476 * x * x * x) - (1.37418593 * x * x) + (2.09137015 * x) - 0.16748867;
    else if (t < 25000.0)
        xyY->y = (3.0817580 * x * x * x) - (5.87338670 * x * x) + (3.75112997 * x) - 0.37001483;
    else
        xyY->y = 0.0;

    if (!(xyY->x == 0.0 && xyY->y == 0.0))
        xyY->Y = 1.0;
    else
        xyY->Y = 0.0;
}

The valid temperature range isn't as wide but it avoids storing a LUT. I'm starting to use the lcms2 API for datatypes now, because of what comes next. Finally, BQ Octantis proposed a formula on the Cloudy Nights astronomy forum. Reference provided in the comment:

// Returns a valid xyY for 2000K <= t, otherwise xyY = { 0.0 }
// Uses BQ Octantis's 6th order best fit for D65 values down to 2000K
// (https://www.cloudynights.com/topic/849382-generating-a-planckian-ccm/)
static void BQ_temp_to_xyY(cmsCIExyY *xyY, cmsFloat64Number t) {
    if (t < 2000.0) {
        xyY->x = xyY->y = xyY->Y = 0.0;
    } else {
        xyY->x = -1.3737e-25 * pow(t,6.0) + 3.0985e-22 * pow(t, 5.0) + 1.5842e-16 * pow(t, 4.0)
                - 4.0138e-12 * pow(t, 3.0) + 4.3777e-08 * pow(t, 2.0) - 2.4363e-04 * t + 8.7309e-01;
        cmsFloat64Number x = xyY->x;
        xyY->y =  2.5110e+01 * pow(x, 5.0) - 5.9883e+01 * pow(x, 4.0) + 5.5545e+01 * pow(x, 3.0)
                - 2.7667e+01 * pow(x, 2.0) + 8.1167e+00 * x - 7.0462e-01;
        xyY->Y = 1.0;
    }
}

This is a single 6th order spline, but requires expensive function calls to pow() to compute.

Step 5: Convert xyY to XYZ: This is trivial using lcms2:

cmsxyY2XYZ(&XYZ, &xyY);

where xyY is the black body color you calculated at step 4.

Step 6: Adapt the chromaticity: the Planckian locus formulas are based on D65 whereas the lcms2 built in XYZ profile is based on D50 so we do a chromatic adaptation:

const cmsCIEXYZ D65 = {0.95045471, 1.0, 1.08905029};
const cmsCIEXYZ D50 = {0.964199999, 1.000000000, 0.824899998};
cmsCIXYZ XYZ, XYZ_adapted;
cmsAdaptToIlluminant(&XYZ_adapted, &D65, &D50, &XYZ);

Step 7: convert to working color space RGB: Finally, create a cmsHTRANSFORM from XYZ to your chosen working color space. Assume your working profile is cmsHPROFILE profile...

float xyz[3], rgb[3] = { 0.f };
xyz[0] = (float) XYZ_adapted.X;
xyz[1] = (float) XYZ_adapted.Y;
xyz[2] = (float) XYZ_adapted.Z;

cmsHPROFILE xyzprofile = cmsCreateXYZProfile();
cmsHTRANSFORM transform = cmsCreateTransformTHR(com.icc.context_single,
            xyzprofile,
            TYPE_XYZ_FLT,
            profile,
            TYPE_RGB_FLT,
            INTENT_RELATIVE_COLORIMETRIC,
            cmsFLAGS_NONEGATIVES);
cmsCloseProfile(profile);
cmsCloseProfile(xyzprofile);
cmsDoTransform(transform, &xyz, &rgb, 1);
cmsFloat64Number maxval = max(max(rgb[0], rgb[1]), rgb[2]);
rgb[0] /= maxval;
rgb[1] /= maxval;
rgb[2] /= maxval;

Alternatives: There are also straightforward lookup tables that convert directly from B-V values to sRGB R, G and B values. If all you care about is sRGB then they're absolutely fine, however while the method sketched out above is a little more complex it exposes all the workings and also allows you to generate R, G and B values (or even CMYK values, I guess) that are correct for whatever color space you are working in.