Repérage des couleurs CIE XYZ

1. Fonctions colorimétriques

La commission internationale de l'éclairage (CIE) a défini en 1931 un système de 3 couleurs primaires virtuelles permettant de représenter toutes les couleurs visibles par trois composantes trichromatiques (X,Y,Z) positives ([1] [2]).

Les fonctions colorimétriques

x_{c} (λ), y_{c} (λ), z_{c} (λ)

représentent les composantes trichromatiques des radiations monochromatiques. Elles se déduisent des fonctions colorimétriques du système CIE 1931 RGB, lesquelles ont été obtenues expérimentalement avec les trois primaires monochromatiques de longueur d'onde 436, 546 et 700 nanomètres.

La fonction y_c(λ) est identique à la sensibilité spectrale visuelle relative. Elle vaut 1 pour la longueur d'onde du maximum de sensibilité visuelle, 555 nm.

Le fichier ciexyz31.txt contient les valeurs de ces fonctions tous les nanomètres de 360 à 830 nm.

Connaissant ces fonctions, il est possible de déterminer les composantes trichromatiques XYZ d'une lumière dont on connait la composition spectrale. Le module Python présenté ici effectue les différents calculs.

Le fichier CieXYZ.py définit une classe CIEXYZ. Le constructeur lit le fichier de données.

import re
from math import * 
class CIEXYZ:
    def __init__(self,filename):
        f=open(filename)
        self.Lambda=[]
        self.xc=[]
        self.yc=[]
        self.zc=[]
        self.S=[]
        self.Sxc=[]
        self.Syc=[]
        self.Szc=[]
        I=0
        for line in f.readlines():
            li = line.strip()
            t = re.split('[\s|,]+',li)
            self.Lambda.append(float(t[0]))
            self.xc.append(float(t[1]))
            self.yc.append(float(t[2]))
            self.zc.append(float(t[3]))
            self.S.append(1.0)
            self.Sxc.append(float(t[1]))
            self.Syc.append(float(t[2]))
            self.Szc.append(float(t[3]))
            I += float(t[2])
        self.K = 1.0/I
        f.close()

La fonction suivante fournit les valeurs des fonctions colorimétriques pour une longueur d'onde donnée, avec une interpolation linéaire des données de la table :

    def fColor(self,L):
        if ((L<360)|(L>829)):
            return [0,0,0]
        else:
            i=int(floor(L-360.0))            
            x=self.xc[i]+(L-float(i)-360.0)*(self.xc[i+1]-self.xc[i])
            y=self.yc[i]+(L-float(i)-360.0)*(self.yc[i+1]-self.yc[i])
            z=self.zc[i]+(L-float(i)-360.0)*(self.zc[i+1]-self.zc[i])
            return [x,y,z]

Représentation graphique des fonctions colorimétriques :

from pylab import *
from CieXYZ import CIEXYZ
cie=CIEXYZ("ciexyz31.txt")
clf()
xlabel('lambda (nm)')
plot(cie.Lambda,cie.xc,c='r',label='xc')
plot(cie.Lambda,cie.yc,c='g',label='yc')
plot(cie.Lambda,cie.zc,c='b',label='zc')
title(u'Fonctions colorimétriques')
legend()
grid(True)

plotA.pdf

print(cie.fColor(450))
--> [0.3362, 0.037999999999999999, 1.7721100000000001]

2. Composantes et coordonnées trichromatiques

Une lumière quelconque est représentée par ses composantes trichromatiques X,Y,Z. Afin de séparer la luminance des caractéristiques chromatiques, ont définit les coordonnées trichromatiques par :

\begin{matrix} x = \frac{X}{X + Y + Z} \\ y = \frac{Y}{X + Y + Z} \end{matrix}

x et y définissent la couleur et Y la luminance.

Inversement, si les coordonnées chromatiques x et y et la luminance Y sont connues, on calcule les composantes trichromatiques par :

\begin{matrix} X = \frac{x}{y} Y \\ Z = \frac{1 - x - y}{y} Y \end{matrix}

Les fonctions suivantes calculent les coordonnées trichromatiques x et y à partir des composantes trichromatiques X,Y,Z et inversement:

    def XYZ2xy(self,X,Y,Z):
        return [X/(X+Y+Z),Y/(X+Y+Z)]
    def xyY2XYZ(self,x,y,Y):
        return [x/y*Y,Y,(1-x-y)/y*Y]

La fonction suivante effectue le calcul des coordoonées chromatiques (x,y) des radiations monochromatiques :

    def xyColor(self):
        x=[]
        y=[]
        for k in range(len(self.xc)):
            xy=self.XYZ2xy(self.xc[k],self.yc[k],self.zc[k])
            x.append(xy[0])
            y.append(xy[1])
        return [x,y]

Le diagramme de chromaticité consiste à représenter les coordonnées (x,y) dans un repère orthogonal. En particulier, on peut représenter sur ce diagramme les couleurs monochromatiques (spectrum locus) :

xy=cie.xyColor()
clf()
xlabel('x')
ylabel('y')
plot(xy[0],xy[1],c='k')
title(u"Diagramme de chromaticité")
grid(True)

plotB.pdf

3. Les illuminants

La couleur des objets, par réflexion ou par transmission, dépend de la source de lumière utilisée.

3.a. Corps noir

Le corps noir est un illuminant théorique dont la luminance spectrale est donnée par la loi de Planck.

La répartition spectrale relative S(λ) d'un illuminant a par convention la valeur 100 pour λ = 560 nm. Pour le corps noir de température T, on aura donc :

S (λ) = 100 (\frac{560}{λ})^{5} (\frac{exp (\frac{A}{560 T}) - 1}{exp (\frac{A}{λ T}) - 1})

La constante A est :

A = 1.4388 \cdot 1 0^{7} n m \cdot K

La lumière solaire directe en milieu de journée est approximativement celle du corps noir à 5800 K. La lumière diffusée par temps couvert a une température plus élevée. Ces variations de la température de couleur en fonction de la météo et de l'heure sont bien connues des photographes.

La fonction suivante calcule la répartition spectrale du corps noir :

    def corpsNoir(self,L,T):
        return 100*pow(560/L,5)*(exp(2.569e4/T)-1)/(exp(1.4388e7/(L*T))-1)

3.b. Illuminants standard

L'illuminant standard A est défini par le rayonnement du corps noir à 2855.5 K. Il est réalisé par une lampe à filament de tungstène dont le courant est convenablement régulé.

La CIE a défini des illuminants plus réalistes que le corps noir pour la lumière du jour (Daylight). Par exemple, l'illuminant D65 correspond à la lumière du jour par temps couvert. L'illuminant D50 est utilisé dans les arts graphiques (imprimerie, photographie, etc). On peut citer aussi les illuminants de la famille F (F1 à F12) qui correspondent aux lampes à fluorescence.

Les fonctions suivantes lisent les répartitions spectrales des illuminants A et D65 à partir d'un fichier de données (IlluminantA.txt et IlluminantD65.txt):

    def readA(self,filename):
        f=open(filename)
        self.illumALambda=[]
        self.illumASpectre=[]
        for line in f.readlines():
            li = line.strip()
            t = re.split('[\s|,]+',li)
            self.illumALambda.append(float(t[0]))
            self.illumASpectre.append(float(t[1]))
        f.close()
    def readD65(self,filename):
        f=open(filename)
        self.illumD65Lambda=[]
        self.illumD65Spectre=[]
        for line in f.readlines():
            li = line.strip()
            t = re.split('[\s|,]+',li)
            self.illumD65Lambda.append(float(t[0]))
            self.illumD65Spectre.append(float(t[1]))
        f.close()

Voyons le tracé de ces répartitions spectrales :

cie.readA("IlluminantA.txt")
cie.readD65("IlluminantD65.txt")
clf()
xlabel('lambda (nm)')
ylabel('S')
plot(cie.illumALambda,cie.illumASpectre,c='r',label='A')
plot(cie.illumD65Lambda,cie.illumD65Spectre,c='b',label='D65')
L=linspace(300,800,200)
corpsNoir6504=zeros(200)
for p in range(200):
    corpsNoir6504[p]=cie.corpsNoir(L[p],6504)
plot(L,corpsNoir6504,c='g',label='T=6504 K')
legend()
title('Illuminants')
grid(True)

plotC.pdf

La température proximale de l'illuminant D65 est 6504 K : c'est la température du corps noir dont le rayonnement s'approche le plus.

La répartition spectrale de l'illuminant D65 est obtenue à partir de mesures de la lumière du jour effectuées par intervalle de 10 nm. La table fournit les valeurs interpolées tous les 1 nm. La fonction suivante permet de calculer la valeur de S(λ) pour une longueur d'onde quelconque :

    def SpectreD65(self,L):
        L=floor(L)
        if ((L<300)|(L>830)):
            return 0
        else:
            i=int(L-300)
            return self.illumD65Spectre[i]

4. Calcul des composantes trichromatiques

Un objet est caractérisé par son facteur spectral de luminance R(λ). Il s'agit du coefficient de réflexion ou de transmission selon le cas. Pour une observation directe de l'illuminant, on posera simplement R(λ) = 1.

Les composantes trichromatiques de cet objet éclairé par l'illuminant de répartition spectrale S(λ) sont :

\begin{matrix} X = K \int S (λ) x_{c} (λ) R (λ) d λ \\ Y = K \int S (λ) y_{c} (λ) R (λ) d λ \\ Z = K \int S (λ) z_{c} (λ) R (λ) d λ \\ K = \frac{1}{\int S (λ) y_{c} (λ) d λ} \end{matrix}

Par convention, le facteur de luminance de l'illuminant observé directement est Y = 1. Sachant que R <= 1, le facteur de luminance d'un objet est inférieur ou égal à 1.

Pour l'intégration numérique, on utilise la méthode des trapèzes avec un pas de Δλ = 1 nm, de 360 à 830 nm.

Par exemple, le calcul de X se fait par la somme suivante :

X = K \sum_{λ = 360}^{830} S (λ) x_{c} (λ) Δ λ

Compte tenu de la normalisation par le facteur K, l'intervalle Δλ peut être omis dans le calcul.

La fonction suivante permet de sélectionner l'illuminant utilisé pour les calculs ultérieurs et de calculer la constante de normalisation K. Les produits de la répartition spectrale par les fonctions colorimétriques sont également calculés. La répartition spectrale de l'illuminant est enregistrée dans un tableau de 360 à 830 nm. Si le corps noir est choisi (CN), la température doit être fournie. Remarque : si cette fonction n'est pas appelée, l'illuminant par défaut est E, défini par S(λ)=1.

    def setIlluminant(self,nom,T):
        if nom=="CN":
            for k in range(471):
                L=float(360+k)
                self.S[k] = self.corpsNoir(L,T)
        elif nom=="D65":
            for k in range(471):
                self.S[k] = self.illumD65Spectre[60+k]
        elif nom=="A":
            for k in range(471):
                L=float(360+k)
                self.S[k] = self.corpsNoir(L,2855.5)
        elif nom=='E':
            for k in range(471):
                self.S[k] = 1
        I = 0
        for k in range(471):
            self.Sxc[k] = self.S[k]*self.xc[k]
            self.Syc[k] = self.S[k]*self.yc[k]
            self.Szc[k] = self.S[k]*self.zc[k]
            I += self.Syc[k]
        self.K = 1.0/I

La fonction suivante permet de définir l'illuminant à partir d'une fonction spectrale donnée en argument.

    def setFonctionIlluminant(self,f):
        for k in range(471):
            L=float(360+k)
            self.S[k] = f(L)
        I = 0
        for k in range(471):
            self.Sxc[k] = self.S[k]*self.xc[k]
            self.Syc[k] = self.S[k]*self.yc[k]
            self.Szc[k] = self.S[k]*self.zc[k]
            I += self.Syc[k]
        self.K = 1.0/I

La fonction suivante calcule les composantes trichromatiques d'un objet dont le facteur spectrale de luminance est donné par une fonction :

    def spectralF2XYZ(self,Rf):
        X=0
        Y=0
        Z=0
        for k in range(471):
            L=float(360+k)
            R = Rf(L)
            X += self.Sxc[k]*R
            Y += self.Syc[k]*R
            Z += self.Szc[k]*R
        return[self.K*X,self.K*Y,self.K*Z]

La fonction suivante traite le cas particulier d'une lumière monochromatique (largeur de raie de 1 nm).

    def XYZMonochrome(self,L):
        if (L<360)|(L>830):
            return [0,0,0]
        k=int(floor(L-360))
        X = self.Sxc[k]
        Y = self.Syc[k]
        Z = self.Szc[k]
        return[self.K*X,self.K*Y,self.K*Z]

Exemple : composantes trichromatiques puis coordonnées trichromatiques du corps noir à 5000 K :

cie.setIlluminant("CN",5000)
def Rf(L):
    return 1
XYZ=cie.spectralF2XYZ(Rf)
xy5000=cie.XYZ2xy(XYZ[0],XYZ[1],XYZ[2])

print(XYZ)
--> [0.98149534412352402, 1.0, 0.86256604654545821]

print(xy5000)
--> [0.3451034310805281, 0.35160985036429865]

On place les illuminants sur le diagramme de chromaticité :

clf()
xlabel('x')
ylabel('y')
xy=cie.xyColor()
plot(xy[0],xy[1],c='k',label='Spectrum locus')
plot([xy5000[0]],[xy5000[1]],c='r',marker='o',label='T=5000 K')
cie.setIlluminant("D65",0)
XYZ=cie.spectralF2XYZ(Rf)
xyD65=cie.XYZ2xy(XYZ[0],XYZ[1],XYZ[2])
plot([xyD65[0]],[xyD65[1]],c='b',marker='o',label='D65')
cie.setIlluminant("A",0)
XYZ=cie.spectralF2XYZ(Rf)
xyA=cie.XYZ2xy(XYZ[0],XYZ[1],XYZ[2])
plot([xyA[0]],[xyA[1]],c='g',marker='o',label='A')
legend()
grid(True)

plotD.pdf

5. Grandeurs énergétiques et lumineuses

La fonction colorimétrique y_c(λ) correspond à la sensibilité spectrale relative du système visuel. Elle permet donc d'effectuer les conversions entre grandeurs énergétiques et grandeurs visuelles associées. Les deux fonctions suivantes effectuent ces conversions :

    def ener2lum(self,e,L):
        xyz=self.fColor(L)
        return e*683.0*xyz[1]
    def lum2ener(self,lu,L):
        xyz=self.fColor(L)
        return lu/(683.0*xyz[1])

6. Retour sur le le système CIE RGB 1931

la fonction suivante calcule les fonctions colorimétriques expérimentales RGB qui ont servi à la détermination des fonctions colorimétriques XYZ.

    def rgbColorim(self):
        r=[]
        g=[]
        b=[]
        for k in range(len(self.xc)):
            r.append(0.41846*self.xc[k]-0.15866*self.yc[k]-0.08283*self.zc[k])
            g.append(-0.09117*self.xc[k]+0.25242*self.yc[k]+0.01571*self.zc[k])
            b.append(0.00092*self.xc[k]-0.00255*self.yc[k]+0.17860*self.zc[k])
        return [r,g,b]

Références

[1] R. Sève, Physique de la couleur, (Masson, 1996)

[2] D. Dupont, D. Steen, J.F. Bloch, Métrique de la couleur, La couleur : lumière, vision et matériaux, sous la direction de M. Elias et J. Lafait (Belin, 2006)