How to get XY value from ct in Philips Hue?

Question

How to get XY value from ct.

Ex: ct = 217, I want to get x="0.3127569", y= "0.32908".

I'm able to convert XY value into ct value using this below code.

        float R1 = [hue[0] floatValue];
        float S1 = [hue[1] floatValue];
        
        float result = ((R1-0.332)/(S1-0.1858));
        
        NSString *ctString = [NSString stringWithFormat:@"%f", ((-449*result*result*result)+(3525*result*result)-(6823.3*result)+(5520.33))];
        
        float micro2 = (float) (1 / [ctString floatValue] * 1000000);
        
        NSString *ctValue = [NSString stringWithFormat:@"%f", micro2];
        ctValue = [NSString stringWithFormat:@"%d", [ctValue intValue]];
        if ([ctValue integerValue] < 153) {
            ctValue = [NSString stringWithFormat:@"%d", 153];
        }

Now I want reverse value, which is from ct to XY.

Answer 1

On Phillips HUE
2000K maps to 500 and 6500K maps to 153
given in ct as color temperature but can be thought as actually being Mired.
Mired means micro reciprocal degree wikipedia.

ct is possibly used because it is not 100% Mired. Quite sure Phillips uses a lookup table as a lot CIE algorithms do because there are just 347 indexes in this range from 153 to 500.

The following is not a solution, it's just simple concept of a lookup table. And as the CIE 1931 xy to CCT Formula by McCamy suggests found here it is possible to use a lookup table to find x and y as well. A table can be found here but i am not sure if that is the right lookup table.

reminder so the following is not a solution, but to find an reverse algo the code may help.

typedef int Kelvin;
typedef float Mired;

Mired linearMiredByKelvin(Kelvin k) {
    if (k==0) return 0;
    return 1000000.0/k;
}

-(void)mired {
    Mired miredMin = 2000.0/13.0; // 153,84 = reciprocal 6500K
    Mired miredMax = 500.0;       // 500,00 = reciprocal 2000K
    
    Mired lookupMiredByKelvin[6501]; //max 6500 Kelvin + 1 safe index
    //Kelvin lookupKelvinByMired[501]; //max 500 Mired + 1 safe index
    
    // dummy stuff, empty unused table space
    for (Kelvin k = 0; k < 2000; k++) {
        lookupMiredByKelvin[k] = 0;
    }
    //for (Mired m = 0.0; m < 154.0; m++) {
    //    lookupKelvinByMired[(int)m] = 0;
    //}
    
    for (Kelvin k=2000; k<6501; k++) {
        Mired linearMired = linearMiredByKelvin(k);
        float dimm = (linearMired - miredMin) / ( miredMax - miredMin);
        Kelvin ct = (Kelvin)(1000000.0/(dimm*miredMax - dimm*miredMin + miredMin));
        lookupMiredByKelvin[k] = linearMiredByKelvin(ct);
        if (k==2000 || k==2250 || k==2500 || k==2750 ||
            k==3000 || k==3250 || k==3500 || k==3750 ||
            k==4000 || k==4250 || k==4500 || k==4750 ||
            k==5000 || k==5250 || k==5500 || k==5750 ||
            k==6000 || k==6250 || k==6500 || k==6501 )
            fprintf(stderr,"%d %f %f\n",ct, dimm, lookupMiredByKelvin[k]);
    }
}

at least this is proof that x and y will not sit on a simple vector.

CCT means correlated colour temperature and like the implementation in the question shows can be calculated via n= (x-0.3320)/(0.1858-y); CCT = 437*n^3 + 3601*n^2 + 6861*n + 5517. (after McCamy)

but a cct=217 is out of range of above link'ed lookup table.

following the idea in this git-repo from colour-science and ported to C it could look like..

void CCT_to_xy_CIE_D(float cct) {
    //if (CCT < 4000 || CCT > 25000) fprintf(stderr, "Correlated colour temperature must be in domain, unpredictable results may occur! \n");
    float x =  calculateXviaCCT(cct);
    float y = calculateYviaX(x);
    NSLog(@"cct=%f x%f y%f",cct,x,y);
}

float calculateXviaCCT(float cct) {
    float cct_3 = pow(cct, 3); //(cct*cct*cct);
    float cct_2 = pow(cct, 2); //(cct*cct);
    if (cct<=7000)
    return -4.607 * pow(10, 9) / cct_3 + 2.9678 * pow(10, 6) / cct_2 + 0.09911 * pow(10, 3) / cct + 0.244063;
    return -2.0064 * pow(10, 9) / cct_3 + 1.9018 * pow(10, 6) / cct_2 + 0.24748 * pow(10, 3) / cct + 0.23704;
}
float calculateYviaX(float x) {
    return -3.000 * pow(x, 2) + 2.870 * x - 0.275;
}

CCT_to_xy_CIE_D(6504.38938305); //proof of concept
//cct=6504.389160 x0.312708 y0.329113

CCT_to_xy_CIE_D(217.0); 
//cct=217.000000 x-387.131073 y-450722.750000
// so for sure Phillips hue temperature given in ct between 153-500 is not a good starting point

//but
CCT_to_xy_CIE_D(2000.0);
//cct=2000.000000 x0.459693 y0.410366

this seems to work fine with CCT between 2000 and 25000, but maybe confusing is CCT is given in Kelvin here.

Answer 2

EDIT

This has been through so many revisions and ideas. To keep it simple I edited most of that out and just give you the final result.

This fits your function perfectly except for a region in the middle (temp from 256 to 316) where it deviates a bit.

The problem with your function is that it has approximately infinite solutions, so to solve it nicely you need more constraints, but what? Ol Sen's reference https://www.waveformlighting.com/tech/calculate-color-temperature-cct-from-cie-1931-xy-coordinates discusses it in some detail and then mentions that you want a Duv to be zero. It also gives a way to calculate Duv and so I added that to my optimiser and voila!

Nice and smooth. The optimiser now solves for x and y that both satisfies your function and also minimises Duv.

To get it to work nicely I had to scale Duv quite a bit. That page mentions that Duv should be very small so I think this is a good thing. Also, as the temp increases the scaling should to help the optimiser.

Below prints from 153 to 500.

#import <Foundation/Foundation.h>

// Function taken from your code
// Simplified a bit
int ctFuncI ( float x, float y )
{
    // float R1 = [hue[0] floatValue];
    // float S1 = [hue[1] floatValue];
    float result = (x-0.332)/(y-0.1858);
    float cubic  = - 449 * result * result * result + 3525 * result * result - 6823.3 * result + 5520.33;
    float micro2 = 1 / cubic * 1000000;
    int ct = ( int )( micro2 + 0.5 );
    
    if ( ct < 153 )
    {
        ct = 153;
    }
    
    return ct;
}

// Need this
// Float version of your code
float ctFuncF ( float x, float y )
{
    // float R1 = [hue[0] floatValue];
    // float S1 = [hue[1] floatValue];
    float result = (x-0.332)/(y-0.1858);
    float cubic  = - 449 * result * result * result + 3525 * result * result - 6823.3 * result + 5520.33;

    return 1000000 / cubic;
}

// We need an additional constraint
// https://www.waveformlighting.com/tech/calculate-duv-from-cie-1931-xy-coordinates
// Given x, y calculate Duv
// We want this to be 0
float duv ( float x, float y )
{
    float f = 1 / ( - 2 * x + 12 * y + 3 );
    float u = 4 * x * f;
    float v = 6 * y * f;

    // I'm typing float but my heart yells double
    float k6 = -0.00616793;
    float k5 =  0.0893944;
    float k4 = -0.5179722;
    float k3 =  1.5317403;
    float k2 = -2.4243787;
    float k1 =  1.925865;
    float k0 = -0.471106;
    
    float du = u - 0.292;
    float dv = v - 0.24;

    float Lfp = sqrt ( du * du + dv * dv );
    float a = acos( du / Lfp );
    float Lbb = k6 * pow ( a, 6 ) + k5 * pow( a, 5 ) + k4 * pow( a, 4 ) + k3 * pow( a, 3 ) + k2 * pow(a,2) + k1 * a + k0;

    return Lfp - Lbb;
}

// Solver!
// Returns iterations
int ctSolve ( int ct, float * x, float * y )
{
    int iter = 0;
    float dx = 0.001;
    float dy = 0.001;

    // Error
    // Note we scale duv a bit
    // Seems the higher the temp, the higher scale we require
    // Also note the jump at 255 ...
    float s = 1000 * ( ct > 255 ? 10 : 1 );
    float d = fabs( ctFuncF ( * x, * y ) - ct ) + s * fabs( duv ( * x, * y ) );
    
    // Approx
    while ( d > 0.5 && iter < 250 )
    {
        iter ++;

        dx *= fabs( ctFuncF ( * x + dx, * y ) - ct ) + s * fabs( duv ( * x + dx, * y ) ) < d ? 1.2 : - 0.5;
        dy *= fabs( ctFuncF ( * x, * y + dy ) - ct ) + s * fabs( duv ( * x, * y + dy ) ) < d ? 1.2 : - 0.5;

        * x += dx;
        * y += dy;

        d = fabs( ctFuncF ( * x, * y ) - ct ) + s * fabs( duv ( * x, * y ) );
    }

    return iter;
}

// Tester
int main(int argc, const char * argv[]) {
    
    @autoreleasepool
    {
        // insert code here...
        NSLog(@"Hello, World!");
        float x, y;
        int sume = 0;
        int sumi = 0;
                
        for ( int ct = 153; ct <= 500; ct ++ )
        {
            // Initial guess
            x = 0.4;
            y = 0.4;

            // Approx
            int iter = ctSolve ( ct, & x, & y );
            
            // CT and error
            int ctEst = ctFuncI ( x, y );
            int e = ct - ctEst;
            
            // Diagnostics
            sume += abs ( e );
            sumi += iter;
            
            // Print out results
            NSLog ( @"want ct = %d x = %f y = %f got ct %d in %d iter error %d", ct, x, y, ctEst, iter, e );
        }
        
        NSLog ( @"Sum of abs errors %d iterations %d", sume, sumi );

    }
    
    return 0;
}

To use it, do as below.

// To call it, init x and y to some guess
float x = 0.4;
float y = 0.4;

// Then call solver with your temp
int ct = 217;

ctSolve( ct, & x, & y ); // Note you pass references to x and y

// Done, answer now in x and y

Answer 3

a bit more compact answer and functions to convert back and forth..
beware there are rounding issues because McCamy's formula relies and mathematical assumptions. And so the backward calculation does also.

if you want to find more results search directly for "n= (x-0.3320)/(0.1858-y); CCT = 437*n^3 + 3601*n^2 + 6861*n + 5517" there are plenty of different methods to convert back and forth.

so here Phillips-Hue @[@x,@y] to Phillips-ct,
Phillips-ct to CCT,
CCT to x,y

void CCT_to_xy_CIE_D(float cct) {
    //if (CCT < 4000 || CCT > 25000) fprintf(stderr, "Correlated colour temperature must be in domain, unpredictable results may occur! \n");
    float x = calculateXviaCCT(cct);
    float y = calculateYviaX(x);
    fprintf(stderr,"cct=%f x%f y%f",cct,x,y);
}
float calculateXviaCCT(float cct) {
    float cct_3 = pow(cct, 3); //(cct*cct*cct);
    float cct_2 = pow(cct, 2); //(cct*cct);
    if (cct<=7000.0)
    return -4.607 * pow(10, 9) / cct_3 + 2.9678 * pow(10, 6) / cct_2 + 0.09911 * pow(10, 3) / cct + 0.244063;
    return -2.0064 * pow(10, 9) / cct_3 + 1.9018 * pow(10, 6) / cct_2 + 0.24748 * pow(10, 3) / cct + 0.23704;
}
float calculateYviaX(float x) {
    return -3.000 * x*x + 2.870 * x - 0.275;
}
int calculate_PhillipsHueCT_withCCT(float cct) {
    if (cct>6500.0) return 2000.0/13.0;
    if (cct<2000.0) return 500.0;
    //return (float) (1 / cct * 1000000); // same as..
    return 1000000 / cct;
}
float calculate_CCT_withPhillipsHueCT(float ct) {
    if (ct == 0.0) return 0.0;
    return 1000000 / ct;
}
float calculate_CCT_withHueXY(NSArray *hue) {
    float x = [hue[0] floatValue]; //R1
    float y = [hue[1] floatValue]; //S1
    //x = 0.312708; y = 0.329113;
    float n = (x-0.3320)/(0.1858-y);
    float cct = 437.0*n*n*n + 3601.0*n*n + 6861.0*n + 5517.0;
    return cct;
}

// MC Camy formula n=(x-0.3320)/(0.1858-y); cct = 437*n^3 + 3601*n^2 + 6861*n + 5517;
-(void)testPhillipsHueCt_backAndForth {
    NSArray *hue = @[@(0.312708),@(0.329113)];
    
    float cct = calculate_CCT_withHueXY(hue);

    float ct = calculate_PhillipsHueCT_withCCT(cct);
    NSLog(@"ct %f",ct);
    
    CCT_to_xy_CIE_D(cct); // check
    
    CCT_to_xy_CIE_D(6504.38938305); //proof of concept

    CCT_to_xy_CIE_D(2000.0);
    
    CCT_to_xy_CIE_D(calculate_CCT_withPhillipsHueCT(217.0));
}