Introduction¶
Le cours d'aujourd'hui reprend, avec des détails supplémentaires, un cours d'« étude de cas » que j'ai donné à des M1 d'ingénierie mathématique (où les codes sont en R au lieu de Python).
Nous allons nous occuper aujourd'hui avec de vraies données recueillies dans une étude de la dynamique des ions calcium $Ca^{2+}$ dans des neurones de l'hypothalamus de souris. Ces neurones libèrent des hormones peptidiques dérivées de la Pro-opiomélanocortine (POMC) dont certains sont fortement impliqués dans la régulation de l'appétit. Comme les problèmes d'obésité deviennent préoccupant dans nos sociétés, les physiologistes s'intéressent de plus en plus à ces problèmes -- comprenez qu'ils peuvent obtenir de l'argent (essentiellement de sources publics) pour faire « tourner » leurs labos quand ils travaillent sur ces problèmes.
Mes collègues du laboratoire de Peter Kloppenburg à l'Université de Cologne (Allemagne) travaillent sur ces questions en collaboration avec les généticiens du laboratoire de Jens Brüning (même université). Ils étudient l'effet du régime alimentaire imposé aux souris (basses calories, niveau normal, hautes calories) sur l'activité des neurones synthétisant la POMC. Plus particulièrement ils étudient l'effet du régime sur la dynamique des ions calcium -- que nous savons par ailleurs avoir une influence prépondérante sur l'activité de ces neurones. Les données que nous allons utiliser ici ont été enregistrées par Andreas Pippow lors de sa thèse dans le labo de Peter Kloppenburg. Nous allons utilser deux jeux de données :
- un jeu de données de calibration du capteur CCD utilisé dans les expériences ; ces données sont disponibles en deux formats :
- un format propre à
numpyà l'adresse : http://xtof.disque.math.cnrs.fr/data/CCD_calibration.npz ; - un format plus général, HDF5 -- mais qui nécessite l'installation du paquet h5py -- à l'adresse : http://xtof.disque.math.cnrs.fr/data/CCD_calibration.hdf5 ;
- un format propre à
- un jeu de données de mesure de fluorescence liées au calcium dans un neurone synthétisant la POMC ; ce jeu est disponible en trois formats :
- un format
textà l'adresse : http://xtof.disque.math.cnrs.fr/data/POMC.txt ; - un format propre à
numpyà l'adresse : http://xtof.disque.math.cnrs.fr/data/Data_POMC.npz ; - au format
HDF5à l'adresse : http://xtof.disque.math.cnrs.fr/data/Data_POMC.hdf5.
- un format
Si vous souhaitez en savoir plus sur les raisons qui poussent les physiologistes à étudier la dynamique des ions calcium en général, consulter mon cours sur le sujet.
Chargement des données « POMC » dans Python¶
Nous procédons de la « manière habituelle » après avoir démarré IPython :
%pylab
from __future__ import print_function, division, unicode_literals, absolute_import
Si nous souhaitons avoir les figures dans notre notebook IPython nous utilisons aussi (avec une commande additionnelle pour améliorer la résolution des graphiques que j'ai trouvée sur un excellent site) :
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
Nous allons maintenant voir successivement comment charger les données dans les trois formats disponibles.
Lecture des données au format text¶
from urllib2 import urlopen
u = urlopen("http://xtof.disque.math.cnrs.fr/data/POMC.txt")
Les 6 premières lignes du fichier contiennent une description des données :
u.readline()
u.readline()
u.readline()
u.readline()
u.readline()
u.readline()
Les 168 lignes suivantes contiennent le temps (en secondes) auquel chacune des 168 expositions a été effectuée. Nous stockons ces temps dans un vecteur :
temps_POMC = np.zeros((168,))
for i in range(168):
temps_POMC[i] = float(u.readline())
Le capteur CCD est une matrice de 60 x 80 pixels. Les 60 x 80 x 168 lignes suivantes contiennent les lectures en sortie de capteur de pour chaque pixel lors de chacune des 168 expositions. Nous les stockons dans un array à trois dimensions :
pile_POMC = np.zeros((60,80,168),dtype=int)
for k in range(168):
for j in range(80):
for i in range(60):
pile_POMC[i,j,k] = int(u.readline())
Une fois la lecture finie, nous fermons le « fichier » :
u.close()
À ce stade, c'est une bonne idée de regarder les données pour s'assurer qu'aucune grosse bêtise n'a été commise lors de leur chargement. Une façon rapide de faire est de calculer la moyenne (au cours du temps) du signal pour chaque pixel et de regarder ce que ça donne.
La moyenne par pixel s'obtient simplement avec la méthode mean :
pile_POMC_moyenne = pile_POMC.mean(2)
On visualise alors avec la fonction imshow de pyplot :
imshow(transpose(pile_POMC_moyenne),origin='lower')
set_cmap('gray')
colorbar()
Lecture des données au format npz¶
La version npz des données créée avec la fonction savez_compressed de numpy doit d'abord être téléchargée, ce qui se fait simplement avec :
import urllib
urllib.urlretrieve("http://xtof.disque.math.cnrs.fr/data/Data_POMC.npz","Data_POMC.npz")
Le fichier Data_POMC.npz contient deux arrays, le premier avec le temps des expositions, le second avec la pile d'images. On charge les données de la façon suivante :
POMC = load("Data_POMC.npz")
temps_POMC_npz = POMC['arr_0']
pile_POMC_npz = POMC['arr_1']
POMC.close()
Ici une rapide manipulation nous permet de vérifier que nous avons chargé les mêmes données :
pile_POMC_npz_moyenne = pile_POMC_npz.mean(2)
imshow(transpose(pile_POMC_npz_moyenne),origin='lower')
colorbar()
Lecture des données au format hdf5¶
Une troisième version des données est disponible en format HDF5 ; un format de stockage hiérarchique -- l'organisation interne du fichier est similaire à celle des fichiers et répertoires sur votre ordinateur -- qui permet d'annoter les jeux de données -- il est ainsi possible de sauvegarder les informations relatives à l'acquisition des données avec les données elles-mêmes. Un module Python, le module h5py permet de manipuler ces fichiers sans peine.
Nous procédons comme dans le cas précédent en téléchargeant les données :
urllib.urlretrieve("http://xtof.disque.math.cnrs.fr/data/Data_POMC.hdf5","Data_POMC.hdf5")
Nous continuons en chargeant le module puis nous ouvrons le fichier en mode lecture :
import h5py
pomc = h5py.File("Data_POMC.hdf5","r")
Pour obtenir la liste des jeux de données dataset et des groups (équivalents aux répertoires d'un système de fichier), on emploie la fonction list :
list(pomc)
Nous retrouvons le vecteur des temps, time, et la pile d'images, stack. Pour extraire ces deux jeux de données nous employons :
temps_POMC_hdf5 = pomc['time'][...]
pile_POMC_hdf5 = pomc['stack'][...]
Une fois les données extraites, nous fermons le fichier :
pomc.close()
Comme précédemment, nous vérifions que l'importation s'est « bien » passée avec :
pile_POMC_hdf5_moyenne = pile_POMC_hdf5.mean(2)
imshow(transpose(pile_POMC_hdf5_moyenne),origin='lower')
colorbar()
Représentation des données¶
Le plus souvent, les physiologistes vont montrer des piles d'images sous forme de films. C'est possible de le faire dans Python et si vous souhaitez le faire vous trouverez sur le oueb les indications nécessaires. Je préfère néanmoins une représentation qui montre sous forme matricielle, les éléments de la matrice correspondant aux pixels du capteur CCD, le décours temporel de la fluorescence. Pour obtenir ce type de représentation, nous allons définir une fonction :
def plotSignal(stack,lw=1):
n_x, n_y, n_t = stack.shape
amp_min = np.min(stack)
amp_max = np.max(stack)
amp_diff = np.ptp(stack)
x_domain = np.arange(n_t)/n_t
y_domain = (0,n_y)
for r_idx in range(n_x):
for c_idx in range(n_y):
y_min = n_x - r_idx - 1
sig = stack[r_idx,c_idx,:]
Y = (sig-amp_min)/amp_diff + y_min
X = x_domain + c_idx
plt.plot(X,Y,lw=lw,color='black')
plt.axis('off')
Nous appelons notre fonction sur notre pile en regardant ce qui se passe du coté du soma :
plotSignal(pile_POMC[20:41,30:51,:],lw=0.5)
Ce graphe est nettement plus lisible en mode interactif (sans la commande inline).
Quelques références sur la contruction des graphes¶
Comme discuté en cours, voici quelques « bons bouquins » qui discutent en détails de la représentation graphique des données :
- Sémiologie graphique. Les diagrammes, les réseaux, les cartes de Jacques Bertin : comme c'est écrit par un Français, c'est forcément plutôt théorique ;
- les bouquins d'Edward Tufte : Visual Display of Quantitative Information, Envisioning Information, Visual Explanations et Beautiful Evidence ; bien écrits et très beaux même si le style de l'auteur qui manque parfois de modestie m'irrite un peu ;
- les bouquins de William Cleveland : Visualizing Data (si vous devez n'en posséder qu'un, je vous suggère celui-ci) et The Elements of Graphing Data.
Vous pouvez aussi consulter avec profit le site de Michael Friendly. Enfin, il y a un module très intéressant en cours de développement pour Python : Bokeh -- je ne vous en ai pas parlé dans le cours car il faut installer un module supplémentaire, mais il vient avec la distribution Anaconda.
La variabilité inhérente aux données de fluorescence¶
L'examen du décours temporel du signal d'un des pixels centraux nous donne :
plot(temps_POMC,pile_POMC[29,39,:],lw=1)
xlabel("Temps (s)")
ylabel("Compte (ADU)")
xlim([525,550])
grid()
Que cherche-t-on ?¶
À partir des données que nous venons de voir, nous pourrions vouloir estimer des paramètres comme :
- l'amplitude au pic ;
- la constante de temps de retour à l'« équilibre » ;
- le niveau de base ;
- le décours temporel complet -- c.-à-d. une fonction à strictement parler.
Nous pourrions aussi considérer un modèle liant la dynamique calcique -- le décours temporel de la concentration de calcium libre dans la cellule -- à l'intensité de fluorescence : $$\frac{\mathrm{d}Ca_t}{\mathrm{dt}} \left(1 + \kappa_{F}(Ca_t) + \kappa_{E}(Ca_t) \right) + \frac{j(Ca_t)}{v} = 0 \, , $$ où $Ca_t$ représente la $[Ca^{2+}]_{libre}$ au temps t, $v$ est le volume du neurite -- à l'intérieur duquel les effets de diffusion peuvent être négligés -- et $$j(Ca_t) \equiv \gamma (Ca_t - Ca_{station.}) \, ,$$ modélise l'extrusion du calcium -- $Ca_{station.}$ est la $[Ca^{2+}]_{libre}$ stationnaire -- $$\kappa_{F}(Ca_t) \equiv \frac{F_{total} \, K_{F}}{(K_{F} + Ca_t)^2} \quad \mathrm{et} \quad \kappa_{E}(Ca_t) \equiv \frac{E_{total} \, K_{E}}{(K_{E} + Ca_t)^2} \, ,$$ où $F$ désigne le fluorophore et $E$ le(s) tampon(s) endogène(s). Quand nous travaillons avec ce type de modèle, nous supposons que les paramètres du fluorophore (le Fura) : $F_{total}$ et $K_F$ ont été calibrés ; nous pourrions vouloir estimer :
- le paramètre d'extrusion : $\gamma$
- les paramètres du tampon endogène : $E_{total}$ et $K_E$.
en utilisant une équation reliant l'intensité de fluorescence au calcium : $$Ca_t = K_{F} \, \frac{S_t - S_{min}}{S_{max} - S_t} \, ,$$ où $S_t$ est la fluorescence mesurée au temps $t$, $S_{min}$ et $S_{max}$ sont des paramètres calibrés correspondant à la fluorescence en absence de calcium et en calcium « saturant ».
- la variabilité de notre signal implique que nos paramètres estimés vont fluctuer d'une réplication à l'autre;
- formellement, nos paramètres estimés sont modélisés par des variables aléatoires et il n'est pas suffisant de résumer une variable aléatoire par un nombre ;
- si nous ne pouvons pas obtenir la loi d'échantillonage complète de nos paramètres, nous voulons au moins fournir un intervalle à l'intérieur duquel la vraie valeur du paramètre se trouve avec une probabilité donnée ;
- dit autrement : une analyse sans intervalle de confiance n'est pas une analyse puisque ses résultats ne sont pas reproductibles -- si j'écris que ma constante de temps est 25,76 ms, la probabilité d'obtenir 25,76 ms en répliquant l'expérience est essentiellement 0 ; si j'écris que la vraie constante de temps a une probabilité de 0,95 de se trouver dans l'intervalle [24,0 ; 26,5], je peux comparer ce résultat avec celui obtenu lors d'une réplication.
L'importance de la prise en compte de la variabilité¶
Considérons un modèle de génération de données très simple : $$Y_i \sim \mathcal{P}(f_i)\, , \quad i=0,1,\ldots,K \; ,$$ où $\mathcal{P}(f_i)$ désigne une loi de Poisson de paramètre $f_i$ : $$\mathrm{Pr}\{Y_i = n\} = \frac{(f_i)^n}{n!} \exp (-f_i)\, , \quad \mathrm{pour} \quad n=0,1,2,\ldots $$ et $$f_i = f(\delta i| f_{\infty}, \Delta, \beta) = f_{\infty} + \Delta \, \exp (- \beta \, \delta i)\; ,$$ $\delta$ est un pas de temps et $f_{\infty}$, $\Delta$ et $\beta$ sont les paramètres du modèle.
Nous illustrons ce modèle avec une petite simulation en commençant par définir les paramètres de celle-ci :
beta_vrai = 1.0
f_infini = 100
Delta = 900
X = linspace(0,5*beta_vrai,51)
Theo = Delta*exp(-X*beta_vrai)+f_infini
Nous simulons 51 observations de 51 v.a. de Poisson dont les paramètres sont les 51 valeurs contenues dans le vecteur Theo :
np.random.seed(20061001)
Echantillon = np.random.poisson(Theo)
Sur une figure ça donne :
plot(X,Echantillon,'o')
xlabel("Temps (s)")
ylabel("Observation")
plot(X,Theo,'r')
plot(X[[3,30]],Theo[[3,30]],'sk')
plot([X[3],X[3]],[0,Theo[3]],'--k')
plot([0,X[3]],[Theo[3],Theo[3]],'--k')
plot([X[30],X[30]],[0,Theo[30]],'--k')
plot([0,X[30]],[Theo[30],Theo[30]],'--k')
text(0.1,730,r'$y_1$')
text(1.5,110,r'$y_2$')
Exemple de données simulées suivant le modèle précédent. Nous allons supposer à présent que $f_{\infty}$ et $\Delta$ sont connus et que $(t_1,y_1)$ et $(t_2,y_2)$ sont donnés. Nous souhaitons estimer $\beta$.
Deux estimateurs¶
Nous allons considérer et comparer deux estimateurs de $\beta$ :
- l'estimateur classique des moindres carrés : $$\tilde{\beta} = \arg \min \tilde{L}(\beta) \; ,$$ où $$\tilde{L}(\beta) = \sum_j \quad \big( y_j - f(t_j \mid \beta) \big)^2 \; ;$$
- l'estimateur des moindres carrés appliqué à la racine carrée des observations : $$\hat{\beta} = \arg \min \hat{L}(\beta) \; ,$$ où $$\hat{L}(\beta) = \sum_j \quad \big( y_j^{1/2} - f(t_j \mid \beta)^{1/2} \big)^2 \; .$$
Nous effectuons une étude empirique comme suit :
- nous simulons $10^5$ expériences telles que : $$(Y_1,Y_2) \sim \big(\mathcal{P}(f(0.3|\beta_0), \mathcal{P}(f(3|\beta_0)\big) \; ,$$ avec $\beta_0 = 1$ ;
- pour chaque paire simulée, $(y_1,y_2)^{[k]}$ ($k=1,\ldots,10^5$), nous minimisons $\tilde{L}(\beta)$ et $\hat{L}(\beta)$ pour avoir : $(\tilde{\beta}^{[k]},\hat{\beta}^{[k]})$ ;
- nous construisons les histogrammes de $\tilde{\beta}^{[k]}$ and $\hat{\beta}^{[k]}$.
Définition des fonctions requises par la méthode de Newton pour le premier estimateur¶
Nous souhaitons minimiser : $$\tilde{\beta} = \arg \min \tilde{L}(\beta) \; ,$$ où $$\tilde{L}(\beta) = \big( y_1 - f(t_1 \mid \beta) \big)^2 + \big( y_2 - f(t_2 \mid \beta) \big)^2 \; .$$
C'est-à-dire que nous cherchons une racine de la dérivée, dérivée que nous pouvons calculer à la main, où avec SymPy :
import sympy as sy
sy.init_printing()
x_1,x_2,y_1,y_2,Delta,beta,f_infini = sy.symbols('x_1,x_2,y_1,y_2,Delta,beta,f_infini',real=True)
L_tilde = (y_1 - Delta*sy.exp(-beta*x_1)-f_infini)**2 + (y_2 - Delta*sy.exp(-beta*x_2)-f_infini)**2
L_tilde
G_tilde = sy.diff(L_tilde,beta)
G_tilde
Une version directement utilisable pour une définition de fonction (numérique) de Python est obtenue avec :
print(G_tilde)
G_prime_tilde = sy.diff(G_tilde,beta)
G_prime_tilde
print(G_prime_tilde)
La définition de la fonction constructrice des fonctions nécessaires à la méthode de Newton suit :
def fait_g_dg_tilde(x_1,y_1,x_2,y_2,Delta=900,f_infini=100):
def g(beta):
return 2*Delta*x_1*(-Delta*exp(-beta*x_1) - f_infini + y_1)*exp(-beta*x_1) + 2*Delta*x_2*(-Delta*exp(-beta*x_2) - f_infini + y_2)*exp(-beta*x_2)
def dg(beta):
return 2*Delta**2*x_1**2*exp(-2*beta*x_1) + 2*Delta**2*x_2**2*exp(-2*beta*x_2) - 2*Delta*x_1**2*(-Delta*exp(-beta*x_1) - f_infini + y_1)*exp(-beta*x_1) - 2*Delta*x_2**2*(-Delta*exp(-beta*x_2) - f_infini + y_2)*exp(-beta*x_2)
return (g,dg)
Nous créons les fonctions dont nous avons besoin :
g_tilde, dg_tilde = fait_g_dg_tilde(X[3],Echantillon[3],X[30],Echantillon[30])
Nous réutilisons notre définition de la méthode de Newton du cours du 9 avril :
def newton(devinette_initiale,
cible,
cible_derivee,
tolerance=1e-6,
iter_max=100):
pos = devinette_initiale
valeur = cible(pos)
idx = 0
while idx <= iter_max and abs(valeur) > tolerance :
pos -= valeur/cible_derivee(pos)
idx += 1
valeur = cible(pos)
return (pos,valeur,idx)
newton(1.0,g_tilde,dg_tilde)
Définition des fonctions requises par la méthode de Newton pour le deuxième estimateur¶
Nous souhaitons minimiser : $$\hat{\beta} = \arg \min \hat{L}(\beta) \; ,$$ où $$\hat{L}(\beta) = \big( \sqrt{y_1} - \sqrt{f(t_1 \mid \beta)} \big)^2 + \big( \sqrt{y_2} - \sqrt{f(t_2 \mid \beta)} \big)^2 \; .$$
Nous utilisons SymPy car à la main c'est assez « lourd » :
L_chapeau = (sy.sqrt(y_1) - sy.sqrt(Delta*sy.exp(-beta*x_1) + f_infini))**2 + (sy.sqrt(y_2) - sy.sqrt(Delta*sy.exp(-beta*x_2) + f_infini))**2
G_chapeau = sy.diff(L_chapeau,beta)
print(G_chapeau)
G_prime_chapeau = sy.diff(G_chapeau,beta)
print(G_prime_chapeau)
def fait_g_dg_chapeau(x_1,y_1,x_2,y_2,Delta=900,f_infini=100):
def g(beta):
return Delta*x_1*(sqrt(y_1) - sqrt(Delta*exp(-beta*x_1) + f_infini))*exp(-beta*x_1)/sqrt(Delta*exp(-beta*x_1) + f_infini) + Delta*x_2*(sqrt(y_2) - sqrt(Delta*exp(-beta*x_2) + f_infini))*exp(-beta*x_2)/sqrt(Delta*exp(-beta*x_2) + f_infini)
def dg(beta):
return Delta**2*x_1**2*(sqrt(y_1) - sqrt(Delta*exp(-beta*x_1) + f_infini))*exp(-2*beta*x_1)/(2*(Delta*exp(-beta*x_1) + f_infini)**(3/2)) + Delta**2*x_1**2*exp(-2*beta*x_1)/(2*(Delta*exp(-beta*x_1) + f_infini)) + Delta**2*x_2**2*(sqrt(y_2) - sqrt(Delta*exp(-beta*x_2) + f_infini))*exp(-2*beta*x_2)/(2*(Delta*exp(-beta*x_2) + f_infini)**(3/2)) + Delta**2*x_2**2*exp(-2*beta*x_2)/(2*(Delta*exp(-beta*x_2) + f_infini)) - Delta*x_1**2*(sqrt(y_1) - sqrt(Delta*exp(-beta*x_1) + f_infini))*exp(-beta*x_1)/sqrt(Delta*exp(-beta*x_1) + f_infini) - Delta*x_2**2*(sqrt(y_2) - sqrt(Delta*exp(-beta*x_2) + f_infini))*exp(-beta*x_2)/sqrt(Delta*exp(-beta*x_2) + f_infini)
return (g,dg)
g_chapeau, dg_chapeau = fait_g_dg_chapeau(X[3],Echantillon[3],X[30],Echantillon[30])
newton(1.0,g_chapeau,dg_chapeau)
La simulation¶
n_rep = int(1e5)
beta_tilde = np.zeros((n_rep,3))
beta_chapeau = np.zeros((n_rep,3))
np.random.seed(20110928)
for rep_idx in range(n_rep):
Y = np.random.poisson(Theo[[3,30]])
g_tilde, dg_tilde = fait_g_dg_tilde(X[3],Y[0],X[30],Y[1])
beta_tilde[rep_idx,:] = newton(1.0,g_tilde,dg_tilde)
g_chapeau, dg_chapeau = fait_g_dg_chapeau(X[3],Y[0],X[30],Y[1])
beta_chapeau[rep_idx,:] = newton(1.0,g_chapeau,dg_chapeau)
Une fois la simulations -- avec les ajustements associés -- terminée, nous vérifions qu'aucun des ajustements ne s'est arrêté parce-que le nombre maximal d'itérations (100) a été atteint :
[any(beta_tilde[:,2] == 100),any(beta_chapeau[:,2] == 100)]
Nous pouvons maintenant construire les histogrammes des estimateurs (et ainsi obtenir une estimation de la loi d'échantillonage de ceux-ci). Les histogrammes sont construits avec la fonction hist de pylab. On peut aussi rajouter des versions théoriques de ces lois d'échantillonage -- la dérivation de ces versions théoriques se trouve dans Joucla et col. (2010).
def Ffct(beta): return 900 * exp(-X[[3,30]]*beta) + 100
def dFfct(beta): return -X[[3,30]]*900 * exp(-X[[3,30]]*beta)
sd0 = sqrt((sum(dFfct(1.0)**2*Ffct(1.0))/sum(dFfct(1.0)**2)**2))
sd1 = sqrt(1.0/sum(dFfct(1.0)**2/Ffct(1.0)))
beta_vecteur = linspace(0.6,1.6,501)
truc_inutile = hist([beta_tilde[:,0],beta_chapeau[:,0]],bins=50,normed=True)
plot(beta_vecteur,exp(-0.5*(beta_vecteur-1)**2/sd0**2)/sd0/sqrt(2*pi),color='blue',lw=2)
plot(beta_vecteur,exp(-0.5*(beta_vecteur-1)**2/sd1**2)/sd1/sqrt(2*pi),color='green',lw=2)
Chaque histogramme est construit avec 50 classes (bins). $\hat{\beta}$ (vert) est clairement meilleur que $\tilde{\beta}$ (bleu) car sa variance est plus faible.
Le « bruit » d'un capteur CCD¶
Le capteur CCD¶
Source: L. van Vliet et col. (1998) Digital Fluorescence Imaging Using Cooled CCD Array Cameras (figure 3).
Les sources de bruit d'un capteur CCD¶
- le bruit photonique ou « bruit de grenaille » (shot noise en anglais) vient du fait que mesurer une intensité de fluorescence, λ, implique un comptage de photons dont la loi est une loi de Poisson -- à part « changer les lois de la Physique », il n'y a rien à faire pour diminuer ce « bruit » -- : $$\mathrm{Pr}\{N=n\} = \frac{\lambda^n}{n!} \exp -\lambda\, , \quad n \, = \, 0,1,\ldots\, , \quad \lambda > 0\; ;$$
- le bruit thermique vient du fait que l'agitation thermique peut faire tomber des électrons dans les puys de potentiels, ce bruit suit aussi une loi de Poisson mais il peut être rendu négligeable en refroidissant la caméra ;
- le bruit de lecture vient de la conversion d'un nombre de photo-électrons en tension équivalente, il suit une loi normale dont la variance est indépendante de la moyenne (tant qu'on n'effectue pas de lectures à trop haute fréquence) ;
- le bruit de numérisation qui vient de la conversion d'une grandeur continue, la tension, en une valeur sur une grille, il est négligeable dès qu'on dispose de plus de 8 bit de codage.
Un modèle simple de CCD¶
Nous allons facilement pouvoir obtenir un modèle simple de CCD prenant en compte les deux sources de bruits principales. Mais pour y arriver nous avons besoin de réaliser qu'en général un grand nombre de photons est détecté, ce qui va nous permettre d'approcher la loi de Poisson par une loi normale dont la variance est égale à la moyenne : $$\mathrm{Pr}\{N=n\} = \frac{\lambda^n}{n!} \exp -\lambda \approx \mathcal{N}(\lambda,\lambda) \; .$$ Écrit plus « explicitement » : $$ N \approx \lambda + \sqrt{\lambda} \, \epsilon \; ,$$ où $\epsilon \sim \mathcal{N}(0,1)$ (suit une loi normale centrée, réduite). Le bruit de lecture s'ajoute après et suit une loi normale de moyenne nulle et de variance $\sigma_{L}^2$. On ajoute donc à notre variable aléatoire $N$ une nouvelle v.a. statistiquement indépendante $L \sim \mathcal{N}(0,\sigma_{L}^2)$ ce qui nous donne : $$M \equiv N+L \approx \lambda + \sqrt{\lambda+\sigma_{L}^2} \, \epsilon \; ,$$
$$M \equiv N+L \approx \lambda + \sqrt{\lambda+\sigma_{L}^2} \, \epsilon \; ,$$où j'ai utilisé le fait que la somme de deux v.a. normales indépendantes est une v.a. normale dont la moyenne est la somme des moyennes et la variance la somme des variances.
Enfin, la capacité des puys de photo-électrons de la caméra utilisée est de 35000 (c'est une valeur typique) et nous allons numériser ce nombre de photons, après lecture, sur 12 bit (soit 4096 niveaux), nous devons donc appliquer un « gain » $G$ /plus petit que un/ si nous voulons pouvoir représenter fidèlement (sans saturation) un puy presque plein. Nous observons ainsi : $$Y \equiv G \cdot M \approx G \, \lambda + \sqrt{G^2 \, (\lambda+\sigma_{L}^2)} \, \epsilon \; .$$
Calibration d'un capteur CCD¶
D'après ce qui précède, nos observations $Y$ suivent la loi : $$Y \sim G \, \lambda + \sqrt{G^2 \, (\lambda+\sigma_{L}^2)} \, \epsilon \; ,$$ où $G$ est le gain du capteur, $\sigma_{L}^2$ est la « variance de lecture » et $\epsilon$ suit une loi normale centrée réduite. La valeur de $G$ et $\sigma_{L}^2$ est une spécification du capteur, mais l'expérience prouve que les fabriquants ont tendance à être trop optimistes sur les performances de leurs produits. C'est donc une bonne idée de mesurer ces paramètres par des expériences de calibration. Cela va aussi nous permettre de vérifier que notre modèle très simple ci-dessus s'applique.
- notre problème devient donc : comment tester $Y \sim G \, \lambda + \sqrt{G^2 \, (\lambda+\sigma_{L}^2)} \, \epsilon$ ? ou comment faire varier $\lambda$ ;
- considérons un pixel de notre capteur « regardant » un volume fixe d'une lamelle fluorescente, nous avons deux façons de manipuler le paramètre λ :
- faire changer l'intensité $i_{e}$ de la lumière excitant notre colorant ;
- faire varier la durée $\tau$ d'exposition.
Nous pouvons en effet écrire le λ de notre loi comme suit : $\lambda = \phi v c i_{e} \tau \, ,$ où $v$ est le volume de solution « vu » par le pixel, $c$ est la concentration de colorant, $\phi$ est le rendement quantique -- un paramètre qui décrit l'efficacité d'excitation du colorant, les pertes de photons de la lumière excitatrice et de la lumière émise, etc.
En pratique, il est plus facile de faire varier le temps d'exposition τ et c'est ce que nous avons fait dans les expériences qui suivent. J'ai ainsi demandé à mes collaborateurs du laboratoire de Peter Kloppenburg à l'Université de Cologne de : choisir 10 temps d'expositions ; pour chacun des 10 temps effectuer 100 expositions successives ; pour chacune des 10 x 100 expositions, enregistrer la valeur $y_{ij}$ prise par la v.a. $Y_{ij}$ au pixel $i,j$ de la CCD.
La raison pour laquelle on introduit une v.a. $Y_{i,j}$ par pixel et qu'il est très difficile d'avoir une illumination ($i_e$), un volume $v$ et un rendement quantique ($\phi$) uniforment. On a donc pour chaque pixel : $$Y_{i,j} \sim G \, p_{i,j} \tau + \sqrt{G^2 \, (p_{i,j} \tau+\sigma_{L}^2)} \, \epsilon_{i,j}\; ,$$ où $p_{i,j} = c \phi_{i,j} v_{i,j} i_{e,i,j}$.
- si notre modèle est correct alors pour chaque pixel $i,j$ nous devons avoir, pour un temps d'exposition donné, une moyenne empirique : $$\bar{y}_{i,j} = \frac{1}{100} \sum_{k=1}^{100} \quad y_{i,j,k} \approx G \, p_{i,j} \tau \; ;$$
- une variance empirique : $$S_{i,j}^2 = \frac{1}{99} \sum_{k=1}^{100} \quad (y_{i,j,k}-\bar{y}_{i,j})^2 \approx G^2 \, (p_{i,j} \tau+\sigma_{L}^2) \; ;$$
- et le graphe de $S_{i,j}^2$ en fonction de $\bar{y}_{i,j}$ devrait être une droite de pente $G$ et d'ordonnée en zéro $G^2 \sigma_{L}^2$.
urllib.urlretrieve("http://xtof.disque.math.cnrs.fr/data/CCD_calibration.npz","CCD_calibration.npz")
CCD_calibration_npz = load("CCD_calibration.npz")
CCD_calibration_piles = CCD_calibration_npz['arr_0']
CCD_calibration_npz.close()
L'objet que nous venons de créer, CCD_calibration_piles, est une liste de 10 piles d'images -- une pile par temps d'exposition -- de dimension 60 x 80 x 100 (le capteur fait 60 x 80 pixels et nous avons 100 expositions par pile). La première pile correspond à un temps de 10 ms et la dernière à un temps de 100 ms.
L'image de la première exposition de 10 ms est ainsi :
imshow(transpose(CCD_calibration_piles[0][:,:,0]),origin='lower')
colorbar()
Les décours temporels de trois pixels voisins sont ici :
subplot(311)
plot(CCD_calibration_piles[0][31,41,:])
ylabel("ADU")
grid()
subplot(312)
plot(CCD_calibration_piles[0][31,40,:])
ylabel("ADU")
grid()
subplot(313)
plot(CCD_calibration_piles[0][31,39,:])
xlabel("Unités de temps (1 unité = 100 ms)")
ylabel("ADU")
grid()
- comme les données vont être analysées comme si les $Y_{i,j,k}$ étaient IID (pour $i$ et $j$ fixés), mais qu'elles ont en fait été enregistrées successivement, il est fortement recommandé de vérifier qu'elles se comportent comme si elles étaient IID ;
- le petit exemple de la figure précédente montre qu'il n'y a pas de dérive ;
- il faut aussi regarder les fonctions de corrélations.
Le graphe de la variance en fonction de la moyenne pour chaque pixel et pour chaque temps d'exposition est obtenu avec :
%config InlineBackend.figure_format = 'png'
ADU_m = [pile.mean(2) for pile in CCD_calibration_piles]
ADU_v = [pile.var(2) for pile in CCD_calibration_piles]
for i in range(10):
scatter(ADU_m[i].flatten(),ADU_v[i].flatten(),0.05)
xlabel(r'$\overline{ADU}$')
ylabel("Var(ADU)")
Ici, je retourne au format png sinon la taille explose !