29
Photogrammétrie numérique à l'attention des DC2 Raphaële Héno – Dias v1 octobre 2008 Photogrammétrie numérique – DC2 – 2008 1/29

à l'attention des DC2 Raphaële Héno – Dias v1 octobre ...darkyoup.free.fr/tipe/Photogram...pdf · Caméra numérique professionnelle, ... Satellite à 800km de la terre,

  • Upload
    lyminh

  • View
    219

  • Download
    0

Embed Size (px)

Citation preview

Photogrammétrie numérique

à l'attention des DC2

Raphaële Héno – Dias

v1 octobre 2008

Photogrammétrie numérique – DC2 – 2008 1/29

I. INTRODUCTION

Définition : technique permettant de mesurer en 3D sur des images (par essence en 2D) . On simule ce qui se passe en vision naturelle où notre cerveau fusionne les perspectives rendues par chacun des deux yeux ; c'est à dire que l'on prend deux images de la scène d'intérêt, sous deux angles différents.

Rappels : Photogrammétrie analogique (la photogrammétrie historique). Elle utilise des films argentiques et du matériel composé de mécanique et d'optique. Exemples de tels appareils : Planimat (Zeiss), B8 (Wild)

Principe : [Orientation interne] on reconstruit mécaniquement les deux faisceaux perspectifs (centrage des films et introduction de la distance principale), et [Orientation relative] on leur donne une position relative similaire à la position du moment de la prise de vues, en jouant sur les paramètres de translation et de rotation de chacun des deux faisceaux, enfin [Orientation absolue] on met le modèle à l'échelle souhaitée, et on fait correspondre ses horizontales à celles du terrain.

La phase d'orientation relative, la plus délicate, consiste à annuler la parallaxe transversale en 6 points du modèle [les 6 points de Von Grüber].

Photogrammétrie numériqueElle utilise des images numériques et du matériel informatique.

Origine des images numériques :Films argentiques scannésImages de caméras numériques professionnellesImages d'appareil photo numérique amateurImages satelliteImages radarImages médicales, échographies, etc.

Caractéristiques des images numériques :Nombre de pixels en ligne et en colonneTaille du pixel au sol [résolution spatiale]Colorimétrie [résolution spectrale]Temps de revisite, pour les images satellite [résolution temporelle]Emprise d'une image au sol [fauchée]

Matériel :Caméra numérique professionnelle, pour faire des images aériennes jusqu'à qq cm de résolution solAppareil photo numérique amateur, pour faire des images «terrestres », ou aériennesSatellite à 800km de la terre, pour des résolutions sol jusqu'à 50cm

Les images aériennes sont prises d'un ballon, d'un ULM, d'un cerf-volant, d'un hélicoptère, d'un avion, d'un drone [liste des « vecteurs »]

Photogrammétrie numérique – DC2 – 2008 2/29

Exemples d'images satelliteRésolutionspatiale

fauchée Résolution spectrale

Résolutiontemporelle

Stéréo

Météosat 2,5 km 12500 km X Ts les jours non

Landsat 15 à 60m 185 km X non

Spot 5 5m 60 x 120 km panchro oui

Ikonos 70cm 11 km panchro 3 jours oui

Quickbird 70cm 16 km panchro 3 jours oui

Exemples de caméras professionnelles :les argentiques : RMKTOP (Zeiss) ; RC10, 20, 30 (Leica), de grand format (24x24cm), à focale fixe (88mm, 153mm, 200mm, ou 300mm)/les numériques : DMC (Zeiss-Intergraph), ADS40 (Leica), Ultracam (Vexcel), de grand format (combinaision de plusieurs matrices CCD) et focale fixe.Les numériques amateurs : longue liste : nikon, canon, etc. Zoom, (focale variable), distorsion

Le matériel de la photogrammétrie numérique :un ordinateur (PC, mac, ...), un écran cathodique, éventuellement un écran plat, une souris de bureau, une souris 3D, un stylo, une tablette graphique, une dalle stéréo, des lunettes stéréo.

http://www.stereo3d.com/3dhome.htm

Visu stéréo par analgyphes

L'image gauche est colorée en bleu, l'image droite est colorée en rouge. Des lunettes aux verres bleu et rouge permettent de sélectionner ce que voit chaque oeil.

Photogrammétrie numérique – DC2 – 2008 3/29

Visu stéréo via des lunettes clignotantes

Des lunettes clignotantes sont branchées sur l'ordinateur, et clignottent à la même fréquence F que l'écran (nécessairement supérieur à 100Hz). A l'instant T0, l'image gauche apparaît sur l'écran, et le verre droit s'obscurcit. A T0+1/F, l'image droite apparaît sur l'écran, et le verre gauche s'obscurcit, et ainsi de suite de telle sorte que le cerveau a l'impression de voir en permanence l'image gauche avec l'oeil gauche et l'image droite avec l'oeil droit.

Visu stéréo via des lunettes polarisées

On utilise des lunettes polarisées, un verre filtrant les ondes selon un angle θ, l'autre filtrant les ondes inclinées selon un angle θ+∏/2. On remarque que les deux verres posés l'un sur l'autre bloquent la lumière. On superpose à l'écran une dalle branchée à l'ordinateur. La dalle change de polarisation (alternativement θ et θ+∏/2) à une fréquence F. A l'instant T0, l'image gauche apparaît sur l'écran, et la dalle prend la polarisation θ. A T0+1/F, l'image droite apparaît sur l'écran, et la dalle prend la polarisation θ+∏/2, et ainsi de suite de telle sorte que le cerveau a l'impression de voir en permanence l'image gauche avec l'oeil gauche et l'image droite avec l'oeil droit.

Photogrammétrie numérique – DC2 – 2008 4/29

II. CALCULS MATRICIELS DE BASE

voir cours de Jean Debord : (http://www.unilim.fr/pages_perso/jean.debord/math/matrices/matrices.htm)

Puis, résolution des systèmes d'équations linéaires surabondants (plus d'équations que d'inconnues à déterminer) par la méthode des moindres carrés.

Le système doit être posé sous la forme

A x = b

où A est une matrice de dimensions n,mx la matrice des inconnues, de dimensions m,let b les mesures, de dimensions n,l

La méthode des moindres carrés donne une solution qui colle au mieux à toutes les mesures. Le résidu de la solution trouvée par rapport aux mesures vaut : Ax - b

Photogrammétrie numérique – DC2 – 2008 5/29

III. ÉQUATION DE COLINÉARITÉ

Formule d’image correspondant à la perspective centrale

On cherche la fonction qui donne l’image m d’un point M, c'est à dire la fonction f telle que m=f(M)

Comment se fait un image, une photo ? Une source lumineuse, le soleil par exemple, illumine le terrain. Comsidérons un point M du terrain : les rayons lumineux sont réfléchis dans toutes les directions, et un des rayons réfléchis rentre dans la caméra en F, et imprime son image en m.

m est l’intersection de la droite SM et du plan P.

SmSM λ= (mise en équation de la colinéarité ; λ est le facteur d’échelle).

0. =CmK ( K est la normale au plan P)

Photogrammétrie numérique – DC2 – 2008 6/29

M

SZ

Y

F

m C

X

Expression en coordonnées :

1) SmSM λ= on écrit a priori (M-S) = λ (m-S), mais les coordonnées de m ne sont pas connues dans le système terrain. Il faut écrire :R(M-S) = λ (m-F), où R est la matrice de passage du système terrain au système fond de chambre., et F est le point de coordonnées fond de chambre pyx ppappa ,, .

2) 0. =CmK 0=Kmt

1) R(M-S) = λ (m-F)est équivalent à :Kt R(M-S) = λ Kt (m-F)

Grâce à 2), on a KF

SMKRt

t )( −−=λ

D’où )()(

SMRSMKR

KFFm t

t

−−

−= (équation de la perspective)

C'est l'équation de colinéarité, ou la formule d'image pour les images modélisables par une perspective, ou la fonction de localisation de ces images.

F K R S M3 paramètres constante 3 paramètres 3 paramètres c’est la

variableL’équation compte 3+3+3 =9 paramètres

Le géoréférencement des images consiste à déterminer pour chaque image sa formule d'image, c'est à dire la position S et l'orientation R. Pour les déterminer (3 paramètres pour S, 3 paramètres pour R, en supposant que F est donné par le certificat d'étalonnage), il faut mesurer au moins 3 points d'appui (puisqu'1 point d'appui mesuré donne 2 équations).Il existe des moyens externes pour connaître R et S (centrale à inertie et GPS). Si l'avion est équipé de tels systèmes, on peut obtenir des images directement géoréférencées.

Développement de l'équation de colinéarité : l'équation des 11 paramètres

cx

F cyp

Photogrammétrie numérique – DC2 – 2008 7/29

coordonnées du PPA

distance principale

F est souvent donné par le certificat d'étalonnage de la caméra.

XS

S YS

ZS

XM

M YM

ZM

x= )()(

333231333231

131211131211

SSSMMM

SSSMMMc ZrYrXrZrYrXr

ZrYrXrZrYrXrPx++−++++−++

réduction au même dénominateur :

x=

)()()()()()()(

333231333231

133312321131133312321131

SSSMMM

ScScScMcMcMc

ZrYrXrZrYrXrYprxrYprxrXprxrZprxrYprxrXprxr

++−++−−−−−−−+−+−

y= )()()()()()()(

333231333231

133312321131233322322131

SSSMMM

ScScScMcMcMc

ZrYrXrZrYrXrYpryrYpryrXpryrZpryrYpryrXpryr

++−++−−−−−−−+−+−

Autre formulation (changement de variable):

x= 4321

4321

cZcYcXcaZaYaXa

MMM

MMM

++++++

y=4321

4321

cZcYcXcbZbYbXb

MMM

MMM

++++++

Ici, il faut mesurer 6 points pour déterminer les 11 coefficients de la formule d'image de chaque image du couple. Ces mesures donneront un système d'équations linéaires que l'on résoudra par la méthode des moindres carrés.

Cette formule paraît plus simple, mais elle est plus dangereuse car les paramètres supplémentaires (11 au lieu de 9, voire souvent 6 si la caméra est connue) peuvent absorber des fautes de mesure.

En conclusion, la formule des 11 paramètres ne sert qu'à trouver une solution approchée de R et S.

Précision des mesures sur un couple stéréoscopique :

SigmaXY = ½ pixel-sol

SigmaZ = ½ pixel-sol x H/B

(en effet, la précision d'appréciation de la profondeur dépend de la configuration stéréoscopique : plus le rapport Base sur Hauteur de vol est fort, plus les rayons perspectifs homologues s'intersecteront franchement, plus la précision en Z sera forte).

Photogrammétrie numérique – DC2 – 2008 8/29

formule dite des 11 paramètres –si on divise haut et bas par c4- (et non pas 12) car le système est invariant par changement d’échelle)

Application NumériqueSoit un couple d'images numériques à 60cm, prises avec une caméra IGN 4000 x 4000 pixels (pixel capteur : 9 microns ; focale 60mm). Le recouvrement stéréoscopique est de 60%. Quelles sont les précisions planimétrique et altimétrique attendues ?

SigmaXY = ½ pixel-sol = 30cmSigmaZ = ½ pixel-sol x H/B

avec H = focale x pixel-sol / pixel-capteur = 4000 met B = (1-R) x 4000 x pixel-sol = 960

d'où sigmaZ=1.25m

Photogrammétrie numérique – DC2 – 2008 9/29

IV. CORRECTIONS À APPORTER AU MODÈLE

Pour atteindre la plus grande précision, la formule énoncée précédemment ne suffit pas. Il convient de prendre en compte certains paramètres physiques. On écrira :

m = perspective de M + corrections

Les différentes corrections à apporter portant sur :

1. La courbure de terre2. La réfraction atmosphérique3. La distorsion Optique4. Les déformations du film

Correction n°1 : Courbure de terreL’espace cartographique dans lequel on exprime les coordonénes terrain, n’est pas un référentiel euclidien. En effet, la terre est ronde, l’espace Z = constante n’est donc pas un plan ! On doit donc trouver la transformation du système cartographique vers un système euclidien.

Passer de (X,Y)euclidien à (X,Y)carto est le résultat d’une projection cartographique. Pour le Z, il faudrait procéder à des mesures de nivellement pour tenir compte des champs de potentiels.

La mise en œuvre de cette méthode demanderait l’implémentation d’un géoïde fin du globe et une base de données de toutes les projections possibles. Aussi, plutôt que de les recenser toutes, on estimera qu’elles sont toutes localement équivalentes à une similitude près; on peut donc prendre n’importe laquelle. On décide d’utiliser la projection « stéréographique oblique » (c’est une projection conforme) de la sphère de courbure moyenne (les chantiers de photogrammétrie couvrant généralement des zones relativement petites, on peut assimiler le géoïde à une sphère.) sur le plan tangent au chantier.

Photogrammétrie numérique – DC2 – 2008 10/29

carte

A: centre du chantier

O

B

I

h

I ’

hJ ’

J

x

carte

A: centre du chantier

O

B

I

h

I ’

hJ ’

J

x

I est le point à transformer ; J est le résultat. Pour calculer les corrections à apporter aux coordonnées de I pour passer dans un repère tridimensionnel, on utilise le fait que I’ et J’ sont inverses dans une inversion de pôle C’ et de puissance 4R2

On démontre en utilisant les relations trigonométriques, que

J-I =

xhR

xR

xR

3

2

24

2

D’où la formule :

−+−−

−−−−

+

=

])()[(

))((2))((2

20

20

00

00

YYXXCZZYYC

ZZXXC

ZYX

ZYX

cartocarto

cartocarto

cartocarto

cartoeuclidien

avec C=1/2R (R = 6378 km)

Attention: Dans des conditions courantes, le terme de correction altimétrique est de loin le plus significatif.

Attention cependant, pour des chantiers de grande envergure, les corrections planimétriques ne peuvent plus être ignorées.

Application numériqueSoit une zone de 25 km de long sur 20 km de large, avec des altitudes comprises entre 200 et 1000m.

1. Donner des coordonnées 3D possibles pour les points A, B, C, D et O, A, B, C et D étant les 4 coins de la zone, O en étant le barycentre.

2. Calculer, en appliquant les formules vues en cours, les erreurs Dz et Dxy que l’on commet si l’on néglige la courbure de terre.

Photogrammétrie numérique – DC2 – 2008 11/29

Xcarto Ycarto Zcarto DeltaX DeltaY DeltaZA 0 0 1000 -1,57 -1,25 -20,09B 25000 0 623 0,12 -0,10 -20,09C 25000 20000 305 0,02 0,02 -20,09D 0 20000 507 -1,04 0,83 -20,09Origine 12500 10000 200 0,00 0 0

Correction n°2 : le rayon MF n’est pas une droite

Avant d’atteindre la caméra, les rayons lumineux traversent l’atmosphère dont la densité diminue avec l’altitude. Cette variation de la densité (indice de réfraction), provoque une courbure des rayons.

L’indice de réfraction varie dans le visible en fonction de la température et de la pression seulement : n=f(P,T°)=f(altitude). On sait qu’à une altitude supérieure à quelques centaines de mètres (altitude facilement atteinte pour la photogrammétrie), les variations de l’indice sont faibles.

n se modélise bien par la formule suivante : 1+ Ze610.105.000278.0

−− (formule donné par l'OACI, l'organisation de l'aviation civile internationale)

Etant donnée par ailleurs la loi de Descartes pour la réfraction (n1 sin i1 = n2 sin i2), on peut calculer la vraie trajectoire du rayon lumineux, et calculer un DeltaZ menant au point M'(schéma ci-dessus) duquel le rayon lumineux paraît être issu.

Correction n°2 ' : prise en compte de la pressurisation de l’avion

Lorsque l’avion est pressurisé, le rayon lumineux est réfracté une seconde fois à l’entrée de la cabine, où la pression est plus élevée qu’à l’altitude de vol. On peut modéliser cette réfraction de manière analogue à la précédente.

La plupart des ouvrages se trompent donc en affirmant qu’il ne faut pas corriger la réfraction atmosphérique si l’avion est pressurisé !

Correction n°3 : prise en compte des défauts de l’optique

Si l’optique des caméras utilisées était une Optique de Gauss (c’est à dire qui vérifie les trois conditions suivantes : lentille mince, rayon proche de l’axe, ouverture faible), alors on aurait bien, comme le reproduit la perspective mathématique : tout rayon passant par le centre non dévié.

Photogrammétrie numérique – DC2 – 2008 12/29

M

F

m’: point théoriquepoint sur la photo : m

M’

y

En réalité, comme il est impossible de placer le diaphragme (qui sert à régler l'ouverture de la lentille) au centre de l'ensemble des lentilles qui composent l'objectif, il y a de la distorsion soit en barillet, soit en coussinet.

Distorsion en barillet Distorsion en coussinet

La distorsion se modélise par une correction radiale et symétrique autour du Point Principal de Symétrie, le PPS (c'est intersection du fond de chambre et de l'axe optique). La distorsion dr se modélise souvent sous la forme : dr = 753 crbrar ++

avec r, distance radiale au PPS :r= ( ) ( ) 22

ppspps yyxx −+−

dx=drr

x

dy=drr

y

dr est exprimé en microns sur les caméras argentiques, en pixels sur les caméras numériques. L'ordre de grandeur de la distorsion sur les caméras argentiques de photogrammétrie est de qq microns, et de plusieurs dizaines de pixels sur les caméras numériques grand public.

Les objectifs à longue focale donnent plutôt une courbe de distorsion positive, tandis que les grands champs (courte focale), donne une courbe de distorsion négative.

Le Point Principal d’Auto-collimation (PPA), projection orthogonale du point d’entrée de l’appareil sur le fond chambre. est nécessaire pour déterminer la formule d’image, qui est la relation entre les coordonnées « caméra » et les coordonnées « terrain » en trois dimensions.

ÉTALONNAGE GÉOMÉTRIQUE D'UNE CAMÉRA

C'est la détermination de ses paramètres intrinsèques, c'est à dire F (xppa, yppa, f) et xpps, ypps, et les coefficients a, b, c du polynôme de distorsion.

Pour ce faire:

Photogrammétrie numérique – DC2 – 2008 13/29

● on photographie un polygone d'étalonnage (site équipé de nombreuses ciblettes dont les coordonnées sont précisément déterminées).

● On mesure les ciblettes sur toutes les images● un logiciel d'étalonnage pose les équations de colinéarité et résout, pour chaque image, les

inconnues habituelles, mais aussi F (xppa, yppa, f), xpps, ypps, et les coefficients a, b, c du polynôme de distorsion.

Correction n°4 : prise en compte de la déformation du film

L’équation de colinéarité a été écrite en supposant que l’intersection du rayon lumineux se fait avec un plan (1), et que ce « plan » (on verra s’il s’agit réellement d’un plan) est stable dans le temps (2).

(1) Planéité de l’émulsion photographique ?→ oui, si l’on utilise du matériel de photogrammétrie car les chambres métriques (à distinguer des appareils photos amateurs…)sont équipées de pompe à vide pour assurer la planéité du film. Ce dispositif est efficace. Seuls des conditions extrêmes (fort vent de sable par exemple) peuvent le mettre à mal si une poussière vient se glisser entre le film et le fond de chambre.

Notons que si l’on utilise un appareil photo quelconque, sans dispositif d’aspiration, on aura des problèmes de déformation, d’autant plus importants que le format est grand.

(2) Stabilité temporelle de l’émulsion ? → non

schéma du film :

La gélatine en absorbant l’humidité, gonfle, et entraîne la déformation des deux couches de gélatine. Le film, de par sa fabrication, est dimensionnellement instable.

Notons qu'on pourrait limiter la déformation en climatisant l’appareil et la pièce où il se trouve.

Comment faire pour tenir compte de cette déformation, c’est à dire pour connaître l’état du film au moment de la prise de vues ?

Photogrammétrie numérique – DC2 – 2008 14/29

La déformation atteint 410 − en facteur d’échelle (c’est à dire 25 microns sur 25 cm de film !)

la caméra est équipée de 8 petits collimateurs qui vont projeter une image sur le film. On les appelle les repères de fond de chambre. Leurs coordonnées sont connues au moment où l’on prend la photo grâce au certificat d’étalonnage de la caméra (qui devra être réactualisé en cas de choc).

gélatine antihalo

polyester

gélatine + sels d’argent

temps

Orientation interne : on va mesurer les repères de fond de chambre au moment de l’exploitation des clichés, et calculer la meilleure affinité (transformation choisie pour modéliser la déformation du film et le changement de repère entre le système « fond de chambre » et le système «image»).

ligne= cybxa +×+×colonne= fyexd +×+×

On a 6 paramètres à trouver par la méthode des moindres carrés à partir de n équations (n = 2xnombre de repères). Le résidu, c’est à dire l’écart entre le modèle et les mesures réelles, est de l’ordre de 0.3 à 0.5 pixel (une valeur supérieure laissera supposer que le certificat d’étalonnage utilisé n’est pas à jour). Pour mieux tenir compte des déformations du film, il faudrait avoir des informations non seulement sur le bord du cliché (les repères de fond de chambre), mais aussi à l’intérieur. On évite bien sûr de placer des repères au milieu du cliché tout simplement pour préserver la zone photographiée !

Récapitulatif : Mise en place puis utilisation des images géoréférencéesIl y a donc deux phases à distinguer : la mise en place des images, qui permet de les géoréférencer, et leur utilisation.

Mise en place des images (cas des images modélisables par une perspective centrale)

1. Orientation interne pour lier les coordonnées « image » aux coordonnées « fond de chambre »

mesure des repères de fond de chambre → détermination d’une affinité (une par image)

2. Orientation « externe » pour lier les coordonnées « terrain » aux coordonnées « fond de chambre »

mesure des points d’appui → détermination des formules d’image (une par image)

Utilisation des images mises en place (valable pour des images de capteur quelconque dès que la formule d’image est connue).

exemple n°1 : stéréorestitution (cf. plus loin boucle temps réel)exemple n°2 : fabrication d’orthoimagesexemple n°3 : calcul de MNT par corrélation automatiqueexemple n°4 : étalonnage d’une caméra

Photogrammétrie numérique – DC2 – 2008 15/29

Exemple d’application de la formule d’image : la saisie 3D sur un appareil de restitution

Photogrammétrie numérique – DC2 – 2008 16/29

2f

m1(x1,y) m2(x2,y2)

ballonnet au sol ?

ouinon

clic suite de la saisie

M(X,Y,Z)M(

M(X,Y,Z)tri

application des paramètres de la projection stéréographique oblique

correction de réfraction atmosphérique (∆Z)

OI1 OI2

- ajout de la distorsion optique-changement de repère fond de chambre → image

m1(x1,y1)image-gauche m2(x2,y2) image droite

1f

périphérique d’entrée 3D

ordinateur

oculaires

mouvement des porte-clichés

ordinateur

)()(

)(1 1111

111 SMR

SMKRKFFMfm t

t

−−

−==

)()(

)(2 2222

222 SMR

SMKRKFFMfm t

t

−−

−==

V. OBTENTION DE LA FORMULE D'IMAGE

Equation de colinéarité )()(

SMRSMKR

KFFm t

t

−−

−=

Les inconnues sont ici R (la matrice de passage du repère terrain au repère fond de chambre) et S (le centre de la perspective dans le système terrain). On suppose F connu grâce au certificat d’étalonnage de la caméra.

Il y a donc 6 inconnues, déterminables par 3 points dont les coordonnées terrain sont connus, et que l’on peut mesurer sur les photos (points d‘appui).

Si F n’est pas donnée dans le certificat d’étalonnage, on a alors 9 inconnues, déterminables par 5 points dont les coordonnées terrain sont connus, et que l’on peut mesurer sur les photos.

Procédure opérateur (sur un couple stéréo)

Mesure de points d'appui (Ground Control Point)Mesure de points de parallaxe (Tie point)Mesure de points de contrôle (Check point)puis CalculOn aurait besoin « d’inverser » la formule pour en trouver les paramètres (R et S). C’est impossible à réaliser par un processus de calcul direct (système non linéaire, car les inconnues se multiplient entre elles) ; on résoudra le système en plusieurs itérations en utilisant la formule de colinéarité linéarisée au voisinage d'une solution approchée.

SOLUTION APPROCHÉE DE R ET S

● soit par un dispositif externe INS/GPS sur l'avion (assez cher !)● soit par calcul

○ en photogrammétrie terrestre : en utilisant la formule des 11 paramètres○ en photogrammétrie aérienne : en supposant le terrain plat et horizontal et l'axe de prise

de vues strictement vertical (la transformation entre le terrain et le fond de chambre devient une simple similitude plane).

LINÉARISATION DE LA FORMULE D'IMAGE

On démontre mathématiquement que l'équation de colinéarité linéarisée sd'écrit sous la forme :dm = a1 dF + a2 dM + a3 dS + a4 dR

On écrit :mobservé = mcalculé + dmoù mcalculé est la valeur calculée par la formule de colinéarité avec les valeurs initiales de R et S, tandis que mobservé , décomposé en xobservé et yobservé, contient les coordonnées image du point mesuré.Une fois tous les points d'appui mesurés, le logiciel peut donc poser un système d'équations où les inconnues sont les dF, dR, dS.

Photogrammétrie numérique – DC2 – 2008 17/29

Une fois ces variables calculées par les moindres carrés, elles sont ajoutées aux valeurs approchées de R et S, et le calcul est réitéré, et ainsi de suite jusqu'à convergence (c'est à dire quand il n'y a plus de variation entre l'itération n et l'itération n+1).

RÉSULTAT DU CALCULLe calcul est une compensation entre toutes les mesures dont les unités et les ordres de grandeur sont hétérogènes. Il faut pondérer les classes de mesures a priori pour les homogénéiser.

Application NumériqueSoit un couple d'images numériques à 40cm, de B/H 0.6. Les pondérations seront raisonnablement les suivantes :

½ pixel pour la précision image20 cm pour la précision planimétrique30 cm pour la précision altimétrique

Mal pondérer les mesures a priori (par exemple, exagérer les sigma plani et alti) donnera de plus faibles résidus terrain mais de plus fort résidus image, ce qui nuira au confort de saisie.

Comme le système est surdéterminé, la solution trouvée est celle qui « colle » au mieux à toutes les mesures. Cette solution ne colle pas parfaitement à chacune des mesures. La solution trouvée, c'est la position et l'orientation des caméras telles que les rayons perspectifs passent au plus près, à la fois des mesures image, et des mesures terrain.

Il est bon d’utiliser par ailleurs des points de contrôle: ce sont des points d’appui qui ne sont pas utilisés pour le calcul de la compensation, mais pour lesquels on calcule l’écart entre les coordonnées terrain mesurées sur le terrain et les coordonnées terrain calculées par le modèle photogrammétrique obtenu en intersectant dans l'espace à trois dimensions les rayons perspectifs passant par les mesures image et les sommets de prise de vues calculés.

Photogrammétrie numérique – DC2 – 2008 18/29

Point d'appui compensé

Point d'appui

Résidu image

Résidu terrain

VI PROBLEME SPECIFIQUE DE L'ORIENTATION D'UN COUPLE : ÉQUATION DE COPLANÉITÉ

Deux remarques préliminaires

Pendant la prise de vue, tous les rayons liant le terrain aux photos se coupent.

On va chercher une position relative des deux images telle que tous les rayons homologues se coupent. Notons que cette figure est invariante par similitude.

On peut donc séparer le problème en deux parties pour simplifier :

1. Cherchons une position relative des deux images telle que tous les rayons se coupent. 2. Appliquons une similitude à cette objet de façon à ce que le modèle réalisé corresponde au

terrain.

Le problème n°1 (Orientation relative) est un problème à cinq paramètres : les 6 paramètres de l’équation de colinéarité de chaque image, diminués des 7 paramètres de la similitude (3 pour la rotation, 3 pour la translation et 1 pour le facteur d’échelle) → 2x6 –7=5

Photogrammétrie numérique – DC2 – 2008 19/29

axe de vol

m1

S1 S2

m2

Le fait que S1m1 et S2m2 se coupent, c'est à dire que S1m1 et S2m2 sont coplanaires) s’exprime mathématiquement par le produit mixte (produit scalaire du produit vectoriel) qui s’annule :

0),,( 212211 =SSmSmS (équation de coplanéité)

Le problème n°2 (orientation absolue) est un problème à 7 paramètres (similitude dans l’espace). Il s’agit d’appliquer le modèle trouvé lors de l’orientation relative sur le terrain. On cherche simplement à faire un changement de référentiel.

+

=

Z

Y

X

M

M

M

T

T

T

TTT

ZYX

ERZYX

E étant l’échelle, R la rotation dans l’espace, et T la translation. Le système est non linéaire. Il existe deux méthodes pour le résoudre, la première s’adaptant bien à la photogrammétrie aérienne, la seconde à la photogrammétrie terrestre.

Photogrammétrie numérique – DC2 – 2008 20/29

VII AUTOMATISATION DES MESURES : LA CORRELATION AUTOMATIQUEBeaucoup d'étapes de la photogrammétrie sont basées sur la mise en correspondance de détails.

● Mesures de points de liaison (dans l'orientation relative d'un couple, l'aéro d'un bloc)● Mesures des repères de fond de chambre● Saisie 3D (qui revient à positionner le curseur sur le même détail dans l'image gauche et

l'image droite)● fabrication de MNT

On trouve schématiquement deux approches :

1. Extraction de primitives (contours, régions), et appariements (cette approche n'est pas étudiée dans ce cours)

2. Extraction de vignettes photommétriques, et mesure de ressemblance.

MISE EN CORRESPONDANCE DE VIGNETTESOn extrait un voisinage, une vignette autour du point d'intérêt. L'idée est de trouver une vignette qui lui ressemble dans l'autre image. Il faut donc se munir d’un outil permettant de dire si deux vignettes se ressemblent (critère de corrélation), c’est à dire une fonction qui aura des propriétés spéciales quand les images se ressemblent localement (attention, on ne traite pas les images en entier, car elles n’ont aucune chance de se ressembler globalement, ayant été faites de deux points de vue différents).

Remarque: plus la vignette est grande, moins il y a de chances de trouver une vignette ressemblante, mais plus elle est petite, plus on a de chance de faire de mauvais appariements (n’importe quoi peut ressembler à n’importe quoi).

On utilise habituellement le coefficient de corrélation linéaire suivant :

∑ ∑∑

−−

−−=

22 )()(

))((

yyxx

yyxxr

où x est la moyenne des radiométries sur la vignette choisie dans la première image, et y la moyenne des radiométries sur la vignette dans la deuxième image.

En se basant sur les deux développements du produit scalaire ( ∑= ii yxyx . et θcos. yxyx = ), on reconnaît que r peut être géométriquement assimilé au cosinus de l'angle θ que forment les deux vecteurs )( xx − . et )( yy − . Notons que le fait d'utiliser les vecteurs )( xx − et )( yy − .au lieu des vecteurs x et y, c'est à dire de travailler sur des valeurs centrées sur la moyenne, permet de s'affranchir d'une éventuelle différence d'éclairement entre les deux images.

Cette comparaison géométrique permet de dire que 11 ≤≤− r

Si r = -1 (les deux vecteurs sont opposés), la première vignette est le négatif de l’autre ; Si r = 0 (les deux vecteurs sont orthogonaux), les deux vignettes ne se ressemblent pas du tout, enfin si r = 1 (les deux vecteurs sont confondus), les deux vignettes sont parfaitement identiques.

Photogrammétrie numérique – DC2 – 2008 21/29

Chercher deux détails homologues revient donc à chercher r proche de 1.

Pour chercher l’homologue d’un point de l’image n°1 dans l’image n°2, il faut calculer r, le coefficient de corrélation entre V1 et toutes les vignettes V2 de même taille disponibles dans la zone de recherche.

Le problème de la corrélation est un problème mal posé :● il n'y a pas forcément de solution (pb des parties cachées)● la solution n'est pas forcément unique (cas des détails répétitifs : marquages au sol en milieu

urbain, ou rangs de vigne en zone rurale par exemple)

Il existe des techniques de validation interne d'un résultat de corrélation :1. un seuil d'acceptation (par exemple, on peut exiger qu'une vignette soit considérée

homologue si r est supérieur à un seuil, 0.8)2. un critère de morphologie de la surface de corrélation (on exigera que le maximum de

corrélation soit marqué par un pic, ce qui n'arrive pas par exemple dans les zones homogènes)

3. une fois un candidat trouvé, on le prendra comme point de départ et on cherchera son homolgue dans la première image, en espérant tomber sur le point de départ.

Le principe de recherche implique que les fausses corrélations seront particulièrement fréquentes sur les détails parallèles à la base stéréoscopique.

AUTOMATISATION DES TRAITEMENTS DE BASEL'orientation interne, qui concerne les films argentiques scannés, consiste à mesurer sur l'image les repères de fond de chambre, dont la position est par ailleurs très bien connue dans le repère "fond de chambre", pour établir le changement de repère entre l'image et le fond de chambre.

Sachant que les repères sont en général noir et blanc, de forme connue (description bitmap), et situés soit aux coins de l'image soit au milieu des côtés, il est facile d'automatiser leur mesure. En effet, la vignette de référence sera tout simplement un repère type.

La mesure des points de liaison pour l'aérotriangulation s'automatise plus ou moins facilement, selon les images, en restreignant la zone de recherche en fonction du tableau d'assemblage (recouvrements moyens ou trajectographie GPS).

Les MNT, Modèles Numériques de Terrain (description analytique du relief topographique), sont produits :

● soit en interpolant des courbes de niveau existantes● soit par saisie photogrammétrique manuelle

Photogrammétrie numérique – DC2 – 2008 22/29

V1 V2

● soit directement par levé sur le terrain (exemple GPS temps réel)● soit par laser aéroporté● soit par corrélation automatique

FORMAT DES MNTSoit une grille régulière en coordonnées terrain, soit une grille irrégulière (Triangle Irrégular Network).

Extraction des MNT par corrélation automatique (SMI, ou Stéréo Matching from Image space)

On extrait sur l'image droite, par corrélation automatique, les homolgues des vignettes centrées sur des points d'intérêt dasn l'image gauche. Les images étant géoréférencées, connaître deux détails homologues permet d'en déduire les coordonnées terrain de ce détail : il s'agit en effet de l'intersection dans l'espace à trois dimensions des rayons perspectifs passant par ces points image et les sommets de prise de vues.

Rééchantillonnage épipolairePour restreindre la zone de recherche, on pourra faire un rééchantillonnage épipolaire, c'est à dire créer de nouvelles images (même aspect que les images brutes, mais géométrie particulière).

Definition : Les droites épipolaires d'un couple sont l'intersection des deux plans image et du faisceau de plans passant par la base.

Soit d1 et d2 deux droites épipolaires. On est sûr qu'un point de d1 aura son homologue sur d2.

Le rééchantillonnage épipolaire consiste à créer de nouvelles images telles deux droites épipolaires sont parallèles au bord de l'image et de même indice. Ces nouvelles images correspondent à ce qu'aurait donné une prise de vue à axes parallèles et perpendiculaires à la base.

Attention, le rééchantillonnage épipolaire n'est possible qu'après la mise en place du couple.

L'intérêt d'un tel rééchantillonnage, est que la recherche d'un homologue se fait désormais sur une ligne, et non plus sur un ensemble de lignes.

Extraction des MNT par corrélation automatique (SMO, ou Stéréo Matching from Object space)

On part du terrain. On cherche à déterminer l'altitude d'un point de coordonnées X,Y, connaissant un intervalle d'altitudes probables. Les images étant géoréférencées, on projette par photogrammétrie tous les points de coordonnées (X,Y,Zi) dans les deux images, Zi étant la discrétisation de l'intervalle Zmin-Zmax. Pour chaque point projeté, on extrait une vignette dans les images gauche et droite, et l'on calcule le coefficient de corrélation entre ces vignettes. L'altitude Zi ayant donné meilleur coefficient de corrélation sera gardée.

Photogrammétrie numérique – DC2 – 2008 23/29

VIII ORTHOIMAGE

Images quelconques : projection perspective du terrain sur un plan

Carte : projection orthogonale du terrain sur un plan

Les différences entre la carte et l'image sont les suivantes :● la sémiologie● la projection géométrique● les variations d'échelle sur les images brutes, dues au relief et à l'inclinaison del'axe de prise

de vues● les dévers de bâtiment sur les images brutes, sachant que sur une carte on ne voit que le haut

des bâtiments

Définition : Une orthoimage est une image qui a la même géométrie que la carte, mais les couleurs de l'image. C'est une image géoréférencée, c'est à dire que tous ses pixels sont connus en coordonnées terrain. C'est une image (ou une mosaïque d'images) corrigée des déormations dues :

● au relief● à l'inclinaison del'axe de prise de vues● à la distorsion optique

Ses caractéristiques sont● emprise carto● résolution spatiale (taille du pixel sol)● le système de coordonnées de référence● résolution spectrale (panchro, RVB, IR, etc)● précision ● actualité

Il est impératif de dissocier les notions de précision et de résolution spatiale.

ExempleSoit une orthoimage de résolution 50cm, de précision 10m.Connaissant les coordonnées image (ligne, colonne) d'un point de contrôle, on compare ses coordonnées sur l'orthoimage (X = Xmin + colonne x résol / Y = Ymin + ligne x résol) avec ses coordonnées de référence. Si la précision annoncée de l'orthoimage est de 10m, on ne s'étonnera pas de trouver un tel écart.

Photogrammétrie numérique – DC2 – 2008 24/29

Réduction d’échelleRéduction d'échelle

PRODUCTION D'UNE ORTHOIMAGESachant qu'une orthoimage est une image corrigée des déformations dues au relief, à l'inclinaison de l'axe de prise de vues, et à la distorsion optique, il est impératif de connaître le relief, l'orientation de la caméra, et les paramètres intrinsèques à la caméra. Ces informations sont données par, respectivement, un MNT, la fonction de localisation de l'image brute, et le certificat d'étalonnage de la caméra.

On part de l'orthoimage à produire. L'orthorectification consiste à trouver, pour chaque pixel, quelle couleur de l'image brute utiliser pour le "peindre".

On considère un pixel de coordonnées terrain (X,Y). Grâce au MNT, on lui associe un Z. On peut alors appliquer la formule d'image, et connaître les coordonnées fond de chambre de l'image de ce point. A cette étape, on ajoute la distorsion optique aux coordonnées fond de chambre obtenues. Ensuite, si l'on travaille avec des images de caméra numérique, les coordonnées seront directement des coordonnées image (puisqu'alors le référentiel fond de chambre est identique au référentiel image). Si l'on travaille avec des clichés scannés, on transformera les coordonnées fond de chambre en coordonnées image grâce aux paramètres de l'orientation interne.

Comme les coordonnées image ainsi calculées n'ont aucune chance d'être des coordonnées "rondes", on interpole la radiométrie dans l'image brute, soit

● interpolation au plus proche voisin (on prend alors la radiométrie du pixel le plus proche du point calculé) (rapide / effet d'aliasing)

● interpolation bilinéaire (on prend une combinaison linéaire des 4 pixels les plus proches du point calculé, sachant que plus un pixel est proche du point calculé, plus sa radiométrie contribuera au résultat final) (assez rapide / effet flou)

● interpolation bicubique (on prend une combinaison des 12 ou 16 pixels les ples proches du point calculé) (long / esthétique correcte)

Photogrammétrie numérique – DC2 – 2008 25/29

(x,y) = g (X,Y,Z)(x,y) = g (X,Y,Z)x,y fond de

chambrex,y fond de

chambre

(l,c) = f(x,y)(l,c) = f(x,y)

ligne /colonne

radiométrieradiométrie

X,Y,Z terrainX,Y,Z terrain

X,Y

MNTMNT

PRODUCTION D'UNE MOSAIQUEPour couvrir une emprise donnée, soit on fait une orthorectification image par image, puis on assemble ces orthos, soit on met en entrée du processus tout un bloc d'images en définissant des règles de choix de l'image sur laquelle on ira chercher la radiométrie d'un pixel donnée (X,Y). En effet dans ce cas, si l'on part d'une prise de vues stéréoscopique, un point terrain est forcément vu sur au moins deux images.

La règle fréquemment utilisée dans ce cas est la règle du plus proche nadir. Elle consiste à prendre la radiométrie sur l'image où le point est le plus proche du nadir. Ceci évitera que les objets non pris en compte par le MNT (bâtiments par exemple) se voient trop. En effet, le schéma ci-dessous montre que les objects seront d'autant moins "étalés" sur l'image qu'ils seront proche du centre (dr est petit quand r est petit).

On s'attachera à rendre les lignes de mosaïcages les moins visibles possibles.

EGALISATION RADIOMETRIQUESur une zone étendue, les images peuvent avoir été prises à des dates différentes : elles seront alors radiométriquement différentes. Par ailleurs, sur une même mission, le phénomène de hotspot, dû à la position relative du soleil et de la caméra (l’énergie reçue par le capteur augmente lorsque la direction de visée se rapproche de la direction du soleil), crée des systématismes esthétiquement génants. Il faut corriger ce phénomène, soit par une égalisation physique, soit par des corrections pragmatiques. On procédera ensuite à un réhaussement en couleur et en dynamique.

Application numériqueTaille en Go de l’orthophoto couleur au pas de 50 cm couvrant la France métropolitaine ? 550 000 000 km2, un pixel de 0,25m2 d'où le nombre de pixels, puis le nombre d'octets.

LIVRAISON DES ORTHOIMAGESOn livre une image, ou une image par morceaux (dalles), accompagnée d'un fichier de géoférérencement pour un calage direct dans un SIG (le format de ces fichiers texte est propre à chaque logiciel, et contient généralement :le nom du fichier image, Xmin, Xmax, Ymax, Ymax, taille du pixel sur le terrain, nombre de lignes de l'image, nombre de colonnes de l'image).

Photogrammétrie numérique – DC2 – 2008 26/29

H

∆ Z

R∆ R

rdr

dr =∆ ZH

r

dr

r=

∆ RR

=∆ ZH

QUALITÉ DES ORTHOIMAGESOn distingue qualité géométrique et esthétique.La qualité géométrique est liée à l'atteinte de la précision spécifiée, la continuité des réseaux.La qualité esthétique est liée à la qualité des images en entrée, à l'égalisation radiométrique, auréhaussement, à l'absence de défaut visible.

ORTHOIMAGES VRAIES OU ORTHOIMAGES INTEGRALESUne orthoimage « vraie » est une image débarrassée des déformations dues :

● à l'inclinaison de la caméra● à la distorsion optique● au relief du sol et du sursol (bâtiments, végétation...)

Pour produire une orthoimage « vraie », il existe au moins trois solutions :● prise de vues à forts recouvrements latéral et longitudinal, de façon à ne pouvoir conserver que la partie centrale des images où les bâtiments ne sont pas trop déformés, et processusd'orthorectification classique● prise de vues par l'ADS40 (caméra à barrette de Leica), de façon à n(utiliser que l'image au nadir.● processus d'orthorectification classique, en remplaçant le MNT par un MNS (description du sol et du sursol).

FABRICATION D'UNE ORTHOIMAGE VRAIE

Attention cependant, il ne suffit pas de remplacer le MNT par un MNS dans le processus d'orthorectification. En effet, ne pas gérer explicitement les parties cachées de l'image brute pourrait générer des doublons dans l'orthoimage.

Photogrammétrie numérique – DC2 – 2008 27/29

(x,y) = g (X,Y,Z)(x,y) = g (X,Y,Z)x,y fond de

chambrex,y fond de

chambre

(l,c) = f(x,y)(l,c) = f(x,y)

ligne /colonne

radiométrieradiométrie

X,Y,Z terrainX,Y,Z terrain

X,Y

MNS

ATTENTION !

(x,y) = g (X,Y,Z)(x,y) = g (X,Y,Z)x,y fond de

chambrex,y fond de

chambre

(l,c) = f(x,y)(l,c) = f(x,y)

ligne /colonne

radiométrieradiométrie

X,Y,Z terrainX,Y,Z terrain

X,Y

MNS

(x,y) = g (X,Y,Z)(x,y) = g (X,Y,Z)x,y fond de

chambrex,y fond de

chambre

(l,c) = f(x,y)(l,c) = f(x,y)

ligne /colonne

radiométrieradiométrie

X,Y,Z terrainX,Y,Z terrain

X,Y

MNS

ATTENTION !

En effet, dans le processus d'orthorectification, un pixel de l'ortho en sortie peut, combiné au MNS, avoir les coordonnées (X1,Y1,Z1) L'utilisation de la formule d'image ne permet pas de savoir que ce point n'est pas vu sur l'image. Si on n'ajoute aucun test, le logiciel mettra dans ce pixel la radiométrie du point haut, de coordonnées (X2,Y2,Z2) qui se projette photogrammétriquement au même endroit sur l'image brute (voir schéma ci-dessus).

Pour gérer les occlusions, on utilise un Z-buffer :

Image brute : chaque pixel contient une radiométrie

Z-buffer : image de mêmes dimensions que l'image brute, dont les pixels sont remplis, non pas par des radiométries, mais par les Z du MNS projetés photogrammétriquement dans l'image. Quand deux points se projettent au même endroit (cas des points 1 et 2 dans l'exemple ci-dessus), on ne garde que le Z le plus haut.

Photogrammétrie numérique – DC2 – 2008 28/29

X2,Y

2,Z

2

X1,Y

1,Z

1

On inclut ensuite, dans la boucle d'orthorectification, après avoir appliqué la formule d'image, un test qui contrôle que le Z du point considéré n'est pas plus petit que le Z stocké dans le Z-buffer à cet endroit. Le cas échéant, on ne remplit pas le pixel par une radiométrie, mais par du noir. Pour livrer l'orthoimage finale, on combinera plusieurs orthoimages où les parties cachées se complètent.

Photogrammétrie numérique – DC2 – 2008 29/29

(x,y) = g (X,Y,Z)(x,y) = g (X,Y,Z)x,y fond de

chambrex,y fond de

chambre

(l,c) = f(x,y)(l,c) = f(x,y)

ligne /colonne

radiométrieétrie

X,Y,Z terrainX,Y,Z terrain

X,Y

MNS

Z du Zbuffer

1 2

3

45

Si Z du Zbuffer > Zterrain

6

Si Z du Zbuffer < Zterrain

(x,y) = g (X,Y,Z)(x,y) = g (X,Y,Z)x,y fond de

chambrex,y fond de

chambre

(l,c) = f(x,y)(l,c) = f(x,y)

ligne /colonne

radiométrieétrie

X,Y,Z terrainX,Y,Z terrain

X,Y

MNS

Z du Zbuffer

1 2

3

45

Si Z du Zbuffer > Zterrain

6

Si Z du Zbuffer < Zterrain