How to find cct of LEDs in cie 1960 graph?

Hello,

I found numerous similar questions to this one but none gave the answer I need.

As the title says, I am trying to find cct of LEDs using the algorithm presented here (slides 11-14). The problem I am having is finding the orthogonal projection onto the black-body line. I do not know how to program it. The graph with the LEDs(red) is shown below.

Is it even possible to find cct of every LED in these positions?

Thank you

Answer 1

[John D'Errico]

I don't have the data you show for those curves. However, an orthogonal projection of a point onto a curve in space is not that difficult. Well, it is possible to do, at least. ;-)

The point is that an orthogonal projection for a point onto a curve merely requires finding the point on said curve that lies at a minimal distance. To project one curve onto another is merely the act of doing the above many times. In fact, I've posted a simple utility on the file exchange which does exactly that for you. It will take me more time to create some data for this, than it will to show how to use the utility though. Oh well.

So, first a couple of curves. Lets see, I'll just create two elliptical arcs in the plane.

t1 = linspace(-.5,2,30)';
xy1 = [cos(t1),sin(t1)]*[-1 2;-2 1] + [.3 .3];
plot(xy1(:,1),xy1(:,2),'ko-')
axis equal

hold on
t2 = linspace(0.75,2.75,100)';
xy2 = [cos(t2),sin(t2)]*.6 + [-1 1.6];
plot(xy2(:,1),xy2(:,2),'r.')
grid on

Yes, easier to hack up some data than it is to actually find that black body curve online, and I'm just too lazy to think today.

Now, for any point on the red line, I'll project it onto the curve in black. The orthogonal line of projection will be in blue here.

[xyproj,distance] = distance2curve(xy1,xy2,'spline');

plot([xy2(:,1),xyproj(:,1)]',[xy2(:,2),xyproj(:,2)]','b-')

As you can see, there can easily be a discontinuous jump for that projection, but what can you expect given two nonlinear arcs? I doubt you care to know what the distance for that projection was, but the second output argument tells you that information.

A spline interpolant was used for the black curve, so I solve the equations to find the actual minimal distance to a curved arc of a spline in the code. Not that difficult, but you need to understand splines.

Find distance2curve on the file exchange for free download, here:

https://www.mathworks.com/matlabcentral/fileexchange/34869-distance2curve

[Leonardo Max Golusin]

The reason why I need orthogonal projection is because the algorithm for finding cct requires the line to the desired point to be orthogonal to the planckian locus. Otherwise, your function would work perfectly.

[John D'Errico]

Then you needed to be far more clear about what you wanted. In fact, your question is still not remotely clear.

The closest point to a curved line will have a projection that is ALWAYS orthognal to the curved line. If it was not, then there would be another point on the curved line that was closer. Do I need to prove this statement? As such, if that is your only requirement, then distance2curve does exactly what you want. The only time there is an issue is when the curve to be projected onto is limited in extent, so it ends at some point. Then the closest point on the curve may merely be an end point. And of course, if the curve itself is not differentiable, then if it has a cusp on it somewhere that will again be a problem. However, I explicitly told you to use the spline option for distance2curve.

Be VERY careful, because if you do not scale the axes properly to have equal units, then even two lines that intersect at right angles will not appear to do so. Do I need to prove that statement too? I hope not.

So what is your question? Be clear.

plot([0 0;1 1],[1 0;0 1],'-')
axis([0 1 -1 2])

Do those lines look orthogonal to you? I hope not, yet they are provably so. Differential axis scaling will cause this.

Answer 2

[Cloud Yan]

Please ref a python project colour-science

The ohno's method is shown as:

https://colour.readthedocs.io/en/develop/_modules/colour/temperature/ohno2013.html#uv_to_CCT_Ohno2013

An typical method to get CCT is Roberson's method.

It is a pretty clear code, and you can easily modify to matlab codes.

Here is an example https://frudawski.de/ciexy2cct/

The paper is downloaded here https://opg.optica.org/josa/abstract.cfm?uri=josa-58-11-1528

Here is my matlab function:

function CCT = uv_to_CCT_Robertson1968v(ut,vt)
% input：CIE 1960 chromaticity coordinates ut and vt,
%        ut is a column vector, so as vt.
% output：Correlated Colour Temperature (CCT), a column vector.
%     Define the data for a *Roberston (1968)* iso-temperature line.
%     r Temperature T=1/r*MK^(-1) in reciprocal mega-kelvin degrees.
%     u chromaticity coordinate of the temperature T.
%     v	chromaticity coordinate of the temperature T.
%     t Slope of the u-v chromaticity coordinate.
ruvt =[
0   0.18006 0.26352 -0.24341 % ∞ K
1	0.18012	0.26375	-2.4443E-01
5	0.18035	0.26468	-2.4878E-01
10	0.18066	0.26589	-2.5479E-01
15	0.18098	0.26715	-2.6145E-01
20	0.18133	0.26846	-2.6876E-01
25	0.18169	0.26980	-2.7674E-01
30	0.18208	0.27119	-2.8539E-01
35	0.18250	0.27261	-2.9470E-01
40	0.18293	0.27407	-3.0470E-01
45	0.18340	0.27557	-3.1538E-01
50	0.18388	0.27709	-3.2675E-01
55	0.18440	0.27864	-3.3881E-01
60	0.18494	0.28021	-3.5156E-01
65	0.18551	0.28181	-3.6500E-01
70	0.18611	0.28342	-3.7915E-01
75	0.18674	0.28504	-3.9400E-01
80	0.18740	0.28668	-4.0955E-01
85	0.18808	0.28832	-4.2580E-01
90	0.18880	0.28997	-4.4278E-01
95	0.18954	0.29162	-4.6047E-01
100	0.19032	0.29326	-4.7888E-01
112.5	0.19238	0.29737	-5.2811E-01
125	0.19462	0.30141	-5.8204E-01
137.5	0.19704	0.30536	-6.4082E-01
150	0.19962	0.30921	-7.0471E-01
162.5	0.20236	0.31291	-7.7398E-01
175	0.20525	0.31647	-8.4901E-01
187.5	0.20827	0.31988	-9.3024E-01
200	0.21142	0.32312	-1.0182E+00
212.5	0.21469	0.32619	-1.1135E+00
225	0.21807	0.32909	-1.2168E+00
237.5	0.22155	0.33182	-1.3291E+00
250	0.22511	0.33439	-1.4512E+00
262.5	0.22876	0.33679	-1.5844E+00
275	0.23247	0.33904	-1.7298E+00
287.5	0.23626	0.34113	-1.8890E+00
300	0.24010	0.34308	-2.0637E+00
312.5	0.24399	0.34488	-2.2559E+00
325	0.24792	0.34655	-2.4681E+00
337.5	0.25190	0.34809	-2.7031E+00
350	0.25591	0.34951	-2.9641E+00
362.5	0.25994	0.35081	-3.2553E+00
375	0.26400	0.35200	-3.5814E+00
387.5	0.26808	0.35308	-3.9483E+00
400	0.27218	0.35407	-4.3633E+00
412.5	0.27628	0.35496	-4.8355E+00
425	0.28039	0.35577	-5.3762E+00
437.5	0.28451	0.35649	-6.0001E+00
450	0.28863	0.35714	-6.7262E+00
462.5	0.29274	0.35772	-7.5799E+00
475	0.29685	0.35823	-8.5955E+00
487.5	0.30095	0.35868	-9.8205E+00
500	0.30505	0.35907	-1.1324E+01
512.5	0.30913	0.35940	-1.3207E+01
525	0.31320	0.35968	-1.5628E+01
537.5	0.31725	0.35992	-1.8847E+01
550	0.32129	0.36011	-2.3325E+01
562.5	0.32531	0.36027	-2.9960E+01
575	0.32931	0.36038	-4.0770E+01
587.5	0.33328	0.36046	-6.1433E+01
600	0.33724	0.36051	-1.1645E+02 % 62 line
625	0.34508	0.36053	1.8978E+02  % 63 line
650	0.35281	0.36044	5.7748E+01
675	0.36044	0.36026	3.5979E+01
700	0.36795	0.36002	2.7095E+01
725	0.37535	0.35972	2.2303E+01
750	0.38262	0.35937	1.9328E+01
775	0.38976	0.35897	1.7318E+01
800	0.39677	0.35855	1.5879E+01
825	0.40366	0.35810	1.4807E+01
850	0.41040	0.35763	1.3985E+01
875	0.41701	0.35715	1.3339E+01
900	0.42349	0.35665	1.2822E+01
925	0.42983	0.35615	1.2403E+01
950	0.43603	0.35564	1.2059E+01
975	0.44209	0.35513	1.1773E+01
1000	0.44801	0.35462	1.1535E+01
1050	0.45945	0.35362	1.1163E+01
1100	0.47034	0.35263	1.0894E+01
1150	0.48069	0.35167	1.0695E+01
1200	0.49052	0.35074	1.0546E+01
1250	0.49984	0.34985	1.0434E+01
1300	0.50865	0.34901	1.0348E+01
1350	0.51698	0.34820	1.0282E+01
1400	0.52485	0.34743	1.0231E+01
1450	0.53226	0.34671	1.0191E+01
1500	0.53924	0.34602	1.0160E+01
1550	0.54580	0.34537	1.0135E+01
1600	0.55196	0.34476	1.0115E+01
1650	0.55775	0.34419	1.0100E+01
1700	0.56316	0.34366	1.0087E+01
1750	0.56823	0.34315	1.0077E+01
1800	0.57296	0.34268	1.0069E+01
1850	0.57738	0.34224	1.0062E+01
1900	0.58150	0.34183	1.0057E+01
1950	0.58533	0.34145	1.0052E+01
2000	0.58888	0.34110	1.0048E+01 % 500 K
];
r = ruvt(:,1);
u = ruvt(:,2);
v = ruvt(:,3);
t = ruvt(:,4);

N = length(vt);% number of test points
% calculate the distance of (ut,vt) to the iso-temperature line: y-v=t(x-u)
% t is slope of the iso-temperature line
% (u,v) is the intersection point of Planckian locus and iso-temperature line
% distance = |vt-v -t(ut-u)|/(t^2+1)^0.5
d = zeros(length(r),N);% 98 x N
for  i = 1:length(vt)
    d(:,i) = ( (vt(i) - v) - t .* (ut(i) - u) )./sqrt(t.^2+1);
end

d(63:end,:) = -1*d(63:end,:); % force to be positive, as the slopes of the 63rd to 98th iso-temperature line is negtive.
[~,I] = min(abs(d));% 1 x N
Ip = 0*I;
d_pre = 0*I;
d_nex = 0*I;
for i = 1:N
    if d(I(i),i)<=0
        Ip(i) = I(i)-1;
        d_pre(i) = d(Ip(i)  ,i);
        d_nex(i) = d(Ip(i)+1,i);
    else
        Ip(i) = I(i);
        d_pre(i) = d(Ip(i)  ,i);
        d_nex(i) = d(Ip(i)+1,i);
    end
end

p = d_pre ./ (d_pre - d_nex); 
CCT = 1.0e6 ./ (( r(Ip+1)' - r(Ip)' ).*p + rp);

end