Data | Fundamental Data |
Return to | Internet Data Base | COLORING Index | or visit | Bill Dawes' Home Page |
Contributed by David Heggie, HeggieD@scot.hw.ac.uk
Scottish College of Textiles, Heriot Watt University
I asked David Heggie if he'd be willing to contribute code, along with his description of the formulas for CIE94*. He sent the following, which calculates CIE94, as well as CIELab, BFD, and CMC differences!
If you need descriptions of CIE76 dE*, or CMC, descriptions can be found in Alan Roberts and Adrian Ford's FAQ.
Bill:
I've included the program at the end of this mailing - bear in mind that it's a program (well, part of a program) I cobbled together to help me out with part of my work, so it's certainly *not* user-friendly! It's written in Turbo Pascal for DOS, but I should imagine with a couple of changes it should be OK elsewhere. You could easily modify it to enable user input - as it stands, you have to modify the program to input the sample / standard L* a* b* values - not exactly easy if it's already compiled, but as I said, it was never intended for this purpose. The actual functions which do the maths are OK 'though. Hope it's some use.
------------------DIFS.PAS-------------------------------
{
DIFFS.PAS (C)1994 David Heggie, SCOT
root - calculates yth root of x.
geth - calculates hue angle (h) for sample a*, b*
cielab - gets a CIELAB difference between L*1, a*1, b*1 and L*2, a*2, b*2
cielab94 - gets CIE94 difference between L*1, a*1, b*1 and L*2, a*2, b*2
bfd - gets BFD(1:1) difference between L*1, a*1, b*1 and L*2, a*2, b*2
cmc - gets CMC(1:1) difference between L*1, a*1, b*1 and L*2, a*2, b*2
}
program diffs;
uses crt;
function root(x, y : double) : double;
begin
root := exp((1/y)*LN(abs(x)));
end;
function cielab(l1,a1,b1,l2,a2,b2 : double) : double;
var
dl, da, db : double;
begin
dl := l2-l1;
da := a2-a1;
db := b2-b1;
cielab := sqrt(sqr(dl)+sqr(da)+sqr(db))
end;
function geth(a, b : double) : double;
begin
if (a=0) and (b=0) then geth := 0;
if (a=0) and (b>0) then geth := 90;
if (a=0) and (b<0) then geth :="270;" if (b="0)" and (a<0 ) then geth :="180;" if (b="0)" and (a>0) then geth := 0;
if (a0) and (b0) then
begin
geth := 360+(arctan(b/a)*180/pi);
if a<0 then geth :="180+(arctan(b/a)*180/pi);" if (a>0) and (b>0) then geth:= arctan(b/a)*180/pi;
end;
end;
function cielab94(l1,a1,b1,l2,a2,b2 : double) : double;
var
c1,c2,de,dl,dc,dh,sl,sc,sh,kl,kc,kh : double;
begin
sl := 1;
kl := 1;
kc := 1;
kh := 1;
c1 := sqrt(sqr(a1)+sqr(b1));
c2 := sqrt(sqr(a2)+sqr(b2));
dl := l2-l1;
dc := c2-c1;
de := cielab(l1,a1,b1,l2,a2,b2);
if sqr(de)>(sqr(dl)+sqr(dc)) then dh := sqrt(sqr(de)-sqr(dl)-sqr(dc)) else
dh := 0;
sc := 1 + 0.045*c1;
sh := 1 + 0.015*c1;
cielab94 := sqrt(sqr(dl/(kl*sl)) + sqr(dc/(kc*sc)) + sqr(dh/(kh*sh)));
end;
function bfd(l1,a1,b1,l2,a2,b2 : double) : double;
const
loge = 0.434294481;
var
c1,c2,h1,h2,lbfd1,lbfd2,avec,aveh,de,deltal,
deltac,deltah,dc,t,g,dh,rh,rc,rt,deltae,yr,yt,xt,zt : double;
begin
if l1 > 7.996969 then
yt := (sqr((l1+16)/116)*((l1+16)/116))*100
else
yt := 100 * (l1 / 903.3);
lbfd1:=54.6*(loge*(ln(yt+1.5)))-9.6;
if l2 > 7.996969 then
yt := (sqr((l2+16)/116)*((l2+16)/116))*100
else
yt := 100 * (l2 / 903.3);
lbfd2:=54.6*(loge*(ln(yt+1.5)))-9.6;
c1:=sqrt(sqr(a1)+sqr(b1));
c2:=sqrt(sqr(a2)+sqr(b2));
deltal:=lbfd2-lbfd1;
deltac:=c2-c1;
de:=cielab(l1,a1,b1,l2,a2,b2);
if sqr(de)>(sqr(l2-l1)+sqr(deltac)) then
deltah:=sqrt(sqr(de)-sqr(l2-l1)-sqr(deltac))
else
deltah:=0;
avec:=(c1+c2)/2;
h1:=geth(a1,b1);
h2 := geth(a2,b2);
aveh:=(h1+h2)/2;
dc:=0.035*avec/(1+0.00365*avec)+0.521;
g:=sqrt(sqr(sqr(avec))/(sqr(sqr(avec))+14000));
t:=0.627+(0.055*cos((aveh-254)/(180/pi))-
0.040*cos((2*aveh-136)/(180/pi))+
0.070*cos((3*aveh-31)/(180/pi))+
0.049*cos((4*aveh+114)/(180/pi))-
0.015*cos((5*aveh-103)/(180/pi)));
dh:=dc*(g*t+1-g);
rh:=-0.260*cos((aveh-308)/(180/pi))-
0.379*cos((2*aveh-160)/(180/pi))-
0.636*cos((3*aveh+254)/(180/pi))+
0.226*cos((4*aveh+140)/(180/pi))-
0.194*cos((5*aveh+280)/(180/pi));
rc:=sqrt((avec*avec*avec*avec*avec*avec)/((avec*avec*avec*avec*avec*avec)+70
000000));
rt:=rh*rc;
bfd:=sqrt(sqr(deltal)+sqr(deltac/dc)+sqr(deltah/dh)+(rt*(deltac/dc)*(deltah/
dh)));
end;
function cmc(l1,a1,b1,l2,a2,b2 : double) : double;
var
c1,c2,h1,h2,avec,aveh,de,dl,dc,dh,sl,sc,sh,t,f : double;
begin
c1:=sqrt(sqr(a1)+sqr(b1));
c2:=sqrt(sqr(a2)+sqr(b2));
dl:=l2-l1;
dc:=c2-c1;
de:=cielab(l1,a1,b1,l2,a2,b2);
if sqr(de)>(sqr(dl)+sqr(dc)) then dh:=sqrt(sqr(de)-sqr(dl)-sqr(dc)) else
dh:=0;
h1:=geth(a1,b1);
h2:=geth(a2,b2);
if (h1>164) and (h1<345) then t:="0.56+abs(0.2*cos(((h1+168)/(180/pi))))" else t:="0.36+abs(0.4*cos(((h1+35)/(180/pi))));" sc:="0.0638*c1/(1+0.0131*c1)+0.638;" sl:="0.040975*l1/(1+0.01765*l1);" if l1<16 then sl:="0.511;" f:="sqrt((c1*c1*c1*c1)/((c1*c1*c1*c1)+1900));" sh:="sc*(t*f+1-f);" cmc:="sqrt(sqr(dl/sl)+sqr(dc/sc)+sqr(dh/sh));" end; var l, a, b,l2,a2,b2, dl, da, db, de : double; begin l :="50.00;" a :="20.00;" b :="20.00;" l2 :="51;" a2 :="21;" b2 :="21;" writeln; writeln ('cielab : ',cielab(l,a,b,l2,a2,b2):7:2); writeln ('cmc : ',cmc(l,a,b,l2,a2,b2):7:2); writeln ('bfd : ',bfd(l,a,b,l2,a2,b2):7:2); writeln ('cie94 : ',cielab94(l,a,b,l2,a2,b2):7:2); writeln; end.
Tagged December 19, 1995