Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján Patkó Tamás Budapest, 2002. november
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
1
A cikk célja egy példán keresztül bemutatni a digitális sztereó képpár kiértékelés lehetséges lépéseit és azok felhasználását. A cikkben bemutatom a légifelvétel feldolgozását, a kamerák pozíciójának meghatározását valamint egy algoritmikus megoldást az egymáshoz tartozó pontpárok meghatározására. A pontpárok meghatározásánál bemutatom a sztereó párosítási problémák alapeseteit, valamint ezek kiküszöbölésére alkalmazott megoldásokat.
1 Bevezetés Az 1980-as évektől a számítástechnika rohamos fejlődése lehetővé tette, hogy az egyre gyorsabb és egyre nagyobb kapacitású számítógépek segítségével lehetővé váljon az olyan problémák megoldása is, mint a sztereó képpárok feldolgozása és a háromdimenziós tér rekonstrukciója képek alapján. Első lépésként a számítógépek még csak a hagyományos analóg térkiértékelő műszerek szerepét vették át. Az automatizmusok
beépítése
a
számítógépes
rendszerekbe
napjainkban
is
folyamatosan tart. Az általam kitűzött cél egy olyan algoritmus kifejlesztése volt, ami lehetővé teszi a képek alapján az automatikus térrekonstrukciót minimális humán beavatkozás igénybevételével. A probléma megoldásához a következő lépések elvégzésre van szükség: - kép belső tájékozásának és paramétereinek meghatározása - kamerapozíciók meghatározása - pontpárok keresése és a háromdimenziós koordináták meghatározása - igényelt térmodell elkészítése A cikk keretein belül az első három lépést fogom részletesen tárgyalni.
2 A kiindulási adatok A cikkben bemutatásra kerülő képek a Hexium Kft. megbízásából készültek. A képek készítéséhez használt kamera típusa: Wild RC-30, a kamera fókusza: 152.70 mm, képméret: 23x23 cm, repülési magasság: kb. 1500 m, képméretarány: kb. 1:10 000, a felvétel időpontja: 2000. Szeptember. A felvételt készítő repülő nagyjából északdéli irányban repülve készítette a két képet. Az először készített kép neve 1084.TIF (bal), a másodszor készített képé 1086.TIF (jobb).
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
2
3 A képek feldolgozása 3.1 A képek belső tájékozása A kameramodellhez szükséges a képek un. belső tájékozása. A belső tájékozáskor meghatározásra kerül a kép középpontja és fizikai felbontása. A képek belső tájékozását a kamera hitelesítési dokumentációja alapján lehet elvégezni. A hitelesítési dokumentáció tartalmazza a keretjelek pontos pozícióját, a kép szimmetria
és
autokollimációs
középpontjának
koordinátáit
és
a
kamera
fókusztávolságát. [Kraus98] A keretjelek elhelyezkedését a következő képek mutatják:
A keretjelek számozása a bal felső saroktól indulva az óra járásával megegyező irányban: 3, 7, 4, 8, 1, 5, 2, 6 A keretjelek koordinátái: 1 2 3 4
X (mm) 106.001 -106.007 -105.998 106.004
Y (mm) -106.000 -106.004 105.997 106.001
5 6 7 8
X (mm) 0.001 -109.997 0.006 110.002
Y(mm) -109.999 -0.003 109.998 0.005
A középpontok koordinátái: (PPA:autokollimációs, PPS szimmetria) PPA PPS
X (MM) -0.001 0.001
Y (MM) -0.005 -0.009
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
3
A kamera fókusztávolsága: 152.700 mm A szkennelt képek mérete kb. 15000x15000 pixel. A fenti adatokból átszámolva a pixelméret 15 µm. A kép pixelméretéből adódik, hogy az autokollimációs és szimmetria-középpontok a pixeles képen ugyanarra a pixelre esnek. (A távolságuk kb. 4.47 µm) A kamera kalibrációja megadja a négy átlós irányban a kamera elrajzolásának mértékét is µm-ben. Sugár mm 10 20 30 40 50 60 70 80 90 100 110 120 130 140 148
1 0.0 -0.7 -1.6 -2.2 -2.4 -2.8 -2.1 -2.2 -0.5 0.0 1.0 2.1 2.2 2.3 2.8
Irány 3 -0.5 -0.7 -1.3 -1.4 -2.0 -1.9 -1.7 -0.6 0.0 0.7 1.6 1.8 1.6 1.8 1.8
Átlag 2 -0.1 -0.8 -1.2 -1.9 -2.3 -1.8 -1.8 -1.4 -1.0 0.1 0.3 0.1 0.7 0.8 1.2
4 -0.2 -0.4 -1.4 -1.6 -1.9. -2.3 -1.9 -1.4 -0.8 0.0 0.8 0.9 1.0 0.5 0.3
-0.2 -0.6 -1.3 -1.7 -2.1 -2.2 -1.8 -1.4 -0.5 0.2 0.9 1.2 1.3 1.3 1.5
Mivel az elrajzolás mértéke nem éri el a pixel méretét, ezért a további számolásoknál nem kell figyelembe venni. A képekre számolt adatok: Kép
Méret (pixel)
Középpont (pixel)
1084
15008,14976
7496,7544
1086
14976,15072
7500,7557
A fenti adatokkal a kép belső tájékozása megtörtént. Az adatok birtokában a kameramodell minden szükséges paramétere ismert.
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
4
3.2 A képek külső tájékozása – kamerapozíciók meghatározása A kamerák pozícióinak meghatározására több módszer ismert az irodalomban. Az adatok meghatározása történhet képenként vagy a képpárra. A képenkénti paraméter meghatározáshoz elengedhetetlen legalább három a világkoordinátarendszerben ismert illesztőpont. Mivel az illesztőpontok megadása hibával terhelt, ezért érdemes minél több illesztőpontot felhasználni a képen egyenletesen elosztva. A képpárokat felhasználó pozíció-meghatározás is kétféle lehet. Az egyik esetben nincs ismert illesztőpont csak a képeken adott pontpárok. Ebben az esetben sajnos a kamerapozíció nem határozható meg egyértelműen ezért mindenképpen szükséges további geometriai információ. [Jankó02] [MVGeom02] A másik esetben legalább három illesztőpont és további pontpárok felhasználásával történik meg a kamerapozíció-meghatározás. Ebben az esetben az illesztőpontok hibájának csökkentésére lehet az ismeretlen koordinátájú pontpárokat felhasználni. 3.2.1 A kamerapozíció meghatározása ismert illesztőpontok alapján. Az egyszerűség kedvéért én az egyképes estet implementáltam. Ha X0,Y0,Z0 az adott illesztőpont világkoordinátái, C'x, C'y, C'z a kamera elforgatott világkoordinátái és R a kamera rotációs mátrixa, akkor egy adott pont kamera koordinátarendszerben vett koordinátáira a következő összefüggés írható fel:
ªX º ª «Y » = « « » « «¬ Z »¼ «¬
R
º ª X 0 º ªC `x º » « Y » − «C `y » »« 0 » « » »¼ «¬ Z 0 »¼ «¬ C `z »¼
, ha az illesztőpont koordinátáinak elforgatottját Rx, Ry, Rz-nek hívjuk, továbbá a centrális vetítés szabályaként felírjuk a képkoordinátákat: x= f
X Y ,y= f Z Z
,ahol f a kamera fókusztávolsága, akkor adott R forgatási mátrix esetén C' kamerapozícióra
egy
túlhatározott
egyenletrendszert
kapunk.
A
világ-
koordinátarendszerben lévő kamerapozíciót C` és az R vektor inverzének szorzataként kapjuk. C=R-1C`
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
5
A programmal a fenti egyenletet négy különböző módon kiszámoltam. (Különböző részeit felhasználva az egyenletrendszernek.) Az eredménypozíciókat feljegyeztem minden lehetséges illesztőpont párosra. Kiszámítottam az így kapott pozíciókra a hibavektorokat. Hibavektornak azt a vektort tekintettem, amit az adott C ponthoz hozzáadva megkapom az adott illesztőponthoz tartozó kameraközéppontot. A teljes hibavektor ezen hibavektorok átlaga. Az algoritmus célfüggvénye, hogy ez az átlag hibavektor legyen minimális. Az algoritmus adott határok közt négy lépésben egyre finomodó módon végigmegy az összes forgatási szögön, és közben keresi a minimális hibavektorral rendelkező kameraközéppontot. A futási eredmény:
Az eredmény a ’legjobb:’ sorban van a következő sorrendben: α, β, γ, Cx, Cy, Cz,
ahol
α, β, γ a tengelyek körüli forgatás fokokban. A koordinátaértékek
méterben értendők. A kameraparaméterek alatti sorban a hiba növelt értékének (+0.99) négyzete, az átlagos hiba és a maximális hiba értéke látható. (Az algoritmus a növelt érték négyzetét minimalizálta.) A fenti értékekből látszik, hogy a bal kamera pozíciója átlagosan 30 cm hibával, míg a jobb kamera pozíciója kb. 46 cm hibával került meghatározásra. Ez a további eredmények szempontjából additív hibának tekintendő.
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
6
3.2.2 A felhasznált illesztőpontok Az előző számolásnál használt
illesztőpontok
meghatározása
egy
négy
tesztterületből álló mérési sorozat eredményeképpen jött ki. A mérések GPS helymeghatározó készülék segítségével készültek centiméteres pontossággal. A mért pontok AutoCad DXF állományba kerültek. Az illesztőpontok meghatározásakor a DXF állományban olyan pontokat kellett keresni, melyek jól beazonosíthatóak mindkét képen. Példa az illesztőpontra: (út menti mérföldkő teteje)
3.3 Pontpárok meghatározása A képpontok háromdimenziós koordinátáinak meghatározásához a következő feladat az egymásnak megfelelő pontpárok megkeresése. Ha az összes pont párját meg kívánjuk keresni, akkor az még a mai számítógépeknek is rengeteg időbe telne, ezért a pontok halmazát szűkíteni kell. Figyelembe kell venni, hogy a pontok száma képenként kb. 225 millió. 3.3.1 Potenciálisan párosítható pontok szűrése Első lépésként ki kell szűrni azon pontokat, melyek párosítása eleve reménytelen feladat. Ezt az adott képekre olyan programmal végeztem, ami megvizsgálta a későbbiekben leírt pontillesztő algoritmussal, hogy az adott képen egymás mellett lévő pontok mennyire hasonlóak. Ha két egymás melletti pont hasonló, akkor a feltételezés az, hogy a későbbi algoritmus úgysem tudna dönteni köztük, tehát a mérés eredménye nagy hibával terhelt lenne. Ebből kiindulva az egymás melletti hasonló pontokat eleve kiszűrtem a párosításból. Ezzel a módszerrel az alábbi
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
7
képrészlet
esetében
a
pontoknak
csak
46.5%-a
maradt
meg
a
további
számításokhoz.
A fenti bal oldali képen a Rudas fürdő -bal képből kivágott- részlete látható, a jobb képen ennek maszkolt verziója. A képen jól kivehető, hogy az úttest, a tető valamint a Duna pixelei kiesnek a további vizsgálatokból, mivel „önhasonlók”. Egy teljes képmaszk (225 millió pixel) előállításának ideje: kb. 3 óra. (Pentium III 1200MHz 256Mbyte RAM) 3.3.2 Az összehasonlító algoritmus Az összehasonlító algoritmus célja két pont és környezetük összehasonlítása a következő ismérvek alapján: -
a számolás végeredménye egy pozitív szám mely hasonló pontok esetén közelít nullához - az összehasonlítás invariáns kell, hogy legyen a két kép közti intenzitáskülönbségekre - az összehasonlítás eredményéből lehetőleg következtetni lehessen a fantompárosításokra - lehetőleg gyorsan számolható legyen Ennek megfelelően egy súlyozott különbségi értéket képeztem az alábbiak szerint:
(
) (
Score = W x , y PLeft , x , y − PLeft − PRight , x , y − PRight
)
x, y
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
8
,ahol W xy a súlyfüggvény értékei az alábbi ábra szerint. 300 250 200 150 100 50 0
Az ablakméret 32x32 pixel. A fenti összefüggéssel és súlyfüggvénnyel a képek közti különbségek ellenére jó párosítást értem el. A súlyfüggvény lényege, hogy az adott pixel domináns, de a környezete is meghatározza a párosítást. 3.3.3 A pontpár keresése A pont párjának keresésénél tovább szűkíthető a keresés köre, ha felhasználjuk az epipoláris geometriából adódó megszorításokat, tehát a pont párját a másik képen csak ott keressük, ahol a vetítésből adódóan elhelyezkedhet. Ez torzításmentes esetben egy egyenes. Az egyenes koordinátáit a következő módon kaphatjuk meg: Vegyük az adott bal oldali pont térbeli –a tárgy felé néző- vektorát úgy, hogy a kamera tengelye felé eső vetülete legyen 1 hosszú: ª xº «f » « L» y VL = « » , ahol x,y a pont koordinátái a képen, fL a bal kamera fókusztávolsága. « fL » « − 1» « » ¬ ¼ Ezután annak az egyenesnek az irányvektorát, mely a bal kamera vetítési középpontjából a képponton keresztülmegy a kamera rotációs mátrixának inverzével való szorzással kaphatjuk meg. L = RR−1VL Ahhoz, hogy a jobb oldali képen megkapjuk az egyenest, ennek a térbeli egyenesnek a képe kell. Ezt a fenti adatok alapján már elő lehetne állítani, de ezen a ponton újabb
szűkítést
tehetünk
a
kereséssel
kapcsolatban.
Tegyünk
egy
olyan
megszorítást, hogy a kamerától csak adott távolságtartományon belüli pontokat kívánunk meghatározni. Ekkor az egyenes két végpontjának koordinátáit a következő térbeli pontok vetítésével kapjuk:
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
9
PNear = C L + LDMin , PFar = C L + LD Max , ahol CL a bal kamera vetítési középpontja, DMin és DMax a távolságlimit értékek. A később bemutatott példákban a minimum és maximum értékek 1400 és 1600 méter. A jobb oldali képen a pont képét az alábbi összefüggés adja meg:
PR = RR ( P − C R ) , ahol RR a jobb kamera forgatási mátrixa és CR a jobb kamera vetítési középpontja.
x=
PR , X − PR , Z
f R és y =
PR ,Y − PR , Z
f R , ahol a Z koordinátával való osztásnál a –1 a Z tengely
irányának ellentétes állása miatt kell. (A tárgy felé néz.) 3.3.4 A pontpár meghatározása Miután megvan az az egyenes amin a pont párjának lennie kell, már csak azt kell eldönteni, hogy van-e a pontnak párja, és ha igen melyik az. Természetesen a pontpár keresésénél a jobb képen kimaszkolt pontokat is ki lehet hagyni a vizsgálatból tovább csökkentve a számításigényt. Felmerülő problémák: -
a pontnak nincs párja a különböző kamerapozíciók miatt (takarásban van) a pont párja a megváltozott fényviszonyok miatt nem hasonlít (idő telt el a két kép között) - a pont párja nem az egyenesen van, mivel a tárgy elmozdult (pl. autók) - a vonal mentén több találat van (periodikus objektum) Szerencsére a pontpárosító algoritmus igen jellegzetes eredményt ad jó párosítás esetén: (Bal 2673,1519)
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
10
Pár nélküli pont esetén: (Bal 2612,1502)
Elmozdult tárgyon lévő pont: (Bal 2732,1398)
Többszörös találat, ismétlődés miatt: (Bal 2031,2783)
A mellékelt párosítási eredmények jellemzően mutatják, hogy az adott pontnak mikor van párja a másik képen.
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
11
A fentiek alapján a következő megszorítást lehet tenni: ha az abszolút minimumhoz képest létezik olyan lokális minimum, mely az abszolút minimum 120%-nál kisebb, akkor a pontnak nem lehet biztosan megtalálni a párját. A párosítás eredményében található „fésű” képződmény abból adódik, hogy a program a vonal ±6 pixeles környezetében végzi a keresést. Abból adódóan, hogy a pont párja így nem feltétlenül helyezkedik el a vonalon, további mérési hiba adódik. A program a pont háromdimenziós koordinátáit úgy számolja ki, hogy -kitérő egyenesek esetén- az egyenesek közti legrövidebb szakasz felezőpontját adja meg. Hibaértéknek visszaadja a két egyenes távolságának felét. 3.4 Régió feldolgozása A fenti algoritmusok használatával lehetőség van adott régiók vagy a teljes kép feldolgozására. A feldolgozás eredményeképpen az alábbi adatok képződnek: -
bal kép képpont koordináták bal kép szín jobb kép képpont koordináták jobb kép szín világkoordináták (3D) mérési hiba érték Az így számolt adatokból a kívánalmaknak megfelelő felület illetve térmodell előállítható. Előnye ezen megadásnak, hogy a nem számolható pontok nem szerepelnek az eredményben, valamint az eredmény minden pontjához tartozik egy hibaérték is. A számolási teljesítmény érzékeltetésére a mintaterület futási eredménye: - régió koordinátái: 2300.1364-3200,2000 a széleket is beleértve - a régió pontjainak száma: 573937 (100%) - maszkolás utáni pontok száma: 266915 (46.5%) - számolási idő: 10 óra 7 perc (607 perc) Pentium III 1200MHz 256MB RAM - kiszámolt pontok száma: 166530 (az eredeti 29%-a, a maszkolt rész 64.4%-a) - sebesség: 4.572 pont/mp az eredmény szempontjából - sebesség: 7.329 pont/mp a maszkolt részen Minta a kimeneti állományból:
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
12
LeftX
LeftY
LeftVal
RightX RightY RightVal
X
2300
1364
69
2530
4264
87
649989,95 238360,40
Y
Z
Error
110,72
0,13
2300
1365
79
2530
4265
96
649989,94 238360,25
110,72
0,13
2300
1366
88
2530
4266
108
649989,92 238360,10
110,72
0,13
2300
1367
101
2530
4267
118
649989,91 238359,95
110,72
0,13
2300
1368
107
2530
4268
128
649989,89 238359,80
110,71
0,14
2300
1369
115
2530
4269
138
649989,88 238359,65
110,71
0,14
2300
1370
120
2530
4270
141
649989,86 238359,50
110,71
0,14
2300
1371
123
2530
4271
144
649989,85 238359,35
110,71
0,14
2300
1372
125
2530
4272
144
649989,83 238359,20
110,71
0,14
2300
1373
124
2530
4273
143
649989,81 238359,05
110,71
0,15
2300
1374
124
2531
4274
143
649989,84 238358,93
110,63
0,08
2300
1375
124
2531
4275
138
649989,82 238358,78
110,63
0,08
2300
1376
122
2531
4276
134
649989,81 238358,63
110,63
0,08
2300
1377
119
2531
4277
131
649989,79 238358,48
110,63
0,08
2300
1378
114
2531
4278
128
649989,78 238358,33
110,63
0,09
2300
1379
111
2531
4279
129
649989,76 238358,18
110,63
0,09
2300
1380
109
2531
4280
127
649989,75 238358,03
110,63
0,09
2300
1381
107
2531
4281
126
649989,73 238357,88
110,63
0,09
2300
1382
107
2530
4281
126
649989,45 238358,01
110,21
0,16
2300
1538
109
2537
4452
105
649990,75 238330,41
116,87
0,02
2300
1543
124
2534
4452
105
649989,41 238331,07
114,65
0,22
2300
1544
123
2534
4452
105
649989,17 238331,21
114,17
0,21
2300
1545
123
2534
4453
104
649989,15 238331,06
114,17
0,22
2300
1546
121
2533
4454
103
649989,10 238330,89
114,24
0,29
2300
1547
120
2533
4455
104
649989,08 238330,74
114,24
0,29
2300
1548
121
2533
4456
105
649989,07 238330,59
114,24
0,29
2300
1549
121
2533
4457
109
649989,05 238330,44
114,24
0,29
2300
1557
172
2542
4519
185
650001,44 238313,91
139,40
0,13
2300
1558
176
2542
4520
188
650001,43 238313,76
139,40
0,13
2300
1559
179
2542
4521
191
650001,41 238313,61
139,40
0,12
Az eredmény szürkeskálás magassági színekre leképezve: (A pont magasságával egyenesen arányos a kép fényessége. Kék pontok jelölik a nem meghatározott pontokat. A kép a nyomtathatóság érdekében nemlineáris módon világositott.)
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
13
4 Összefoglalás A cikkben egy példán keresztül bemutattam a sztereo-légifelvételek automatikus kiértékelésének egy módját. A fenti eredmények lehetővé teszik a felvételekből nyert további adatok feldolgozását térképészeti, rádió hullámterjedési vagy egyéb feladatoknak megfelelően. A
további
használhatóság
biztosítása
érdekében
szükséges
a
számolás
felgyorsítása annyira, hogy az egész kép feldolgozása elfogadható időn belül megtörténjen. Tudván, hogy néhány évente egész Magyarországot törvényi szabályozás alapján végigfényképezik, ezért a képek feldolgozásának idejét olyan kicsire kell szűkíteni, hogy a következő fényképezés előtt még fel lehessen használni az eredményeket. Irodalomjegyzék: Kraus98
Karl Kraus: Fotogrammetria, ISBN 963 85129 9 7
Tertia
kiadó
Budapest
1998
Szirmay95
L. Szirmay-Kalos: Theory of three-dimensional computer graphics, Akadémiai kiadó 1995 ISBN 963 05 6911 6
Jankó02.
Jankó Zsolt: Háromdimenziós színtér-rekonstrukció alapvető algoritmusai, Képfeldolgozók és alakfelismerők III. konferenciája, Domaszék 2002. január 23-25, pp.10-20
MVGeom02: Richard Hartley and Andrew Zisserman: Multiple View Geometry in computer vision, Cambridge University Press, 2002. (reprinted with corrections) További irodalom: Sárközy Ferenc: Geodézia, Tankönyvkiadó, Budapest, 1984. Karl Kraus, Josef Jansa, Helmut Kager: Photogrammetry (Volume 2), Dümmler, Bonn, 1997. Karl Kraus: Photogrammetrie (Band 3), Dümmler, Köln, 2000. Richard I. Hartley: In Defense of the Eight-Point Algorithm, IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 19, No. 6. June 1997. S. Hati, S Sengupta: Robust camera parameter estimation using genetic algorithm, Pattern Recognition Letters 22 (2001) 289-298. Marc Pollefeys: Self-Calibration and Metric 3D Reconstruction from Uncalibrated Image Sequences, Katholieke Universiteit Leuven, 1999.
Háromdimenziós térrekonstrukció sztereó légifelvétel képpár alapján.
14