Exploitation de l’effet de voile atmosphère

Les differents effets du brouillard sur la vision  ont été présentés et modélisés dans le premier chapitre. Trois effets majeurs ont été décrits. Les effets de halo atmosphérique et de voile rétro-diffusé sont plus significatifs la nuit que le jour. L’effet de voile atmosphérique, quant à lui, est spécifique au brouillard diurne. Ce dernier est bien modélisé par la loi de Koschmieder, qui permet d’exprimer la luminance apparente d’un objet en fonction de sa luminance intrinsèque, de sa distance d’observation, de la densité du brouillard et de la luminance du ciel à l’horizon. En exploitant une propriété mathématique de la loi de Koschmieder, il est possible, à l’aide d’une seule caméra, d’instancier ce modèle et d’estimer le coefficient d’extinction du brouillard. La distance de visibilité météorologique est alors directement déduite. Nous proposons, dans ce chapitre, de présenter cette méthode. Voici comment ce dernier est organisé. Tout d’abord, nous présentons notre capteur mono-caméra et l’hypothèse « monde plan » qui lui est étroitement associée, étudions sa portée et sa résolution, et abordons le problème de son calibrage. Puis, à l’aide de l’hypothèse « monde plan », nous étudions d’un point de vue théorique une propriété intéressante de la loi de Koschmieder. Dans un troisième temps, nous montrons comment exploiter en pratique cette dernière. Pour cela, un algorithme d’analyse d’images en trois phases est décrit. Un indice de confiance sur la mesure est construit. Puis une analyse de sensibilité de la méthode conclut le paragraphe. Nous montrons ensuite les points forts et les limites de notre méthode. Pour en combler certaines, nous montrons que nous pouvons nous ramener facilement à l’approche de Pomerleau [Pomerleau, 1997] et rendre les deux approches complémentaires. Ainsi, nous développons une méthode permettant d’estimer la distance de visibilité météorologique dans de nombreuses situations. Finalement, la visibilité météorologique étant correctement estimée, nous présentons une application directe de notre méthode à la restauration du contraste appliquée à des images de brouillard diurne. 

Modélisation du capteur mono-caméra dans son environnement

Dans ce chapitre, la méthode développée repose sur l’utilisation d’une seule caméra embarquée à bord du véhicule. Dans ce paragraphe, nous présentons la caméra utilisée dans les véhicules expérimentaux du livic. Une modélisation, communément appelée hypothèse « monde plan », en vue d’estimer la profondeur des objets dans les images est décrite. Nous étudions ensuite les performances d’un tel capteur, notamment en terme de résolution. Finalement, le problème du calibrage du capteur est rapidement abordé. 

Présentation du capteur utilisé 

Fig. 2.1 – Disposition de la caméra dans le véhicule. Dans ce chapitre, notre capteur est une simple caméra noir et blanc installée derrière le pare-brise du véhicule comme sur la figure 2.1. La figure 2.2 présente la modélisation du capteur dans l’environnement du véhicule. Dans le repère de l’image, la position d’un pixel est donnée par ses coordonnées (u, v). Les coordonnées de la projection du centre optique dans l’image sont désignées par (u0, v0). θ est l’angle entre l’axe optique de la caméra et l’horizontale. vh désigne la position verticale de la ligne d’horizon. Les paramètres intrinsèques de la caméra sont sa longueur focale f, la taille horizontale tpu et verticale tpv d’un pixel. Nous utilisons aussi αu = f tpu et αv = f tpv . De manière usuelle, nous considérons αu ≈ αv = α. 

Modèle de calcul de la profondeur d’un point vu dans l’image 

Etant donné que nous n’utilisons qu’une seule caméra, nous ne pouvons avoir accès à l’infor- ´ mation de profondeur. Nous contournons le problème en faisant l’hypothèse d’une route plane. Ceci nous permet d’associer une distance à chaque ligne de l’image. Dans la suite du paragraphe, nous allons présenter ce modèle de calcul de distance. En utilisant le modèle sténopé pour la caméra, la projection sur le plan image d’un point de coordonnées tridimensionnelles (x,y,z) dans le repère de la caméra s’exprime par :    u = u0 + α x z v = v0 + α y z (2.1) 28 2.1. Modélisation du capteur mono-caméra dans son environnement f d θ X S Y Z C x y z v u v H M plan route plan image f d θ X S Y Z C x y z v u v H M plan route plan image Fig. 2.2 – Modélisation de la caméra dans son environnement. Elle est située à la pseudo-hauteur H dans le repère (S,X,Y,Z) lié à la scène. Ses paramètres intrinsèques sont sa longueur focale f et la taille t d’un pixel. θ est l’angle entre l’axe optique de la caméra et l’horizontale. Dans le repère de l’image, (u, v) désigne la position d’un pixel, (u0, v0) l’intersection entre l’axe optique et le plan image, C le centre optique et vh la position verticale de la ligne d’horizon. D’après la figure 2.2, la ligne horizontale passant par le centre optique fait un angle θ avec l’axe z de la caméra. Ainsi, dans le plan image, la ligne d’horizon s’exprime par : vh = v0 − α tan(θ) (2.2) Grâce à l’équation (2.1), nous déduisons : v − vh α = y z + tan(θ) (2.3) Si l’on se place dans le repère (S,X,Y ,Z) lié à la scène, l’équation (2.3) devient : v − vh α = Y + H Z + tan(θ) (2.4) Un point M situé sur la route à une distance d de l’origine S est paramétré comme suit : M   X −d sin(θ) d cos(θ)  Fig. 2.3 – Illustration de la formule de calcul de profondeur d’un point dans l’image en fonction de l’angle de tangage de la caméra (—) θ =8˚, (– –) θ =5˚. Les paramètres communs utilisés sont H = 2 m, le rapport de la focale de la caméra sur la taille d’un pixel α = 500, et la demi-hauteur de l’image v0 = 144. Horizontalement, on a le numéro de ligne de l’image. Verticalement, on obtient la distance à la caméra en mètres. On en déduit que : v − vh α = H d cos(θ) (2.6) Finalement, la distance d s’exprime donc par : d =    λ (v−vh) si v > vh ∞ si v ≤ vh o`u λ = Hα cos(θ) (2.7) La figure 2.3 représente la distance d calculée par l’équation (2.7) en fonction du numéro de ligne de l’image et ce pour deux angles de tangage différents de la caméra. La formule étant hyperbolique, l’imprécision ou plus précisément la surface couverte par un pixel augmente avec la distance. En conséquence, l’estimation de la distance sera moins précise et stable pour de grandes distances que pour des petites distances. Ceci fait l’objet du paragraphe suivant. 2.1. Modélisation du capteur mono-caméra dans son environnement Fig. 2.4 – Représentation de la distance ∆ couverte par un pixel à la distance d en fonction de l’angle de tangage de la caméra(—) θ = 5˚, (– –) θ = 8˚. Les paramètres communs utilisés sont la hauteur de la caméra H = 2m, le rapport de la focale de la caméra sur la taille d’un pixel α = 500 et la demi-hauteur de l’image v0 = 144.

Précision du calcul de distance en fonction du tangage de la caméra 

Grâce aux équations (2.2) et (2.7), nous pouvons calculer la surface couverte par un pixel à la distance d en fonction de l’angle de tangage de la caméra : ∆(θ, d) = λ (bvh + λ d c − vh) − λ (dvh + λ d e − vh) (2.8) Dans cette formule, bxc désigne la partie entière de x et dxe, l’entier supérieur ou égal à x. De manière usuelle, l’angle de tangage de la caméra est de 8˚, ce qui permet un bon compromis de contraste entre la route et le ciel. En augmentant cette valeur, on fait remonter la ligne d’horizon dans l’image, ce qui nous permet de faire diminuer la surface couverte par un pixel et de gagner en précision pour les grandes distances, ce qui est intéressant dans le cas de notre application. A la vue de la figure 2.4, nous pouvons considérer que l’angle de tangage n’a d’influence que pour des mesures de distance supérieures à 250 m. 

Calibrage du capteur 

Les paramètres intrinsèques de la caméra ayant été établis une fois pour toutes par des méthodes classiques, nous n’aborderons pas ce problème ici. Le principal problème restant à résoudre est l’obtention des paramètres extrinsèques du capteur. Nous proposons une solution simple qui nécessite uniquement une image et le calcul d’un paramètre. En utilisant le modèle mono-caméra présenté précédemment, une simple estimation du paramètre λ permet de calibrer le capteur. Pour effectuer cela, la connaissance de la distance réelle entre deux points et de leurs coordonnées v1 et v2 dans l’image est suffisant (cf. figure 2.5). En effet, en utilisant la formule 2.7, nous établissons l’équation suivante : λ = d1 − d2 ³ 1 v1−vh − 1 v2−vh ´ (2.9)  d ,v d ,v 10m v Fig. 2.5 – Exemple d’image utilisée pour le calibrage des paramètres extrinsèques de la caméra et capturée sur le site de calibrage des pistes de Versailles Satory. v1 et v2 représentent les lignes correspondant à deux objets situés respectivement aux distances d1 et d2. vh correspond à la ligne d’horizon. 2.2 Estimation de la distance de visibilité météorologique Dans ce paragraphe, nous présentons le cœur de notre méthode permettant d’estimer la distance de visibilité météorologique. Dans un premier temps, nous démontrons que, sous certaines hypothèses, la loi de Koschmieder a une propriété mathématique intéressante, ce qui en théorie nous permet d’atteindre notre objectif. Puis, nous présentons comment nous exploitons cette propriété de manière pratique à l’aide, entre autres, d’un algorithme de croissance de région. Nous finissons par définir un indice de confiance sur la mesure effectuée. Enfin, nous concluons le paragraphe par une courte analyse de sensibilité de la méthode présentée.

Mise en évidence d’un point d’inflexion

Dans le paragraphe 1.3.1, nous avons présenté précisément la loi de Koschmieder qui permet de modéliser la luminance apparente d’un objet en fonction de sa luminance intrinsèque, de sa distance d’observation, de la densité du brouillard et de la luminance du ciel à l’horizon. Dans ce 32 2.2. Estimation de la distance de visibilité météorologique qui suit, nous allons étudier les propriétés mathématiques de cette formule et déduire l’existence d’un point d’inflexion détectable sur l’image, sur lequel se fonde notre solution. Après un changement de variable de d en v à partir de l’équation (2.7), l’équation (1.5) devient : L = L0 − (L0 − Lf )(1 − e −k λ v−vh ) (2.10) En dérivant l’équation (2.10) par rapport à v, on obtient : dL dv = kλ(L0 − Lf ) (v − vh) 2 e −k λ v−vh (2.11) Les courbes représentatives de L et de sa dérivée sont tracées sur la figure 2.6 pour différentes valeurs du coefficient d’extinction k. De manière qualitative, plus le brouillard est dense, plus l’objet se confond rapidement avec la luminance du ciel et plus le maximum de la dérivée se réduit et s’éloigne de la ligne d’horizon. En dérivant à nouveau L par rapport à v, on obtient : d 2L dv2 = kA(v)e −k λ v−vh µ kλ v − vh − 2 ¶ (2.12) o`u A(v) = λ(L0 − Lf ) (v − vh) 3 . L’équation d 2L dv2 = 0 a deux solutions. La solution k = 0 n’a aucun intérêt. Ainsi, la seule solution utile est l’équation (2.13) : k = 2(vi − vh) λ = 2 di (2.13) o`u vi représente la position du point d’inflexion et di sa distance à la caméra. Ainsi, on obtient le paramètre k de la loi de Koschmieder dès lors que l’on connaˆıt vi . De plus, l’équation (2.13) possède la propriété remarquable limvi→vh k = 0. Celle-ci peut servir à détecter la présence de brouillard. En effet, si vi est supérieur à vh, alors on détecte du brouillard, sinon, il n’y a pas de brouillard. Grâce aux équations (1.13) et (2.13), nous déduisons la distance de visibilité météorologique Vmet : Vmet = 3λ 2(vi − vh) (2.14) Si vv désigne la ligne de l’image représentative de la distance de visibilité, nous avons : vv = 2vi + vh 3 (2.15) . Finalement, des valeurs de vi et vh, nous déduisons les autres paramètres de la loi de Koschmieder en utilisant Li et dL dv |v=vi , respectivement les valeurs de la fonction L et de sa dérivée en v = vi  Fig. 2.6 – Différentes courbes représentatives de (a) la loi de Koschmieder et (b) de sa dérivée(niveau de gris en fonction numéro de ligne de l’image). Les paramètres utilisés sont λ = 1000, vh = 100, L0 = 50, Lf = 220, k = 0, 01 (—), k = 0, 05 (−·−), k = 0.09 (. . .), k = 0.13 (– –). Le point d’inflexion, sur la dérivée, pour une densité de brouillard donnée, est clairement mis en évidence. 34 2.2. Estimation de la distance de visibilité météorologique 

Mise en œuvre pratique de la méthodologie 

Comme nous venons de le montrer dans le paragraphe 2.2.1, pour estimer la distance de visibilité météorologique, nous avons besoin de mesurer la position verticale de deux objets dans l’image de la route : le point d’inflexion et la ligne d’horizon. Dans ce paragraphe, nous abordons ces deux problèmes de manière pratique par analyse d’images. Ainsi, deux algorithmes indépendants sont présentés pour les résoudre. L’estimation du point d’inflexion est très détaillée. L’estimation de la ligne d’horizon est rapidement abordée, car de nombreuses solutions existent dans la littérature. Tout d’abord, pour estimer la position du point d’inflexion, le problème est de savoir sur quel objet mesurer la variation de luminance. Dans notre contexte, l’objet le plus adéquat est la route. En effet, c’est un objet sombre, toujours présent dans la scène et c’est un lieu de contact avec le ciel. C’est également un objet suffisamment étendu pour y percevoir une variation spatiale de luminance. Pour être cohérent avec le modèle de Koschmieder qui suppose une luminance intrinsèque L0, nous supposons que la route est homogène et que sa luminance n’est affectée que par le phénomène de voile atmosphérique. Par conséquent, l’algorithme présenté dans les paragraphes suivants, recherche dans l’image une surface présentant une variation continue et faible du niveau de gris lorsqu’on se déplace de ligne en ligne. Puisque la route finit par se confondre avec le brouillard, cette surface inclut également le ciel, de luminance Lf à l’infini. Discontinuités de luminance Dans un premier temps, nous extrayons les contours de l’image de façon à mettre en évidence les ruptures importantes de contraste que constituent les bords de voies, les véhicules suivis ou croisés, les arbres… Ceci est réalisé par un filtre de Canny-Deriche [Deriche, 1987]. Nous notons E l’ensemble de ces contours. Pour sélectionner les contours, nous utilisons un seuillage par hystérésis (cf. figure. 2.7). Nous notons tH and tL les seuils haut et bas de celui-ci. tL et tH sont fixés à des valeurs relativement élevées pour éviter la prise en compte de bruit dans la détection des contours et pour éviter d’obtenir une rupture au niveau de la ligne d’horizon. En fait, tL et tH ne peuvent pas être constants car leur valeur dépend directement du niveau de visibilité dans l’image. Par chance, la détection des contours est faiblement sensible aux valeurs de seuil comme le montre la figure 2.8. Fig. 2.7 – Principe du seuillage par Hystérésis : les points (i, j) sont marqués par deux seuils. Le seuil bas tL génère les points ”candidats” mais est très sensible au bruit. Le seuil haut tH est une bonne indication de l’existence d’un contour. Fig. 2.8 – Sensibilité du filtre de Canny-Deriche aux valeurs de seuils. (a) Image originale (b) tL = 12 et tH = 30 (c) tL = 12 et tH = 25 (d) tL = 12 et tH = 20 (e) tL = 9 et tH = 20 (f) tL = 6 et tH = 20. Principe de la croissance de région Par la suite, nous effectuons une croissance de région. Dans ce paragraphe, nous définissons son objectif, ses germes ainsi que ses paramètres fondamentaux. L’objectif de l’algorithme est de trouver une région dans l’image présentant une variation minime de niveau de gris de ligne en ligne quand on la traverse de bas en haut, de façon à être compatible avec la loi de Koschmieder. Dans cette perspective, on choisit les germes de la croissance de région comme étant les pixels d’une ligne du bas de l’image dont le niveau de gris est proche du niveau de gris médian de cette ligne. En effet, compte tenu de la position et des caractéristiques optiques de la caméra, la majorité des pixels de cette ligne représente le revêtement routier. Ainsi, seuls les pixels de la route sont pris en compte, comme sur la figure 2.12a, évitant par exemple de faire croˆıtre certains germes sur un marquage routier. De même, seuls les trois pixels au-dessus du pixel courant (figure 2.9) peuvent être agrégés à la région R. Cette technique nous permet de contourner les objets ne faisant pas partie de la région d’intérêt. Fig. 2.9 – Principe de la croissance de région. Seuls les pixels gris clair peuvent être agrégés au pixel gris foncé. Les autres pixels du voisinage sont ignorés. Il reste à définir le seuil de gradient à considérer pour passer d’une ligne à l’autre de l’image. Du paragraphe 2.2, on peut déduire le gradient vertical maximal Gmax qui existe entre deux lignes successives de l’image. D’un point de vue théorique, s’il n’y a pas de diffusion atmosphérique, Gmax est égal à Lf − L0 au niveau du point d’inflexion. Malheureusement, cette valeur est trop peu contraignante pour la croissance de région. Ainsi, nous avons choisi de limiter la valeur de 36 2.2. Estimation de la distance de visibilité météorologique Gmax à tL. En choisissant un tel seuil, on limite le saut admissible de gradient entre deux lignes successives à une valeur inférieure à celle définie par la loi de Koschmieder. La croissance de région est donc compatible avec la loi de Koschmieder, mais en procédant ainsi, les brouillards qui correspondent à une distance de visibilité supérieure à 400 m ne sont pas détectables. Ce n’est, en fait, pas un problème pour les applications visées pour deux raisons : – A de telles portées l’estimation de la distance est de toute manière très peu précise (voir paragraphe 2.1.3). – La densité du brouillard est tellement faible qu’il n’est pas perturbant pour le conducteur et les instruments de mesure. Toujours afin de privilégier le mouvement vertical lors de la croissance, nous avons choisi de prendre des seuils plus contraignants dans les directions obliques que dans la direction verticale. Nous introduisons donc les notations Gi max o`u i ∈ {−1, 0, 1} (cf. figure 2.11), avec les contraintes suivantes : G −1 max = G 1 max < G0 max ≤ tL (2.18) Des valeurs usuelles de ces différents seuils sont de l’ordre de 4 à 8. Conditions d’agrégation d’un pixel à la région cible Dans ce paragraphe, nous allons définir les quatre conditions d’agrégation d’un pixel P(i, j) à la région d’intérêt R. Elles sont présentées dans leur ordre d’exécution qui est aussi celui d’un temps de calcul croissant. Ceci permet d’éviter de tester les conditions les plus gourmandes en temps de calcul si les premières sont rejetées. Le pixel P(i, j) est agrégé à la région R en construction si quatre conditions sont remplies : – Première condition, le pixel n’appartient pas à la région R. P(i, j) 6∈ R (2.19) – Deuxième condition, le pixel ne fait pas partie d’un contour détecté par le filtre de CannyDeriche. P(i, j) 6∈ E (2.20) – Troisième condition, le pixel a une certaine similarité avec le germe Pg. La similarité est évaluée en calculant la différence (2.21) de niveau de gris entre le pixel étudié et le pixel Pg. P(i, j) − Pg ≤ ρnr min k∈{−1,0,1} G k max (2.21) o`u ρ < 1 et nr désigne le nombre de lignes entre P(i, j) et Pg. Cette condition constitue un garde-fou sur la croissance de région évitant d’agréger des pixels trop différents du germe pour la ligne considérée. Sans cette condition, si l’on considère une tolérance de gradient vertical égale à huit, des pixels noir et blanc pourraient être agrégés en moins de 32 itérations. – Quatrième condition, le pixel est similaire au pixel situé en dessous. Dans ce cas, la similarité est évaluée en utilisant un filtre inspiré de Nagao [Nagao et Matsuyama, 1979]. Ce filtre calcule la moyenne et la dispersion pour neuf masques différents. Il a pour résultat la moyenne du masque possédant la plus petite dispersion. Dans notre cas, le résultat est le médian (2.24) du masque possédant l’étendue en niveaux de gris (2.22) la plus petite (2.23), de façon à être plus robuste au bruit et plus rapide en temps de calcul. Nous avons également choisi de régulariser la forme des masques en employant des masques carrés [Demigny et al., 1993]. Cela nous permet de réutiliser les calculs précédents, quand la croissance passe à la ligne suivante, ce qui n’était pas possible avec les masques originaux de Nagao.