Diplomová práce
Delaunayho triangulace a její aplikace vypracoval: Filip Leifer vedoucí práce: Doc. RNDr. Ing. Miloš Šeda, Ph.D. obor: Inženýrská informatika a automatizace specializace: Informatika 2006
ANOTACE V této práci je rozebrána problematika generování Delaunayho trojúhelníkových sítí. Jsou zde uvedeny jednotlivé algoritmy pro konstrukci Delaunayho triangulace a některé příklady aplikací Delaunayho triangulace na vybrané problémy. Dále rozbor vybraného problému a jeho implementaci spolu s implementací algoritmu pro generování Delaunayho triangulace. ANNOTATION In this thesis are mentoined problems of generating Delaunay´s triangular meshes. There are described individual algorithms for construction Delaunay´s triangulation and some example of application of Delaunay´s triangulation on selected problems. In addition, an analyse selected problem and its implementation together with an implementation of algorithm for generating Delaunay´s triangulation.
Děkuji vedoucímu této diplomové práce Doc. RNDr. Ing. Miloši Šedovi, Ph.D za jeho vstřícný přístup a poskytnuté podklady pro dokončení této práce.
Prohlášení: Prohlašuji, že jsem diplomovou práci vypracoval samostatně a výhradně s použitím uvedené literatury. V Brně dne 29.5.2006 Podpis
Delaunayho triangulace a její aplikace
1
ÚVOD .................................................................................................................. 13 1.1 1.2
2
Bod ............................................................................................................... 15 Přímka .......................................................................................................... 15 Konvexní geometrický útvar........................................................................ 15 Eukleidovská vzdálenost .............................................................................. 15
DELAUNAYHO TRIANGULACE A SOUVISEJÍCÍ POJMY .................... 17 3.1 3.2 3.3
4
Problém triangulace...................................................................................... 13 Cíle práce...................................................................................................... 13
ZÁKLADNÍ GEOMETRICKÉ POJMY ......................................................... 15 2.1 2.2 2.3 2.4
3
Strana 11
Konvexní obálka .......................................................................................... 17 Voroného diagram........................................................................................ 17 Delaunayho triangulace................................................................................ 21
GEOMETRICKÉ ALGORITMY .................................................................... 25 4.1 Algoritmy pro konstrukci konvexní obálky ................................................. 25 4.1.1 SlowConvexHull ...................................................................................... 25 4.1.2 GHull........................................................................................................ 25 4.2 Algoritmy pro konstrukci Delaunayho triangulace...................................... 26 4.2.1 Lokální zlepšování ................................................................................... 26 4.2.2 Inkrementální konstrukce......................................................................... 26 4.2.3 Inkrementální vkládání............................................................................. 28 4.2.4 Guibas–Stolfi............................................................................................ 34 4.2.5 Delaunay Wall.......................................................................................... 38
5
APLIKACE DELAUNAYHO TRIANGULACE............................................ 41 5.1 Výšková interpolace..................................................................................... 41 5.2 Eukleidovský Steinerův minimální strom .................................................... 43 5.2.1 Základní pojmy ........................................................................................ 43 5.2.2 Vlastnosti Steinerova minimálního stromu .............................................. 44 5.2.3 Heuristický algoritmus založený na Delaunayho triangulaci................... 45 5.3 Geometrická podobnostní metrika pro případově založené usuzování........ 46 5.3.1 Reprezentace problému ............................................................................ 46 5.3.2 Geometrický problém............................................................................... 47 5.3.3 Přesný algoritmus pro nalezení optimální cesty....................................... 48 5.3.4 Rychlý algoritmus pro nalezení vhodné cesty.......................................... 49
6
IMPLEMENTACE VYBRANÉHO PROBLÉMU ......................................... 53 6.1 Algoritmus inkrementální konstrukce a jeho časová složitost ..................... 53 6.2 Výpočet acw orientace trojice bodů ............................................................. 54 6.3 Výpočet Steinerova bodu ............................................................................. 54 6.4 Datové struktury........................................................................................... 55 6.4.1 Třída Vertex ............................................................................................. 55 6.4.2 Třída Edge ................................................................................................ 56 6.4.3 Třída Triangle........................................................................................... 57 6.4.4 Třída List .................................................................................................. 58 6.5 Problém v návrhu algoritmu inkrementálního vkládání............................... 58 6.6 Úprava heuristického algoritmu pro řešení EStMT ..................................... 60 6.7 Výsledky získané upraveným algoritmem pro určení EStMT ..................... 61 VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 12
7.
ZÁVĚR................................................................................................................ 63
8.
LITERATURA ................................................................................................... 65
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
1
Úvod
1.1
Problém triangulace
Strana 13
S pojmem triangulace se můžeme setkat v mnoha vědních oborech a technických aplikacích jako je například počítačová grafika a geometrie, matematika, geodézie, simulace a dalších. Jedná se o problémy v nichž řešení spočívá na použití trojúhelníků či sítě trojúhelníků. Například se můžeme setkat s problémem jak vizualizovat naměřená data ať již se jedná o skalární veličiny či body v rovině nebo v prostoru a provést odvození hodnot v místech které měřeny nebyly. Pokud mají data pravidelnou strukturu která je nám známa, není problém odvodit hodnoty okolních bodů či body pospojovat do trojúhelníkové sítě. Avšak v případě kdy jsou data rozptýlena a nemají pravidelnou strukturu nastává problém jak tyto body propojit tak aby takováto síť vznikla. Triangulace se také používá pro teselaci polygonů, kde se jedná o vytvoření trojúhelníkové sítě uvnitř obecně nekonvexního polygonu za účelem dalšího grafického zpracování. Může být také použita jako aproximace úplného grafu který je vhodnější pro řešení některých matematických algoritmů jako je například algoritmus minimální kostry grafu. Speciálně Delaunayho triangulace, které je tato práce věnována má ještě další důležitou vlastnost, která nám zaručuje lepší výsledky při řešení problémů než jakákoliv jiná triangulace.
1.2
Cíle práce
Cílem této práce je především uvedení do problematiky triangulace a souvisejících pojmů jako je konvexní obálka a Voroného diagram. Dále popis vybraných algoritmů pro konstrukci Delaunayho triangulace a rozbor jejich principu a časové složitosti. V druhé části je to uvedení některých aplikací Delaunayho triangulace pro řešení technických problémů. V poslední části je to samotná implementace a rozbor implementace vybraných algoritmů pro aplikaci, která je řešena jako součást Diplomové práce.
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
2
Základní geometrické pojmy
2.1
Bod
Strana 15
Bod je bezrozměrný geometrický útvar. Označuje se obvykle velkým tiskacím písmenem. Obecně je bod dán N souřadnicemi, kde N je rozměr prostoru. Pro případ bodu v rovině (tj. dvourozměrný případ) je bod pevně dán dvěma souřadnicemi. Pro trojrozměrný případ je bod určen třemi souřadnicemi. Obvykle se pro označení souřadnic bodů používá označení dle názvů souřadnicových os. Tedy pro bod v rovině použijeme souřadnice (x, y) a pro bod v prostoru (x, y, z). Dále v textu i v softwarovém projektu budeme tedy zásadně nebude–li řečeno jinak označovat body velkými písmeny.
2.2
Přímka
Přímka je v matematickém slova smyslu množina nekonečně mnoha bodů. Označujeme ji malými písmeny a patří do geometrických útvarů. O bodech náležejících přímce říkáme, že jsou s přímkou ve vztahu incidence. Speciálně pokud zvolíme na přímce dva libovolné body A, B potom dostáváme úsečku AB . Krajní bod který nenáleží úsečce AB leží na některé z polopřímek které jsou dány body A a B. Dále v textu pokud nebude řečeno jinak budeme úsečky (hrany) označovat malým písmenem nebo zápisem AB .
2.3
Konvexní geometrický útvar
Geometrický útvar se nazývá konvexní, jestliže má tuto vlastnost: Jsou–li X, Y libovolné body množiny M, je každý bod úsečky XY také bodem množiny M (úsečka je podmnožinou množiny M). Za konvexní útvar se pokládá také jednobodová množina a prázdná množina. Obr. 1 znázorňuje konvexní a nekonvexní útvar.
Obr.1 Vlevo konvexní útvar, vpravo nekonvexní útvar.
2.4
Eukleidovská vzdálenost
Dále v textu budeme často používat pojem Eukleidovská vzdálenost, která je pro dva libovolné body A, B v rovině definována vztahem (1).
dist ( A, B) = ( Ax − B x ) + ( B x − B y )
(1)
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 17
3
Delaunayho triangulace a související pojmy
3.1
Konvexní obálka
Jak již bylo uvedeno, obecně je podmnožina bodů M roviny nazývána konvexní jestliže jakákoliv přímka určená dvojicí libovolných bodů z této podmnožiny je celá obsažena v M. Konvexní obálka označme ji jako CH(M) množiny M je nejmenší konvexní množina obsahující všechny body M. Pro konvexní obálku konečné množiny M bodů v rovině tedy platí, že CH(M) je unikátní konvexní mnohoúhelník jehož vrcholy jsou body z množiny M a tento mnohoúhelník obsahuje všechny body množiny M. Konvexní obálku množiny bodů je možno vidět na obr. 2 kde body konvexní obálky jsou zobrazeny černě, ostatní šedě.
Obr.2 Konvexní obálka množiny bodů v rovině
3.2
Voroného diagram
Definujme Voroného diagram pro množinu bodů V jako rozdělení roviny do n buněk, jedna pro každý bod z V s vlastností, že libovolný bod Q leží v buňce odpovídající bodu Pi tehdy a pouze tehdy jestliže dist (Q, Pi ) < dist (Q, Pj ) pro každý bod kde Pj ∈ V , j ≠ i . Označme Voroného diagram množiny bodů V jako Vor(P). Buňka Vor(P) diagramu jenž koresponduje s bodem Pi je označena jako V(Pi) a nazveme ji jako Voroného buňka pro bod Pi. Podíváme se blíže na vlastnosti jediné Voroného buňky. Pro dva body P a Q v rovině definujeme půlící osu bodů P a Q jako přímku kolmou na úsečku danou těmito body. Tato osa rozděluje rovinu na dvě poloroviny. Označme otevřenou polorovinu jenž obsahuje P jako h(P,Q) a otevřenou polorovinu jenž obsahuje Q jako h(Q,P). Poznamenejme, že libovolný bod R ∈ h( P, Q ) jen tehdy a pouze tehdy jestliže dist ( R, P ) < dist ( R, Q ) . Z tohoto dostáváme platnost rovnice (2). B
B
B
B
B
B
B
B
V ( Pi ) = ∩1≤ j ≤ n , j ≠i h( Pi , Pj )
(2)
Tedy V(Pi) je průnik n-1 polorovin a z toho důvodu otevřená konvexní mnohoúhelníková oblast ohraničená nejvíce n-1 vrcholy a nejvíce n-1 hranami. Tedy každá B
B
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 18
buňka diagramu je průnik několika polorovin, takže Voroného diagram je rovinné rozdělení jehož hrany jsou úsečky a některé hrany jsou polopřímky.
Obr.3 Voroného buňka pro bod v rovině Věta 1 Nechť V je množina n bodů v rovině. Pokud všechny body jsou kolineární potom Vor(V) obsahuje n-1 rovnoběžných přímek. Jinak je Vor(V) spojitý a jeho hrany jsou buď úsečky nebo polopřímky.
Důkaz. První část věty není těžké dokázat, takže předpokládejme, že ne všechny body jsou kolineární. Nejdříve ukážeme, že hrany Vor(V) jsou buď úsečky nebo polopřímky. Již víme, že hrany Vor(V) jsou části přímek, jmenovitě části půlicích os mezi páry bodů. Předpokládejme pro důkaz sporem, že existuje hrana e z Vor(V) která je přímka. Nechť e je na hranici Voroného buněk V(Pi) a V(Pj) (viz. obr. 4). Nechť Pk ∈ V je bod který není kolineární s Pi a Pj. Půlící osa bodů Pj a Pk není rovnoběžná s hranou e a proto protíná hranu e. Ale potom část hrany e jenž leží uvnitř h(Pk,Pj) nemůže být na hranici V(Pj), protože to je blíž k Pk než k Pj a to je spor. To je k důkazu, že Vor(P) je spojitý. Jestliže se nejedná o tento případ, potom by zde mohla být Voroného buňka V(Pi) rozdělující rovinu na dvě části. Protože Voroného buňky jsou konvexní, V(Pi) by se skládala z části ohraničené dvěma rovnoběžnými přímkami. Prokázali jsme, že hrany Voroného diagramu nemohou být přímky a to je v rozporu. B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
Obr.4 Nyní tedy rozumíme struktuře Voroného diagramu a prozkoumáme jeho složitost, která je dána celkovým počtem jeho vrcholů a hran. Máme–li tedy n bodů a každá Voroného buňka má nejvíce n-1 vrcholů a hran, složitost Vor(V) je nanejvýš kvadratická. Avšak to, že Vor(P) má právě kvadratickou složitost není zcela zřejmé. Je velmi snadné zkonstruovat příklad kde samostatná Voroného buňka má lineární složitost. Následující věta ukazuje, že se VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 19
nejedná o případ kvadratické složitosti a že průměrný počet vrcholů Voroného buněk je větší než šest. Věta 2 Pro počet vrcholů n ≥3 ve Voroného diagramu množiny n bodů v rovině je nejvýše 2n5 a počet hran je nejvýše 3n-6.
Důkaz. Jestliže body jsou všechny kolineární potom věta okamžitě vyplívá z Věty 1, takže předpokládejme, že se nejedná o tento případ. Dokažme větu použitím Eulerovy rovnice jenž nám říká, že pro jakýkoliv spojitý rovinný graf s mv uzly me oblouky a s mf plochami vyhovuje následující rovnice (3). B
mv - me + mf = 2 B
B
B
B
B
B
B
B
B
B
B
(3)
Nemůžeme aplikovat Eulerovu rovnici přímo na Vor(V), protože Vor(V) obsahuje polopřímky a není proto tedy odpovídajícím rovinným grafem. K nápravě této situace přidáme jeden extra „nekonečný“ bod do množiny bodů a spojíme všechny polopřímky Vor(V) k tomuto bodu tak jak je vidět na obr. 5. Nyní máme spojitý rovinný graf na který můžeme aplikovat Eulerovu rovnici. Dostaneme následující vztah (4) mezi počtem vrcholů Vor(V) nv , počtem hran Vor(V) ne a počtem bodů n: B
B
B
B
(nv + 1 ) - ne + n = 2 B
B
B
B
(4)
Obr.5 Přidání „nekonečného“ bodu do Voroného diagramu. Navíc každá hrana v rozšířeném grafu má přesně dva vrcholy, takže pokud sečteme stupeň všech vrcholů dostáváme dvakrát počet hran. Protože každý vrchol včetně v∞ má stupeň nejméně tři dostaneme vztah (5). B
2ne ≥ 3(nv + 1) B
B
B
B
B
(5)
Dohromady s rovnicí (4) to dokazuje tuto větu. Uzavřeme tuto kapitolu charakteristikou hran a vrcholů Voroného diagramu. Víme, že hrany jsou části půlících os dvojice bodů a že vrcholy jsou průsečíky bodů mezi těmito půlicími osami. Je zde tedy kvadratický počet půlících os, kdežto složitost Vor(V) je pouze
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 20
lineární. Proto tedy ne všechny půlící osy dané hranami Vor(V) a ne všechny průsečíky jsou vrcholy Vor(V). Popišme které půlící osy a průsečíky definují rysy Voroného diagramu – udělejme následující definici. Pro bod Q definujeme největší prázdnou kružnici bodu Q s ohledem na V označenou jako Cv(Q) jenž je největší kružnicí se středem v Q takovou, že neobsahuje uvnitř sebe jakýkoliv jiný bod množiny V . Následující věta popisuje vrcholy a hrany Voroného diagramu. B
B
Obr. 6 Kružnice pro bod Q
Obr. 7 Body ležící na kružnici
Obr.8 Voroného diagramy Věta 3 Pro Voroného diagram Vor(P) množiny bodů V platí následující:
1. Bod Q je vrchol Vor(V) tehdy a pouze tehdy když jeho největší prázdná kružnice Cp(Q) obsahuje tři nebo více bodů na svém obvodu.
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 21
2. Půlící osa mezi body Pi a Pj definuje hranu diagramu Vor(V) tehdy a pouze tehdy jestliže je zde bod Q na půlící ose takový, že Cp(Q) obsahuje oba body Pi a Pj na své hranici a neobsahuje žádný jiný bod. B
B
B
B
B
B
B
B
Důkaz (1). Předpokládejme, že existuje bod Q takový, že Cp(Q) obsahuje tři nebo více bodů na svém obvodu. Nechť Pi, Pj a Pk budou tři tyto body. Jestliže vnitřek kružnice Cp(Q) je prázdný Q musí být na hranici každé buňky V(Pi), V(Pj) a V(Pk) a Q musí být vrchol Vor(P). Na druhou stranu, každý vrchol Q diagramu Vor(V) je incidentní s nejméně třemi hranami a proto tedy s nejméně třemi Voroného buňkami V(Pi), V(Pj) a V(Pk). Bod Q musí být stejně vzdálen k Pi, Pj a Pk a nemůže zde být jiný bod poblíž bodu Q, protože jinak se V(Pi), V(Pj) a V(Pk) nemohou protnout v Q. Proto vnitřek kružnice s body Pi, Pj a Pk na jejím obvodu neobsahuje jiný bod. B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
Důkaz (2). Předpokládejme, že máme bod Q s vlastností postavené na předchozí větě. Protože Cp(Q) neobsahuje jiný bod ve svém vnitřku a Pi a Pj jsou na jejím obvodu máme dist(Q,Pi) = dist(Q,Pj) ≤ dist(Q,Pk) pro všechna 1≤ k ≤ n. Z toho vyplívá, že Q leží na hraně nebo bodu diagramu Vor(V). První část věty značí, že Q nemůže být bodem diagramu Vor(V). Proto tedy, Q leží na hraně Vor(V), jenž je definována půlící osou Pi a Pj. Naopak nechť půlící osa bodů Pi a Pj definuje Voroného hranu. Větší prázdná kružnice opsaná jakémukoliv bodu Q jenž náleží hraně a je uvnitř kružnice musí obsahovat Pi a Pj na jejím obvodu a žádný jiný bod. B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
3.3
B
B
B
B
Delaunayho triangulace
Nechť M = (P1, P2, …, Pn) je množina bodů v rovině. K tomu abychom byly schopni formálně definovat triangulaci množiny bodů M, musíme nejprve definovat maximální rovinný útvar S takový, že hrana spojující dva vrcholy nemůže být přidána do S pokud by porušila jeho rovinnost. Jinými slovy jakákoliv hrana s body množiny M, která protíná jednu z existujících hran z S není v S. Triangulace množiny M je nyní definována jako maximální rovinný útvar jehož množina bodů je M. S touto definicí je zřejmé, že triangulace existuje. Každá plocha kromě jedné neohraničené musí být trojúhelník. Ohraničená plocha je vždy polygon a jakýkoliv polygon může být triangulován [2]. A také jakákoliv hrana spojující dva sousledné body na hranici konvexního obalu (tvořící neohraničenou plochu) M je hrana v jakékoliv triangulaci T [4]. To naznačuje, že sjednocením ohraničených ploch T je vždy konvexní obal množiny M, a že neohraničená plocha je vždy součástí konvexního obalu. Počet trojúhelníků je stejný jako v jakékoliv jiné triangulaci M. To také platí pro počet hran. Přesný počet závisí na počtu bodů v M, které jsou na hranici konvexního obalu M. B
B
B
B
B
B
Věta 4 Nechť M je množina n bodů v rovině, každý z nich je nekolineární a nechť k značí počet bodů v M které leží na hranici konvexního obalu M. Potom libovolná triangulace má 2n - 2 - k trojúhelníků a 3n - 3 - k hran.
Důkaz. Nechť T je triangulace množiny bodů M a nechť m značí počet trojúhelníků T. Poznamenejme, že počet ploch triangulace, které označíme nf je m+1. Každý trojúhelník má tři hrany, a neohraničená plocha má k hran. Mimoto, každá hrana je incidentní přesně se (3m + k ) . Eulerova rovnice pro rovinný dvěma vrcholy. Proto, celkový počet hran T je ne = 2 graf nám říká, že n - ne - nf = 2 . Dosaďme do Eulerovy rovnice hodnoty pro ne a nf a dostaneme m = 2n - 2 - k a po dosazení do ne dostaneme ne= 3n - 3 - k . B
B
B
B
B
B
B
B
B
B
B
VUT v Brně, FSI ústav automatizace a informatiky
B
B
B
Delaunayho triangulace a její aplikace
Strana 22
Nechť T je triangulace množiny bodů M a předpokládejme, že má m trojúhelníků. Vezměme v úvahu 3m úhlů trojúhelníků T vzestupně seřazených. Nechť α1, α2, … ,α3m je výsledná posloupnost úhlů potom αi ≤ αj pro i ≤ j. Nazvěme A(T) = (α1, α2, …, α3m) úhlovým vektorem triangulace T. Nechť T’ je libovolná jiná triangulace stejné množiny bodů M a nechť je její úhlový vektor A(T’) = (α‘1, α‘2, …, α‘3m). Řekneme, že úhlový vektor T je větší než úhlový vektor T’ jestliže A(T) je lexikograficky větší než A(T‘) nebo jinými slovy, jestliže existuje index i kde 1 ≤ i ≤ 3m takový, že αj = α’j pro všechna j < i, a αi > α’i. Označme tuto nerovnost jako A(T) > A(T‘). Triangulace T je nazývána úhlově optimální jestliže A(T) ≥ A(T‘) pro všechny triangulace T‘ množiny bodů M. Dále v textu prostudujeme kdy je triangulace úhlově optimální. K tomu abychom to mohly udělat je užitečné znát následující větu, často nazývanou jako Thaletova věta. Označujme úhel definovaný libovolnými třemi body P, Q, R jako ∠PQR . B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
Věta 5 Nechť C je kružnice, l přímka protínající C v bodech A a B a P, Q, R a S jsou body ležící na stejné polorovině dané přímkou l tak jak je vidět na obr. 9. Předpokládejme, že P a Q leží na C, že R leží uvnitř C a že S leží vně C. Potom platí rovnice (6).
∠ARB > ∠APB = ∠AQB > ∠ASB
(6)
Obr. 9 Velikosti úhlů daných úsečkou AB a kružnicí C
Nyní uvažujme hranu e = Pi Pj triangulace T množiny bodů M. Jestliže e není hrana konvexního obalu pak je společná dvěma trojúhelníkům PiPjPk a PiPjPl. Jestliže tvar těchto dvou trojúhelníků je konvexní čtyřúhelník, můžeme získat novou triangulaci T’ odstraněním hrany Pi Pj z T a vložením hrany Pk Pl namísto ní. Tuto operaci nazýváme „edge flip“ neboli „prohození hrany“ nebo také výměna nelegální hrany za legální. B
B
B
B
B
B
B
B
B
B
B
B
Obr.10 Převracení hran „edge flip“
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 23
Rozdíl v úhlovém vektoru triangulace T a T’ je v šesti úhlech α1, …, α6 v A(T), které jsou zaměněny za α‘1, …. , α‘6 v úhlovém vektoru A(T’) tak jak je to vidět na obr. 10. Nazveme hranu e ilegální hranou jestliže platí min α i < min α 'i . B
B
B
B
B
B
B
B
1≤ i ≤ 6
1≤ i ≤ 6
Jinými slovy, hrana je nelegální jestliže můžeme lokálně zvětšit nejmenší úhel převrácením této hrany. Následující věta nutně vyplývá z definice nelegální hrany. Věta 6 Nechť T je triangulace s ilegální hranou e. Nechť T’ je triangulace odvozená z T převrácením hrany e. Potom A(T’) > A(T).
Ukážeme, že není nezbytné počítat úhly α1, …, α6 a α‘1, …, α‘6 ke kontrole, zda daná hrana je legální. Místo toho, můžeme použít jednoduché kritérium založené na následující větě. Správnost tohoto kritéria vyplývá z Thaletovy věty. B
B
B
B
B
B
B
B
B
B
Věta 7 Nechť hrana Pi Pj je společná trojúhelníkům PiPjPk a PiPjPl a nechť C je kružnice B
B
B
B
B
B
B
B
B
B
B
B
procházející body Pi Pj a Pk. Hrana je nelegální tehdy a pouze tehdy jestliže bod Pl leží uvnitř kružnice C. Mimoto také pokud body Pi, Pj, Pk Pl tvoří konvexní čtyřúhelník a neleží na společné kružnici, pak přesně jedna z hran Pi Pj nebo Pk Pl je nelegální hranou. B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
Pozorujme, že kritérium je symetrické na Pk a Pl: Pl leží uvnitř kružnice opsané bodům PiPjPk tehdy a pouze tehdy jestliže Pk leží uvnitř kružnice opsané PiPjPl. Pokud všechny čtyři body leží na kružnici, obě hrany Pi Pj a Pk Pl jsou legální. Poznamenejme, že dva trojúhelníky B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
které přísluší nelegální hraně musí vytvářet konvexní čtyřúhelník, takže je vždy možno prohodit nelegální hranu. Věta 8 Nechť V je množina bodů v rovině, které jsou v obecné pozici. Delaunayho triangulací (DT) nazveme množinu T trojúhelníků takovou kde každému trojúhelníku kružnice opsaná neobsahuje žádný jiný bod množiny V ve svém vnitřku a proto tedy má maximální úhlový vektor. Je tedy úhlově optimální triangulací a proto všechny hrany jsou legální.
Obr.11 Delaunayho triangulace množiny bodů v rovině
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 25
4
Geometrické algoritmy
4.1
Algoritmy pro konstrukci konvexní obálky
Poté co jsme si definovaly geometrický útvar zvaný konvexní obálka uvedeme také některé algoritmy pro její konstrukci. Ve spojení s algoritmy pro konstrukci DT se však nepoužívají jelikož u algoritmů pro DT je možno vždy nalézt hrany triangulace náležící konvexní obálce. 4.1.1 SlowConvexHull
Nejprve uvedeme popis tohoto algoritmu a potom provedeme rozbor tohoto algoritmu. Popis tohoto algoritmu je následující: Algoritmus SlowConvexHull(V)
Vstup: Množina bodů V Výstup: Seznam L bodů konvexní obálky seřazené v acw sledu.
1. Inicializace E jako prázdné množiny. 2. for všechny uspořádané dvojice ( P, Q) ∈ V × V , P ≠ Q 3. do platnost ← true. 4. for ∀R ∈ V , ( R ≠ P ∧ R ≠ Q) . 5. do if R leží nalevo od orientované hrany PQ 6. then platnost ← false. 7. if platnost then přidej orientovanou hranu PQ do E. 8. Z množiny hran E sestav list L vrcholů CH(V), setříděný v acw sledu.
Test na řádce 5 zda bod R leží nalevo od orientované hrany PQ se dá jednoduše provést výpočtem uvedeným v kapitole 6.2. Analýza časové složitosti tohoto algoritmu je snadná. Kontrolujeme n2-n párů bodů. Pro každý pár hledáme n-2 jiných bodů abychom věděli zda leží na pravé straně. To nám dává celkovou složitost O(n3) . Konečný krok má časovou složitost O(n2) . Takže celková časová složitost bude O(n3). Algoritmus s kubickou časovou složitostí je příliš pomalý proto se nehodí pro praktické využití. P
P
P
P
P
P
P
P
4.1.2 GHull
Tento algoritmus používá jednoduchého geometrického výpočtu pro provedení výpočtu konvexní obálky. Nevýhodou tohoto algoritmu je náchylnost na zaokrouhlovací chyby a také to, že se dá použít jen v 2D případě. Tento algoritmus nejprve najde první bod konvexní obálky tzv. pivot. To může být bod s nejmenší y souřadnicí nebo také bod s nejmenší x souřadnicí. V případě shody můžeme použít obou kriterií zároveň. Z uvedeného zápisu algoritmu je zřejmé, že jeho časová složitost je O(n2). P
P
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 26
Algoritmus GHull(V) Vstup: Množina bodů V. Výstup: Seznam L bodů konvexní obálky. 1. Nalezni bod P, který má nejmenší y souřadnici. 2. Vytvoř pomocný bod C napravo od P a B = P. 3. while A ≠ P do 4. Nalezni bod A takový který spolu s bodem B tvoří úsečku jenž svírá s úsečkou BC nejmenší úhel. 5. Ulož B do L, B = A, C = B.
4.2
Algoritmy pro konstrukci Delaunayho triangulace
V této kapitole popíšeme algoritmy pro konstrukci Delaunayho triangulace. První tři z těchto algoritmů nepoužívají metodu Divide & Concquer (Rozděl a Panuj) kde řešení spočívá na rozdělení problému triangulace na menší podproblémy. Tuto metodu zastupují dva algoritmy Guibas–Stolfi a DelaunayWall které mají také nejmenší časovou složitost. 4.2.1 Lokální zlepšování
Tento algoritmus pracuje následujícím způsobem. Mějme dva trojúhelníky které mají společnou právě jednu hranu. Pokud kružnice opsaná jednomu z trojúhelníků obsahuje vrchol druhého trojúhelníka potom hrana společná těmto dvěma trojúhelníkům je neplatná a musí být prohozena tak jak je to znázorněno na obr. 10 v kapitole 3.3. Je dokázáno [4], že toto prohození je vždy možné provést a pro celou množinu trojúhelníků je toto prohazování konečné. To znamená, že algoritmus musí najít konečnou triangulaci, která bude Delaunayho po n krocích a nezacyklí se. Problémem zůstává jak zkonstruovat počáteční libovolnou triangulaci, která se na Delaunayho teprve pomocí tohoto algoritmu převede. Ani časová složitost příliš nenabádá k použití tohoto algoritmu. Jeho časová složitost je pro nejnepříznivější případ O(n2). P
P
4.2.2 Inkrementální konstrukce
Algoritmus inkrementální konstrukce (InDeCo) je velmi jednoduchý a je založen na vlastnosti Delaunayho triangulace nejmenší kružnice opsané. Tento algoritmus pracuje následujícím způsobem. Na počátku zvolíme libovolný bod množiny určené k triangulaci. K tomuto bodu nalezneme nejbližší bod (tj. bod, který má nejmenší Eukleidovskou vzdálenost). Sestrojíme první hranu triangulace z těchto dvou nalezených bodů a pokusíme se vyhledat další bod napravo od této hrany takový, který s danou hranou tvoří trojúhelník, který má nejmenší Delaunayho vzdálenost. Ta je definována následujícím způsobem. Pokud střed kružnice opsané hraně e a bodu P leží ve stejné polorovině jako bod P potom je Delaunayho vzdálenost rovna dd(e,P)=r jinak dd(e,P)= -r kde r je poloměr kružnice opsané. Poloroviny jsou dány právě hranou e tak jak je vidět na obr.12. Pokud bod P nebyl nalezen otočíme orientaci hrany a vyhledání bodu se opakuje, tentokrát nalevo od hrany e. Po nalezení bodu P máme vytvořen první trojúhelník tvořený třemi hranami. Je důležité udržovat tyto hrany v jednom směru orientace. Tím jsme dokončili inicializaci algoritmu. Nyní tyto
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 27
hrany přidáme do AEL (Active Edge List – Seznam Aktivních Hran) vezmeme první hranu z tohoto seznamu a vyhledáme pro ni bod P s nejmenší Delaunayho vzdáleností. Pokud je nalezen vytvoříme dvě nové hrany a přidáme je do AEL ovšem pouze za předpokladu, že přidávaná hrana není již v AEL obsažena. V tomto případě je z AEL odstraněna a přidána do výsledné triangulace. Aktuální hrana pro kterou byl bod hledán je odstraněna z AEL a přidána do výsledné triangulace. Pokud pro aktuální hranu nebyl nalezen žádný bod znamená to, že hrana je součástí konvexního obalu a takovou jednoduše přidáme do výsledné triangulace.
Obr.12 Vlevo případ výpočtu kdy S neleží ve stejné polorovině jako bod P, vpravo případ kdy bod S leží ve stejné polorovině jako P. Celý postup končí ve chvíli kdy je AEL prázdný. Díky udržování orientace hran v jednom směru na hranici aktuální triangulace uložené v AEL testujeme vždy hrany spolu s body které jsou napravo (resp. nalevo) od orientovaných hran. Zjednodušeně tento postup vystihuje obr.13 na kterém jsou hrany hranice aktuální triangulace znázorněny jako orientované a hrany výsledné triangulace jako neorientované.
Obr.13 Principiální znázornění konstrukce DT pomocí inkrementálního algoritmu VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 28
Algoritmus inkrementální konstrukce
Vstup: Množina bodů V. Výstup: Množina hran E Delaunayho triangulace. 1. Inicializuj AEL jako prázdný seznam, vyber náhodný bod A ∈ V . 2. Najdi nejbližší bod B k bodu A, sestroj první hranu AB . 3. Nalezni bod C pro hranu AB s nejmenší Delaunayho vzdáleností C = min(dd( AB )). 4. if not C then otoč orientaci hrany AB a opakuj krok 3. 5. if C, A, B nejsou v acw sledu then otoč orientaci hrany AB . 6. Vytvoř hrany BC , CA PřidejDoAEL( AB ), PřidejDoAEL( BC ), PřidejDoAEL( CA ). 7. while not AEL prázdný do 8. Vezmi první hranu AB v AEL, otoč její orientaci a nalezni bod C=min(dd( AB )). 9. if not C then odstraň AB z AEL, přidej AB do E, pokračuj bodem 7. 10. Vytvoř hrany BC , CA PřidejDoAEL( BC ), PřidejDoAEL( CA ). V tomto algoritmu se vyskytuje funkce PřidejDoAEL, která nejprve testuje zda AEL přidávanou hranu již neobsahuje. Pokud ano pak tuto hranu odstraní a přidá do výsledné triangulace, pokud ne je tato hrana do AEL přidána. Funkce PřidejDoAEL ( AB ) 1. if AEL obsahuje AB then odstraň AB z AEL, přidej AB do E, return. 2. Přidej AB na konec seznamu AEL. 3. return. V algoritmu je jak již bylo řečeno udržovat hranici aktuální triangulace v jednom směru. To proto abychom na základě acw testu jednoznačně určili polohu nově hledaného bodu. Jsou–li hrany udržovány ve směru proti hodinovým ručičkám budeme vyhledávat body, které nejsou s těmito hranami v acw sledu nebo hranu otočíme a vice versa. Tento test bude uveden v kapitole implementace. 4.2.3 Inkrementální vkládání
Algoritmus je popsán v [2] následujícím způsobem. Na počátku zvolíme velký trojúhelník, který obsahuje celou množinu bodů V určenou k triangulaci. Body tohoto trojúhelníka jsou označeny P-1, P-2, P-3. To znamená, že musíme vypočíst Delaunayho triangulaci V ∪ Ω kde Ω = ( P−1 , P− 2 , P−3 ) . Po samotném výpočtu a vyřazením bodů P-1 P-2 a P3 společně se všemi incidentními hranami dostaneme Delaunayho triangulaci množiny V. Pro funkčnost musíme zvolit body P-1 P-2 a P-3 dostatečně daleko tak aby nedocházelo ke kolizím s trojúhelníky množiny V. Dále musíme zajistit aby neleželi na jakékoliv kružnici opsané třemi body v V. Popis samotného algoritmu je následující. B
B
B
B
B
B
B
B
B
B
B
B
B
B
VUT v Brně, FSI ústav automatizace a informatiky
B
B
B
B
Delaunayho triangulace a její aplikace
Strana 29
Algoritmus je náhodně vkládající, takže přidává body v náhodném pořadí a udržuje Delaunayho triangulaci množiny bodů aktuální. Vezměme v úvahu přidání bodu Pr . Nejprve najdeme trojúhelník aktuální triangulace který obsahuje Pr uvnitř a přidáme hrany z Pr do vrcholů tohoto trojúhelníka. Jestliže se stane, že Pr padne na hranu e triangulace, přidáme hrany z Pr do protějších vrcholů trojúhelníků které mají právě společnou hranu e . Obr. 14 zobrazuje tyto dva případy vkládání bodu. Po této operaci obdržíme opět triangulaci, která ale nemusí nutně být Delaunayho triangulací. To proto, že po přidání bodu Pr se mohou stát některé již existující hrany nelegálními. B
B
B
B
B
B
B
B
B
B
B
B
Obr.14 Situace které mohou nastat při vkládání nového bodu do existující triangulace Nápravu provedeme voláním procedury LegalizeEdge s každou potenciální nelegální hranou. Tato procedura zlegalizuje nelegální hrany pomocí prohození nelegálních hran. Předtím něž přejdeme k detailům, uvedeme přesný popis hlavního algoritmu inkrementálního vkládání, který je následující:
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 30
Algoritmus Inkrementálního vkládání.
Vstup: Množina bodů v rovině V. Výstup: Delaunayho triangulace množiny V. 1. Nechť P-1 P-2 a P-3 jsou vhodné tři body takové, že množina bodů V je obsažena v trojúhelníku P-1 P-2 P-3 B
B
B
B
B
B
B
B
B
B
B
B
2. Inicializuj T jako triangulaci obsahující jediný trojúhelník P-1 P-2 P-3 . B
B
B
B
B
B
3. Vypočti náhodnou permutaci P1 , P2 , …, Pn množiny V. B
B
B
B
B
B
4. for r = 1 to n 5. do (Vkládání bodu Pr do T ) B
B
6. Najdi trojúhelník obsahující Pr. B
B
7. if pr leží uvnitř trojúhelníka pi pj pk B
B
B
B
B
B
B
B
8. then Přidej hrany z Pr do vrcholů Δ Pi Pj Pk, tak, že se Δ Pi Pj Pk rozdělí na tři nové trojúhelníky. B
B
B
B
B
B
B
B
B
B
B
B
B
B
9. LegalizeEdge(Pr, Pi Pj, T) B
B
B
B
B
B
10. LegalizeEdge(Pr, Pj Pk, T) B
B
B
B
B
B
11. LegalizeEdge(Pr, Pk Pi, T) B
B
B
B
B
B
12. else (Pr leží na hraně. Přidáme hrany tak, že vzniknou čtyři trojúhelníky.) B
B
13. LegalizeEdge(Pr, Pi Pl, T) B
B
B
B
B
B
14. LegalizeEdge(Pr, Pl Pj, T) B
B
B
B
B
B
15. LegalizeEdge(Pr, Pj Pk, T) B
B
B
B
B
B
16. LegalizeEdge(Pr, Pk Pi, T) B
B
B
B
B
B
17. Odstraň body P-1 P-2 P-3 a všechny hrany incidentní s těmito body. B
B
B
B
B
B
18. return T Dále probereme detaily změny triangulace kterou dostaneme po řádce 8 nebo řádce 13 v tomto algoritmu. Triangulace je Delaunayho, jestliže všechny její hrany jsou legální. Proto tedy musíme převracet nelegální hrany dokud není triangulace opět legální. Otázka která zůstává je, které hrany se mohly stát nelegálními při vkládání bodu Pr. Pozorujme, že hrana která byla předtím legální se může stát nelegální pouze tehdy jestliže jeden z trojúhelníků s ní incidentní byl změněn. Takže zkontrolovány musí být pouze hrany nových trojúhelníků. Tuto kontrolu zajistíme použitím procedury LegalizeEdge, která testuje a převrací případné ilegální hrany. Jestliže ale procedura LegalizeEdge převrátí hranu, jiné hrany se tak mohou stát opět nelegálními. Proto tedy procedura LegalizeEdge volá sama sebe rekurzívně a jako parametr má vždy takové další potenciální nelegální hrany. B
B
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 31
Procedura LegalizeEdge(Pr, Pi Pj , T) B
B
1. Vkládaný bod je Pr a Pi Pj je hrana z T která má být prohozena. B
B
2. if Pi Pj je nelegální 3. then Nechť Pi Pj Pk je trojúhelník přilehlý k Pr Pi Pj přes Pi Pj . B
B
B
B
B
B
B
B
B
B
B
B
4. Zaměň Pi Pj za Pr Pk 5. LegalizeEdge(Pr, Pi Pk , T) B
B
6. LegalizeEdge(Pr, Pk Pj , T) B
B
Test na řádce 2 zda hrana je nelegální by mohl normálně být proveden pomocí testu kružnice opsané uvedeném ve větě (5). Nicméně máme zde nějaké komplikace plynoucí z přítomnosti speciálních bodů P-1 P-2 a P-3. Tento problém probereme dále v textu. Prvně dokažme, že tento algoritmus je správný. K zaručení správnosti algoritmu, potřebujeme prokázat, že nelegální hrany budou po všech voláních procedury LegalizeEdge zpracovány. Z kódu LegalizeEdge je zřejmé, že každá nová hrana vytvořena kvůli vložení Pr je incidentní s Pr . Naznačuje to obr. 15. Trojúhelníky které byly odstraněny a nové trojúhelníky jsou vyplněny šedě. Rozhodující pozorování (provedené níže) je, že každá nová hrana musí být legální, takže není potřeba je testovat. Spolu s dřívějším pozorováním, že každá hrana se může stát legální pouze když jeden z jejich incidentních trojúhelníků se změní to prokazuje, že algoritmus testuje jakékoliv hrany které se mohly stát nelegálními. Odtud plyne, že algoritmus je správný. Algoritmus se nemůže dostat do nekonečné smyčky, protože každé převrácení zvětšuje úhlový vektor triangulace. B
B
B
B
B
B
B
B
B
B
Obr.15 Změna triangulace po přidání nového bodu Pr B
Zatímco budujeme Delaunayho triangulaci, tak také budujeme strukturu která je orientovaným acyklickým grafem a nazýváme ji DAG (Data Acyclic Graph). Listy DAG odpovídají trojúhelníkům aktuální triangulace T. Vnitřní uzly DAG odpovídají trojúhelníkům které byly v triangulaci v předchozích stupních ale již byly zničeny. Struktura je postavena následujícím způsobem. Na řádce 2 inicializujeme DAG jako DAG s jediným listem, který odpovídá velkému trojúhelníku P-1 P-2 a P-3 . Nyní dejme tomu, že v nějakém bodě rozdělíme trojúhelník aktuální triangulace na tři nebo dva nové trojúhelníky. Odpovídající změna v DAG je přidání třech nebo dvou nových B
B
B
B
B
B
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Obr.16 Znázornění DAG struktury a odpovídajících změn v ní.
VUT v Brně, FSI ústav automatizace a informatiky
Strana 32
Delaunayho triangulace a její aplikace
Strana 33
listů do DAG a vytvoření listu pro Pi, Pj, Pk v interním uzlu a vytvořením ukazatelů na tyto tři nebo dva listy. Podobně, když vyměníme dva trojúhelníky Pk Pi Pj a Pi Pj Pl trojúhelníky Pk Pi Pl a Pk Pl Pj převrácením hrany, vytvoříme listy pro dva nové trojúhelníky a uzly trojúhelníků Pk Pi Pj a Pi Pj Pl nastaví ukazatele na tyto dva nové listy. Obr. 15 znázorňuje příklad změny v DAG struktuře způsobené přidáním bodu. Všimněme si, že když vytvoříme list ve vnitřním uzlu, vytvoří se nejvíce tři ukazatele. Používání DAG struktury můžeme lokalizovat další libovolný bod Pr, který je přidán do aktuální triangulace. Začneme v kořeni DAG stromu, který koresponduje s počátečním velkým trojúhelníkem P-1 P-2 a P-3. Zkontrolujeme všechny potomky kořene abychom viděli v kterém z nich Pr leží a sestoupíme na odpovídajícího potomka. Potom zkontrolujeme potomka tohoto uzlu, sestoupíme na potomka jehož trojúhelník obsahuje Pr a tak dál, dokud nedosáhneme odpovídajícího listu z DAG. Tento list odpovídá trojúhelníku v aktuální triangulaci který obsahuje Pr. Výstupní stupeň jakéhokoliv uzlu je nejvíce tři, vyhledání odpovídajícího trojúhelníka zabere lineární čas vzhledem k počtu uzlů, respektive na cestě po uzlech které obsahují Pr. Problémem tedy je jak zvolit P-1 P-2 a P-3 a jak implementovat test jestli je hrana legální či nikoliv. Někdo tvrdí, že máme zvolit P-1 P-2 a P-3 velmi vzdáleně, protože nechceme jejich přítomnost kvůli ovlivnění DT(V). Na druhou stranu, nechceme uvádět ani obrovské souřadnice kvůli numerickým omezením. Proto tedy vybereme tyto body tak, že trojúhelník P-1 P-2 a P-3 obsáhne celou množinu V a modifikujme test pro ilegální hrany tak, že bude pracovat tak jako bychom měli tyto body vzdáleny velmi daleko. Přesněji, položme P1 = (3M, 0), položme P2 = (0, 3M) a položme P3 = (-3M, -3M), kde M je absolutní hodnota maximální souřadnice bodu z V. To nám zaručí, že celá množina V bude obsažena v trojúhelníku P-1 P-2 P-3. B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
Obr.17 Volba souřadnic základního trojúhelníku. Nicméně, když kontrolujeme jestli je hrana nelegální či nikoliv nepoužijeme tyto souřadnice. Místo toho předstíráme, že P-1 leží vně všech kružnic definovaných jakýmikoliv třemi body v P, že P-2 leží vně všech kružnic definovaných libovolnými třemi body v V ∪ {P−1 } a že P-3 leží vně všech kružnic definovaných jakýmikoliv třemi body v B
B
B
B
B
B
V ∪ {P−1 , P−2 }. Proto implementujeme část funkce která testuje hrany následovně. Nechť Pi Pj
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 34
je hrana určená k testování. Potom dostáváme pět případů které mohou nastat při testování legálnosti této hrany. Případ 1: Oba indexy i a j jsou negativní (tedy jsou to body z (P-1, P-2, P-3). B
B
B
B
B
B
V tomto případě rozhodneme že je legální, protože hrany velkého trojúhelníka musíme udržet stejné jako v inicializaci. Pro jiné případy, nechť Pk a Pl jsou dva zbylé vrcholy trojúhelníků incidentních s Pi Pj . B
B
B
B
Případ 2: Indexy i, j, k, l jsou všechny kladné.
Toto je normální případ. Žádný z bodů zahrnutý v testu není speciálním bodem. Proto, hrana je ilegální tehdy a pouze tehdy jestliže Pl leží uvnitř kružnice opsané Pi, Pj a Pk . B
B
B
B
B
B
B
B
Případ 3: Přesně jeden z indexů i, j, k, l je negativní.
Nechceme aby speciální bod způsobil zničení jakýchkoliv Delaunayho hran mezi body množiny V. Proto, jestliže i nebo j je negativní (tedy Pi nebo Pj je speciální bod) pak rozhodneme že Pi Pj je ilegální tudíž bude vyměněna hranou Pk Pl . V opačném případě B
B
B
B
rozhodneme že je legální. Případ 4: Přesně dva z indexů i, j, k, l jsou negativní.
Nyní jeden z i, j a jeden z k, l musí být negativní. Nemůže nastat případ kdy oba k, l jsou negativní a případ, že oba indexy i, j jsou negativní to je případ 1 uvedený výše. Jestliže negativní index z i, j je menší než negativní index z k, l potom rozhodneme že Pi Pj je legální, jinak rozhodneme že je nelegální. Případ 5: Přesně tři z incidentních i, j, k, l jsou negativní.
Tento případ se nemůže vyskytnout: i a j nemohou být oba negativní to byl případ 1. Indexy k a l nemohou být oba negativní protože jeden z Pk, Pl, musí být bod Pr který je právě vkládán. B
B
B
B
B
B
4.2.4 Guibas–Stolfi
Tento algoritmus odpovídá standardnímu „rozděl a panuj“ paradigmatu pro řešení problému jeho rekurzivním rozdělením na menší podproblémy a následným sloučením výsledků podproblémů. Vlastnímu rekurzivnímu algoritmu předchází předzpracování dat, kdy jsou body vstupní množiny setříděny vzestupně (nejdříve podle souřadnice x a pokud dojde ke shodě, potom podle souřadnice y). Seřazení bodů nám umožní provádět rozdělení množiny bodů na dvě poloviny přímkou kolmou na osu x v konstantním čase. Seřazenou množinu bodů pak rekurzivně dělíme na dvě poloviny levou a pravou označme L a R dokud nedosáhneme množin které obsahují dva nebo tři body kdy je triangulace triviální (hrana nebo trojúhelník). Následovně spojujeme tato rozdělení do jedné triangulace tak jak se vracíme z nižší úrovně rekurze do vyšší. Princip rozdělení a výpočtu Delaunayho triangulace touto metodou pro 8 bodů v rovině naznačuje obr. 18.
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 35
Obr.18 Principiální znázornění postupu tvorby DT metodou Guibas–Stolfi Triangulace je dokončena spojením dílčích triangulací L a R na nejvyšší úrovni rekurze. Slučování dílčích triangulací je nejnáročnějším krokem algoritmu. Při slučování jsou přidány do triangulace nové hrany spojující L a R (křížové hrany, cross edges) a případně odstraněny některé nevyhovující hrany. Křížové hrany jsou určeny směrem odspodu nahoru v pořadí, ve kterém by protínaly přímku kolmou na osu x a oddělující obě dílčí triangulace. První takovou hranou je hrana spojující nejspodnější body množin L a R (tato hrana je součástí konvexního obalu. Další křížové hrany nalezneme tímto způsobem: B
B
B
B
1. Najdeme „nejlepší možnou“ křížovou hranu vycházející z levého vrcholu nejhornější křížové hrany. 2. Najdeme „nejlepší možnou“ křížovou hranu vycházející z pravého vrcholu nejhornější křížové hrany. 3. Z obou nalezených hran vybereme nejlepšího kandidáta a tuto hranu vložíme do triangulace jako novou křížovou hranu (a tato hrana se stává zároveň novou nejhornější hranou). 4. Algoritmus končí přidáním křížové hrany spojující nejhornější body množin L a R.
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 36
Obr.19 Slučovací fáze algoritmu Guibas–Stolfi Zbývá uvést kritéria pro nejvhodnější kandidáty na vkládané křížové hrany. Označme si prozatímní nejhornější křížovou hranu jako AB , kde bod A ∈ L a B ∈ R . V prvním kroku hledáme křížovou hranu vycházející z bodu A, tedy bod C ∈ R takový, aby splňoval naše kritérium optimality. Postup bude následující (viz obr. 20). Zvolíme C ∈ R jako prvního souseda bodu B napravo od orientované hrany BA . Tento bod budeme považovat za prozatímního nejlepšího kandidáta na koncový bod křížové hrany. Označme souseda bodu B napravo od orientované hrany BC jako C´. Pokud leží bod C´ „pod“ hranou AB (neboli napravo od orientované hrany AB ), bod C´ je „neplatný“ a hledání končí. Pokud bod C´ leží vně kružnice opsané trojúhelníku ABC, hledání končí. Pokud bod C´ leží uvnitř kružnice opsané trojúhelníku ABC, stává se novým nejlepším kandidátem. Hrana AC je odstraněna z triangulace a bod C´ si označíme jako C. Pokračujeme znovu od bodu 2. Postup hledání nejlepší křížové hrany v druhém kroku algoritmu (vycházející z bodu B) je analogický. Tímto jsme nalezli body a určující křížové hrany AC a BD . Pokud je některý z bodů C, D neplatný, vybereme tu křížovou hranu, která ho neobsahuje. V opačném případě vybereme přidávanou křížovou hranu na základě již zmíněného kružnicového testu: pokud je kružnice opsaná trojúhelníku ABC menší než kružnice opsaná trojúhelníku ABD, přidáme hranu AC (v opačném případě přidáme hranu BD ).
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Obr.20 Postup tvorby křížových hran
VUT v Brně, FSI ústav automatizace a informatiky
Strana 37
Delaunayho triangulace a její aplikace
Strana 38
4.2.5 Delaunay Wall
Algoritmus DelaunayWall je dalším ze skupiny algoritmů „rozděl a panuj“ pro konstrukci DT. Na rozdíl od předchozího algoritmu není omezen na rovinnou triangulaci. Princip „rozděl a panuj“ je zde aplikován odlišně od předchozího algoritmu. Algoritmus Delaunay Wall (krátce DeWall) proto používá jiný přístup. Při něm se poněkud zkomplikuje fáze dělení množiny bodů, na druhou stranu odpadnou problémy při slučování. Myšlenka je následující. Množina bodů V je rozdělena dělící nadrovinou α na dvě disjunktní podmnožiny. Trojúhelníky tvořící budovanou triangulaci můžeme rozdělit do tří skupin: 1. Trojúhelníky mající neprázdný průnik s α (tvoří tzv. trojúhelníkovou bariéru) 2. Trojúhelníky ležící v levém poloprostoru určeném nadrovinou α 3. Trojúhelníky ležící v pravém poloprostoru určeném nadrovinou α Trojúhelníkovou bariéru můžeme použít jako vhodnou slučovací triangulaci, pokud každý trojúhelník bariéry je obsažen ve výsledné triangulaci a pokud odstraněním trojúhelníkové bariéry z výsledné triangulace obdržíme dva oddělené shluky trojúhelníků. Princip algoritmu je následující: 1. Zvolíme dělící nadrovinu α 2. Množinu bodů V rozdělíme dělící nadrovinou na dvě disjunktní podmnožiny V1 a V2 a sestrojíme trojúhelníkovou bariéru Σα B
B
B
3. Algoritmus aplikujeme rekurzivně na V1 a V2 , tím získáme triangulace Σ1 a Σ2 B
B
B
B
B
B
B
B
Výsledná triangulace je sjednocením Σ1, Σ2 a Σα. B
B
B
B
B
B
Obr.21 Konstrukce Delaunayho triangulace pomocí algoritmu DeWall
VUT v Brně, FSI ústav automatizace a informatiky
B
B
B
Delaunayho triangulace a její aplikace
Strana 39
Jako dělící nadrovinu obvykle volíme cyklicky nadrovinu kolmou na osy souřadného systému. Získané podmnožiny by měly mít přibližně stejnou mohutnost, což zajistíme např. použitím mediánu dělené množiny bodů. Trojúhelníkovou bariéru (množinu trojúhelníků mající neprázdný průnik s nadrovinou α sestrojíme pomocí jednoduchého inkrementálního algoritmu, který trochu připomíná algoritmus inkrementální konstrukce DT. Při implementaci se používají tři pomocné seznamy hran trojúhelníků – AEL (Active Edge List). Tyto seznamy v každém okamžiku obsahují hrany tvořící vnější hranici dosud nalezené triangulace. Jsou implementovány tak, že při vložení hrany do seznamu se otestuje, jestli tato hrana již v seznamu neexistuje. Tento systém je tedy stejný jako u již uvedeného algoritmu inkrementálního vkládání. Pokud seznam hranu obsahuje, je tato hrana ze seznamu odebrána, v opačném případě je do seznamu vložena. Algoritmus používá tři takové (disjunktní) seznamy: 1. AELα hrany mající průnik s nadrovinou α B
B
2. AEL1 hrany mající všechny vrcholy v V1 B
B
B
B
3. AEL2 hrany mající všechny vrcholy v V2 B
B
B
B
Při rekurzivním dělení se algoritmu konstrukce trojúhelníkové bariéry předává dělená množina bodů V a seznam aktivních hran AEL, který je při prvotní inicializaci algoritmu prázdný. Tento seznam se posléze rozdělí do tří podseznamů AELα, AEL1 a AEL2 podle dělící roviny. Seznam AELα představuje seznam hran hranice trojúhelníkové bariéry majících průsečík s α. Pokud je tento seznam na počátku dělení prázdný, je nalezen první trojúhelník tvořící trojúhelníkovou bariéru. Postup jeho nalezení je následující. Jako první vrchol trojúhelníku je zvolen bod, který leží nejblíže nadrovině α. Jako druhý bod P2 vybereme ten bod opačného poloprostoru (vzhledem k α), který má nejkratší vzdálenost od P1 . Další vrchol trojúhelníku vybíráme tak, aby poloměr kružnice opsané vzniklé hraně byl minimální. Nalezený trojúhelník je vložen do Σα a jeho hrany do příslušných seznamů AEL. Bariéru budeme konstruovat tak, že budeme hledat trojúhelníky sousedící s touto zatím nalezenou hranicí splňující Delaunayho kritérium nejmenší kružnice opsané. Pokud je seznam AELα neprázdný, obsahuje alespoň jednu hranu f. Tuto hranu ze seznamu vyjmeme a hledáme takový bod P, který společně s hranou f vytvoří hledaný minimální trojúhelník. Bude to bod, který minimalizuje tzv. Delaunayho vzdálenost. Ta je definována stejně jako u algoritmu inkrementální konstrukce. Nalezený trojúhelník přidáme do Σα a jeho hrany (kromě hrany f) vložíme v závislosti na jejich poloze vůči dělící nadrovině do seznamů AEL. Tento postup opakujeme do té doby, dokud nebude seznam AELα prázdný. Tím je konstrukce dílčí trojúhelníkové bariéry dokončena. Pokud jsou aktivní seznamy AEL1 a AEL2 neprázdné, algoritmus konstrukce trojúhelníkové bariéry potom rekurzivně aplikujeme na dvojice (P1, AFL1) a (P2, AFL2). Body P minimalizující Delaunayho vzdálenost od hrany f má smysl hledat jen ve vnějším poloprostoru vzhledem k f (tedy poloprostoru neobsahujícím trojúhelník vygenerovaný v předchozím kroku, ke kterému patřila hrana f). Tím zajistíme, že budeme pokračovat v konstrukci bariéry správným směrem. Pokud je hrana f částí konvexního obalu množiny bodů P, neobsahuje příslušný vnější poloprostor žádný bod z P. To znamená, že dále není možné tímto směrem hranici rozšiřovat a že hrana f je zároveň částí vnější hranice triangulace. B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
VUT v Brně, FSI ústav automatizace a informatiky
B
Delaunayho triangulace a její aplikace
5
Aplikace Delaunayho triangulace
5.1
Výšková interpolace
Strana 41
Mějme model části zemského povrchu, dále jen jako terén. Terén je dvourozměrný povrch v trojrozměrném prostoru se speciálními vlastnostmi. Je to graf funkce f : A ⊂ R 2 → R takové, že přiřazuje výšku f(P) každému bodu P v oblasti terénu A. (Země je kulatá, takže v globálním měřítku takto definované terény nejsou dobrým modelem povrchu země. Ale ve větším počtu v malém místním měřítku poskytují docela dobrý model). Samozřejmě neznáme výšku každého bodu na Zemi. Známe pouze výšku těch bodů, pro které máme tuto výšku změřenu. Tím se myslí to, že když mluvíme o nějakém terénu, tak známe pouze hodnoty funkcí v konečné množině N ⊂ A vzorových bodů. Z výšky vzorových bodů máme nějak aproximovat výšku ostatních bodů v této oblasti. Můžeme prostým způsobem přiřadit každému bodu P z množiny A výšku nejbližšího vzorového bodu. Avšak, tím dostáváme terén diskrétní, který nevypadá moc přirozeně. Proto tedy náš přístup k aproximaci terénu bude následující. V předchozí kapitole 3.3 jsme vymezili pojem triangulace množiny bodů P. Předpokládáme, že vzorové body jsou takové, že můžeme utvořit trojúhelníky pokrývající oblast terénu. Poté zvedneme každý vzorový bod do jeho správné výšky, tím zobrazíme (převedeme) každý trojúhelník v rovinné triangulaci na trojúhelník v 3D prostoru.
Obr.22 Získání mnohostěnného terénu z množiny vzorových bodů.
Na obr. 22 je tato procedura zobrazena. To co dostaneme je mnohostěnný terén, graf spojité funkce která je po částech lineární. Mnohostěnný terén můžeme použít jako aproximaci originálního terénu. Otázkou zůstává, jak sestrojit triangulaci množiny vzorových bodů. Obecně, to může být provedeno mnoha způsoby ale jaká triangulace je nejvhodnější pro účel aproximace terénu ? VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 42
Neexistuje jednoznačná odpověď na tuto otázku. Neznáme celý původní terén pouze známe jeho výšku ve vzorových bodech. Do té doby nemáme žádnou jinou informaci a výška v místě vzorových bodů je daná pro jakoukoliv libovolnou triangulaci. Tedy všechny triangulace množiny V se mohou zdát stejně dobré. Nicméně, některé triangulace vypadají přirozeněji něž ostatní. Na obr. 23 je příklad který ukazuje dvě triangulace stejné množiny bodů.
Obr.23 Rozdíl v triangulacích může způsobit značný rozdíl ve výšce bodu.
Obr.24 Příklad modelu terénu získaného pomocí DT Z výšek vzorových bodů dostáváme dojem že vzorové body jsme vzaly z horského hřebene. Triangulace vlevo na obr. 23 odráží toto tušení. Triangulace vpravo ačkoliv je zde jen jedna hrana přehozena, představuje úzké údolí nebo ostrý žlab či horskou řeku v tomto horském hřebeni. Problém triangulace vpravo je, že výška bodu Q je určena dvěma body které jsou relativně vzdálené. Stane se tak protože Q leží ve středu hrany dvou dlouhých a ostrých trojúhelníků. Ostrost těchto trojúhelníků způsobuje problémy, takže je vidět, že triangulace
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 43
která obsahuje velmi malé úhly v trojúhelnících není dobrá pro účely aproximace terénních dat. Proto tedy vhodnou triangulací pro tyto účely je DT která maximalizuje minimální úhel.
5.2
Eukleidovský Steinerův minimální strom
Problém Eukleidovského Steinerova minimálního stromu je v nalezení nejkratší sítě dané množinou pevných bodů v rovině, dovolující přidávání pomocných bodů do této množiny. Problém se řadí do skupiny NP–těžkých problémů. Mějme dánu množinou pevných bodů P={P1,P2,....,Pn} v Eukleidovské rovině, nazývaných terminály, a hledejme nejkratší přímou spojnici těchto bodů. Řešení má podobu stromu, nazývaného Eukleidovský Steinerův minimální strom (EStMT). Na rozdíl od problému minimální kostry, spojnice vrcholů v EStMT nejsou vyžadovány pouze mezi terminály množiny V. Přídavné neterminálové body, nazývané jako Steinerovy body nám umožní dosáhnout nejkratší sítě. Eukleidovská vzdálenost dvou bodů v rovině je dána rovnicí (1), která je uvedena v kapitole 2.4. Aplikace Eukleidovského Steinerova stromu má velké rozpětí, od modelování evoluce druhů v biologii k modelování mazacích filmů pro drátěné mřížky, od navrhování sběru dat k navrhování větracích nebo klimatizačních systémů v budovách a od vytváření návrhu rozložení olejových a benzinových potrubí k vytváření komunikačních sítí, rozvodným sítím elektrického proudu nebo silničních či železničních sítí. B
B
B
B
B
B
5.2.1 Základní pojmy Definice 1 Nechť G = (V,E) je spojitý neorientovaný graf ve kterém každé hraně e je přiřazena reálná nezáporná váha (nebo cena) w(e) a kde V je množinou vrcholů a E je množina hran grafu G.
1. Kostra grafu T je acyklický souvislý podgraf grafu G obsahující každý bod V. 2. Minimální kostra grafu (MSpT) je kostra grafu minimální váhy. Definice 2 Eukleidovský Steinerův poměr, označený jako ρE je supremum přes množinu všech poměrů délky Eukleidovské minimální kostry w(EMSpT(V)) k délce Eukleidovského Steinerova minimálního stromu w(EStMT(V)). Tento poměr vychází z EMSpT a EStMT pro tři body, které tvoří rovnostranný trojúhelník tak jako na obr. 25. Steinerův bod pro případ uspořádání třech bodů takového, jenž tvoří rovnostranný trojúhelník odpovídá středu kružnice opsané tomuto trojúhelníku. A hrany jenž vytváří EStMT pro tento případ jsou poloměry r B
B
Obr.25 a) Graf pro množinu třech bodů vytvářející rovnostranný trojúhelník. b) Eukleidovská minimální kostra. c) Eukleidovský Steinerův minimální strom. kružnice opsané. Označme strany rovnostranného trojúhelníka jako a. Potom váha
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 44
3a 3 = a 3. 3 Steinerovu poměru pro Eukleidovský problém bude tedy vyhovovat následující rovnice (7).
Eukleidovské minimální kostry bude w(EMSpT) = 2a a váha w( EStMT ) = 3r =
ρ E (V ) = sup V ⊂E
2
w ( EMspT (V )) 2a 2 = = ≈ 1 .1547 w ( EStMT (V )) a 3 3
(7)
To znamená, že délka EMSpT není větší než délka EStMT více jak o 15.47% (průměrná hodnota je samozřejmě menší). Proto tedy EMSpT může být použita jako přibližný Steinerův strom a přirozeně se stává měřítkem kterémukoliv aproximačnímu nebo heuristickému algoritmu. Problém konstrukce EMSpT může být tak jednoduše zredukována na vyhledávání váhově minimálních podgrafů ohodnoceného grafu. Pokud máme pouze množinu bodů v Eukleidovské rovině, musíme nejdříve zkonstruovat graf s jejími vrcholy které reprezentují body. Nejednoduší cesta je spojení všech dvojic vrcholů hranami. Váhy hran v tomto úplném grafu odpovídají Eukleidovským vzdálenostem mezi těmito body. Počet hran je dán rovnicí (8). Časová složitost je O(n2) kde n = |P|. P
P
n 2
=
n(n − 2) (8) 2
Problém minimální kostry grafu může být jednoduše řešen pomocí dobře známého Jarníkova (Primova), Kruskalova nebo Borůvkova polynomiálního algoritmu. Všechny tyto algoritmy mají stejný asymptotický čas výpočtu O((|V|+|E|) log |V|) pro graf G=(V,E) ale jelikož časová složitost konstrukce úplného grafu je vyšší, celkový čas algoritmu je O(|V| 2) kde V = P. Jiný přístup k řešení EMSpT je založen na koncepci Delaunayho triangulace. Odkažme se na fakt že, když hledáme hranu s minimální váhou v algoritmu pro minimální kostru grafu, vede to k vyhledávání přes hrany Delaunayho triangulace. Výpočet DT může být proveden v čase O(n log(n)) a proto tedy výpočet EMSpT z grafu DT je efektivnější něž výpočet z úplného grafu. P
P
5.2.2 Vlastnosti Steinerova minimálního stromu
Eukleidovské Steinerovy minimální stromy mají následující vlastnosti. 1. Všechny Steinerovy body jsou incidentní přesně se 3 hranami svírající úhel 120° s každou z nich.(t.j. jsou stupně 3 s ohledem na hrany použité ve stromu). Označme tuto vlastnost jako úhlová a stupňová podmínka. 2. EStMT stromy pro n terminálů mají nejvýše n-2 Steinerových bodů. 3. EStMT stromy jsou sjednocení plných Steinerových stromů. (FST). FST stromy mají o dva Steinerovy body méně než mají terminálů a terminály překlenuté stromem FST stromu jsou stupně 1. Jestliže dva FST stromy sdílí terminál z, potom dvě hrany incidentní se z (jedna z každého FST stromu) vytváří nejméně úhel 120° se všemi VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 45
ostatními. Žádný terminál proto tedy nemůže být ve více než třech FST stromech. Ke správnému určení plných Steinerových stromů dané podmnožiny terminálů, je zapotřebí uvažovat všechny jejich plné topologie. Plná topologie specifikuje jak jsou terminály a Steinerovy body spojeny. 5.2.3 Heuristický algoritmus založený na Delaunayho triangulaci
Jestliže |V| = 3 můžeme přímo zkonstruovat EStMT takto: Nechť V={A,B,C}. 1. Jestliže velikost jednoho z úhlů trojúhelníku ABC je alespoň 120°, potom EStMT se jednoduše skládá z dvou hran přilehlých k tomuto tupému úhlu. 2. Jestliže všechny vnitřní úhly trojúhelníku ABC jsou menší než 120° potom nakreslíme rovnostranný ABD a obkreslíme kružnici okolo tohoto trojúhelníku. Steinerův bod S je dán průsečíkem přímky CD s kružnicí opsanou trojúhelníku ABD tak jak je znázorněno na obr. 26. Můžeme ukázat, že celková délka segmentů AS, BS, CS je rovna délce segmentu CD, který je znám jako Simpsonova přímka pro FST strom přes terminály A, B, C.
Obr.26 Vlevo EStMT pro nerovnostranný trojúhelník, vpravo pro rovnostranný trojúhelník. Z obr. 25 můžeme snadno odvodit, že Euklidovské Steinerovy minimální stromy pro body ležící ve vrcholech rovnostranných trojúhelníků jsou 2 / 3 –krát kratší než jejich odpovídající Euklidovské minimální kostry a rovnicí (1) reprezentují nejlepší zlepšení s ohledem na minimální kostry grafu. Všechny úhly v rovnostranných trojúhelnících jsou rovny 60° a evidentně minimální úhel v takovém trojúhelníku je maximalizován. Z kapitoly 3.3 víme, že jakákoliv Delaunayho triangulace vrcholů V maximalizuje minimální úhel přes všechny trojúhelníky množiny vrcholů V a proto použijeme Delaunayho triangulaci DT(V). Jak jsme se zmínili výše, můžeme použít DT(V) pro výpočet EMSpT. Dosáhnout lepší počáteční aproximace můžeme také tak, jestliže nejprve nahradíme všechny trojúhelníky které mají každý úhel menší než 120° EStMT stromy a potom najdeme EMSpT. Délka grafu získaného touto procedurou může být zmenšena opakováním výběrů nejmenšího úhlu založeného na faktu, že nejmenší úhel je lepší pro dosažení lepšího lokálního zlepšení. Samozřejmě musí být splněno, že tento úhel je menší než 120°. Výměna takových podgrafů jejich EStMT stromy je stejná jako nahrazení uvedené předtím. Algoritmus můžeme tedy popsat následovně:
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 46
1. Najdi Delaunayho triangulaci pro danou množinu V. 2. Nahraď všechny trojúhelníky které mají každý úhel menší než 120° EStMT stromy. 3. Urči EMSpT pro graf nalezený v kroku 2. 4. For každou hranu spojující body X a Y do 5. Nalezni hranu YZ která se setkává s XY v nejmenším úhlu. 6. If úhel je menší než 120° then 8. Umísti nový Steinerův bod Sn do trojúhelníku daného hranou XY a YZ . B
B
9. Odstraň hrany XY a YZ 10. Přidej hrany XS n , YS n a ZS n . 11. Odstraň všechny Steinerovy body stupně 1 dohromady s jejich incidentními hranami.
5.3
Geometrická podobnostní metrika pro případově založené usuzování
Delaunayho triangulace byla využita v [6] . Jedná se o problém používající systém případově založeného usuzování, který slouží pro ukládání, znovu používání a návrh nejlepší možné cesty mezi dvěma body pro reálnou mapu města. Toto usuzování bylo použito pro navigaci robota ve městě které obsahuje přes 18000 sekcí 5000 ulic uložených v 25000 segmentech. 5.3.1 Reprezentace problému
Když systém pro generování plánu (dále jako SGP) generuje plán, výsledná dráha je buď uložena jako případ nebo rozdělena do několika částí a uložena jako sekvence případů. Reprezentace každého případu zahrnuje detailní popis situací odehrávajících se v prováděném čase, zahrnující vysvětlivky jakýchkoliv chyb které nastaly a všechna přeplánování která byla provedena ke korekci problémů. Pro oblast plánování cesty je každý případ také aproximován úsečkou ve dvou rozměrném grafu a úsečkám je dovoleno protínat se pouze v jejich koncových bodech. Tento graf vystupuje jako index v případové knihovně. Když SGP generuje plán tak, že protíná existující segmenty, potom plán a segmenty které protnul jsou rozděleny do menších segmentů. Výsledný graf, který nazýváme případový graf je znázorněn na obr. 27 kde vlevo je mapa a vpravo je naznačen abstraktní způsob v němž jsou cesty uloženy v indexovém souboru. Poznamenejme, že případ 20 (Case 20) zjednodušuje cestu ale skupina ulic na této cestě nemůže změnit koncovou trasu, takže toto zobecnění je akceptovatelné. Jestliže jakýkoliv z těchto segmentů je potřebný pro znovu použití, podobnostní metrika bude odhadovat nejlepší místo uvnitř každého segmentu který bude počátkem použití segmentu. Nazveme segmenty případového grafu jako případové segmenty. Jestliže případový graf obsahuje segmenty které se protínají, potom tyto segmenty musí být v těchto průsečících rozděleny a musí být přidán vrchol ke každému takovému průsečíku. Tento proces je proveden inkrementálně, právě když je každý případ přidán do případové knihovny.
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 47
Obr.27 Vlevo Mapa na níž plné úsečky jsou dříve navštívené ulice a tečkované segmenty jsou nenavštívené ulice. Vpravo reprezentace mapy případovým grafem. K reprezentaci případů je použita aproximace úseček. Více případových segmentů dohromady může popisovat jednu cestu a v jednom segmentu může být obsaženo více ulic.
5.3.2 Geometrický problém
Zabývejme se plánem cesty na naší mapě z nějakého počátečního místa i do cílového místa g. Ačkoli chceme znovu použít případů, preferujme raději neprobádanou oblast k dlouhým a klikatícím se cestám. To je důležité pro nalezení řešitelného kompromisu mezi starou cestou a hledáním nějaké nové cesty, protože přiřazujeme každému případovému segmentu efektivní hodnotu β která je hrubou mírou toho jak moc daný známý případ může být upřednostněn oproti neprozkoumaným oblastem. V oblasti mapy by mohla efektivita částečného případu záviset na nějakých faktorech jako jsou silniční podmínky nebo dopravním ruch. Efektivita vyhovuje nerovnosti β > 0 a může se měnit od případu k případu. Malé hodnoty β odpovídají vhodným ulicím. Hodnoty β > 1 odpovídají nevhodným ulicím. Předpokládejme, že cena projíždění neznámého regionu v rovině je rovna procestované vzdálenosti, kdežto cena projíždění ve známém případě je rovna β krát procestovaná vzdálenost. Zmíněná cesta bude spojitá jednoduchá cesta v rovině. Cesta může zahrnovat nějaké případové segmenty (nebo jejich části) a může také zahrnovat neprozkoumané oblasti. Přiřaďme každé cestě cenu která je sumou cen jejich částí. Problém hledání správné množiny případových segmentů je redukován na geometrický algoritmus ve kterém nalezneme optimální cestu (tj. cesta s nejmenší cenou) z počátečního vrcholu i do cílového vrcholu g v případovém grafu. Případové segmenty nalezené v nejkratší cestě jsou vráceny do SGP, který vyváří detailní plán který používá případů pro účely navigace. Na konci jakéhokoliv řešení problému jsou úspěšná řešení přidána do případové knihovny. S každým případem je spojeno řešení trasy, jenž poskytuje detailní popis plánu a oprávnění pro všechny vytvořená rozhodnutí. Pokud je případová knihovna prázdná nebo obsahuje případy které nejsou aplikovatelné v nových situacích, systém použije své současné doménové znalosti ke konstrukci schůdného řešení a postupně rozšíří případovou knihovnu. Podobně v ojedinělé situaci kdy systém selže pří generování řešení založeného na stejné množině případů, se může
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 48
vrátit přes toto rozhodnutí a použít doménové znalosti k nalezení schůdného řešení pro nový problém, jestliže alespoň jedna taková znalost existuje. 5.3.3 Přesný algoritmus pro nalezení optimální cesty
Pro řešení geometrického problému jej transformujeme do grafového problému. Definujme G jako úplný graf (každý pár vrcholů je spojen hranou) jehož vrcholy jsou vrcholy případového grafu (koncové body případů) plus počáteční a cílové body (i a g) . Každé hraně (u,v) grafu G je přiřazena cena jenž reprezentuje cestu správné jednoduché cesty mezi u a v. Na tento graf aplikujeme časově náročný ale přesný Dijkstrův algoritmus pro nalezení nejlevnější cesty z i do g v grafu G. Složitá část algoritmu je ve výpočtu ceny každé hrany. Tam kde jsou dva vrcholy spojeny případovým segmentem je cena jednoduše β násobeno Eukleidovská vzdálenost mezi vrcholy. Výpočet není vždy jednoduchý. Proč dokazuje obr. 28. Představme si, že yz je hlavní silnice a je rychlejší použít hlavní silnici na část cesty než přímo jít z x do y. Nejlepší místo na sloučení s hlavní silnicí je bod a.
Obr.28 Hledání cesty minimální ceny Na obr.28 je yz případový segment. Optimální cesta mezi x a y je (xa, ay) pro nějaký bod a takový, že závisí na hodnotě β. Cena optimální cesty je rovna délce xa+β násobeno délka ay. Nechť θ reprezentuje úhel xaz. Pozice bodu a je vypočtena z faktu že θ = cos-1β. V okrajovém případě kde β = 0 je xa kolmá na yz. P
P
Obr.29 Čtyři příklady nejkratších cest. Jsou možné i jiné případy. Body kde se spojují tečkované přímky s plnými závisí na cenovém faktoru β případového segmentu. Přiřaďme cenu každé hraně (x,y) grafu G s ohledem na počet možných cest mezi x a y a vezměme cestu s minimální cenou. Předpokládejme, že první cesta je úsečka mezi x a y. Potom pro každý případový segment předpokládáme optimální cestu která používá segment. Několik příkladů je demonstrováno na obr. 29 pro nejkratší cestu mezi x a y. První tři příklady jsou všechny řešeny rovnicí θ = cos-1β. Příklad 2 a 3 jsou pozoruhodné protože oni ukazují, že segment nepotřebuje být spojen s x nebo y aby poskytl nejlevnější cestu. Příklad 4 ukazuje, že segment mezi x a y vždy nezlepšuje optimální cestu. Pro každou hranu (u,v) v G zaznamenáme která jednoduchá cesta z u do v měla minimální cenu tak, že může být později rekonstruována úplná cesta (a množina případů). Třebaže předpokládáme pouze malý počet možných cest pro každou hranu, může být ukázáno že hledání nejlevnější cesty z i do g v grafu G je ekvivalentní s nalezením optimální cesty z i do P
P
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 49
g v případovém grafu. Bohužel, náročnost tohoto algoritmu je neúměrná, protože potřebujeme spojit každý pár vrcholů, dokonce i pokud jsou tyto páry vzdáleny od jiných bodů (ačkoli i a g mohou být velmi vzdálené, optimální cesta mezi nimi může být přímka). Graf G má O(n2) hran, kde n je počet vrcholů. To dává O(m) čas k určení ceny každé hraně kde m je počet případových segmentů. Z toho tedy vyplívá, že celková složitost je O(m3). Nicméně, můžeme navrhnout rychlejší algoritmus, který nalezne aproximaci optimální cesty. P
P
P
P
5.3.4 Rychlý algoritmus pro nalezení vhodné cesty
K odvození rychlého aproximativního algoritmu vezměme výhodu lokálnosti – princip který je nejvíce pravděpodobný pro cestování z vrcholů do jiných vrcholů které jsou poblíž. To nám ušetří výpočetní úsilí, protože můžeme ignorovat interakci mezi vzdálenými vrcholy a případovými segmenty. Delaunayho triangulace je definována kapitole 3.3 jako rovinný graf. Takovýto graf má pro účely návrhu tyto výhody: 1. Poskytuje strukturu která vytváří možnost rychle stanovit váhy cen grafu G. 2. Lokální modifikace triangulace je jednoduchá a rychlá. 3. Je to forma aproximace jejíž vrcholy jsou nejbližší ke každému jinému. Podívejme se na příklad na obr. 30. Ten ukazuje malou množinu bodů a jejich Delaunayho triangulaci. Představme si že každý z bodů a po h jsou koncové body případových segmentů a že i je počáteční bod. V každém směru okolo i triangulace indikuje který případ je nejblíže k i. Jiné případy vně šestiúhelníka se středem v i (takové jako jsou g a h) jsou více vzdálené a je daleko méně pravděpodobné, že budou přímo spojeny k i v koncové řešení cesty. Právě tyto více vzdálené případy v této heuristice ignorujeme.
Obr.30 Delaunayho triangulace množiny bodů Rychlý algoritmus se skládá z několika kroků. Prvně, vytvoříme DT vrcholů případového grafu. Potom vypočteme ceny hran a použijeme Dijkstrův algoritmus pro nalezení množiny případových segmentů. Krok 1: Přizpůsobování Delaunayho triangulace. Triangulace případového grafu vytváří možnost určení nejpravděpodobnějších případů k nimž bude bod spojen a dovolí nám redukovat počet případů zkoumaných pro každý pár vrcholů. Tento čin je však heuristika, to znamená, že ve velmi zvláštních případech jsou některé případové segmenty ignorovány i když by být neměly. K reprezentaci případového grafu použijeme přizpůsobenou DT. Vrcholy případového grafu jako byl na obr. 27 vlevo jsou triangulovány. Potom jsou přidány přídavné body které
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 50
zajistí, že všechny případové segmenty se budou jevit jako hrany v triangulaci (obr. 31) a že další podmínky, které tu nebudeme uvádět jsou splněny.
Obr.31 (a) Případový graf (b) DT případového grafu (c) Přizpůsobená DT případového grafu Krok 2: Cena hran Výpočet ceny pro hranu (x,z) je podobná jako u optimálního algoritmu, ale místo prozkoumávání každého případového segmentu prozkoumáme pouze hrany trojúhelníků v sousedství hrany (x,z). Může být ukázáno, že pouze segmenty které potřebují být prozkoumány jsou ty jejichž průsečík je uvnitř kružnice o průměru velikosti hrany (x,z).
Obr.32 Přizpůsobená DT a cesta nalezená Dijkstrovým algoritmem pro danou množinu hodnot β spojenou s každou hranou. Plné úsečky jsou případové segmenty čárkované úsečky jsou hrany triangulace. Obvykle potřebujeme pouze prozkoumat dva trojúhelníky přilehlé k hraně (x,z). (Výjimky nastávají pro málo případů jako jsou například případy 2 a 3 na obr. 29. Dobrá heuristika je taková která ignoruje tyto případy a prozkoumá pouze dva trojúhelníky). Pro dva trojúhelníky ukazuje některé možnosti obr. 33. K určení ceny hrany (x,z), musí být uvážen pouze malý počet alternativních cest zabírajících pouze konstantní čas na jednu hranu. Krok 3: Dijkstrův algoritmus Jakmile je triangulace a cenové přiřazení kompletní použijeme Dijkstrův algoritmus pro nalezení optimální cesty jako předtím. Můžeme ukázat, že v nejhorším případě cena této
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 51
cesty je menší než 3,5 faktor optimality. V praxi se tento algoritmus ukázal ještě mnohem lepší. Nikdy nenalezl cestu horší než 1,3 krát delší než optimální cesta. Cesta zobrazená na obr. 32 ukazuje cestu vypočtenou Dijkstrovým algoritmem pro daný počáteční a cílový bod a nějakou beta hodnotou přiřazenou každé hraně. Případové segmenty na cestě vrácené tímto algoritmem jsou poslány do SGP, jenž transformuje tento plán tak aby mohl být použitý robotem.
Obr.33 Výpočet ceny pro dva trojúhelníky
.
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 53
6
Implementace vybraného problému
6.1
Algoritmus inkrementální konstrukce a jeho časová složitost
Tento algoritmus je na implementaci nejjednodušší ze všech algoritmů pro konstrukci DT vyjmenovaných v kapitole 4. Při implementaci tohoto algoritmu jsem použil dva dvousměrné lineární seznamy do kterých se ukládají po řadě hrany hranice aktuální triangulace a hrany výsledné triangulace. Každý z těchto seznamů je implementován tak, že přidání nového prvku je provedeno v čase O(1). Uvažujeme–li dle kapitoly 3.3, že počet hran triangulace je ne= 3n - 3 - k a pro každou hranu je vyhledáván bod s nejmenší Delaunayho vzdáleností pro všechny body je pro množinu n vstupních bodů funkce f(n) definována rovnicí. B
B
f (n) = (3n − 3 − k )n = 3n 2 − 3n − kn kde n je počet vstupních bodů a k je počet hran konvexního obalu. Ovšem počet hran konvexního obalu předem neznáme a navíc je tento počet vzhledem k hranám triangulace zanedbatelný proto zjednodušíme tuto funkci na 3n 2 − 3n = 3(n 2 − n) ≈ n 2 − n , kde konstantu zanedbáváme a ze zbylého výrazu nás zajímá dominantní složka. Tedy odhadovaná funkce g(n) = n2 a časová složitost tohoto algoritmu bude tedy O(n2). Praktickým měřením tohoto algoritmu jsem však zjistil, že časová složitost je o něco menší než teoretická a to i bez použití P
P
P
P
Tabulka 1 n 50 100 200 400 800
t [ms] 10 20 120 631 2994
Časová složitost algoritmu InDeCo 3500 3000
Čas [t]
2500 2000 1500 1000 500 0 0
100
200
300
400
500
600
700
Bodů [n]
VUT v Brně, FSI ústav automatizace a informatiky
800
900
Delaunayho triangulace a její aplikace
Strana 54
jakýchkoliv urychlovacích technik. V tabulce 1 jsou výsledky naměřené na počítači s procesorem AMD Duron 1,1 Mhz a operačním systémem Windows XP Professional SP2. V tabulce jsou dva sloupce kde n je počet vstupních bodů a t je naměřený čas. B
6.2
B
Výpočet acw orientace trojice bodů
Mějme trojici bodů (P1, P2, P3). Tato trojice má orientaci v protisměru hodinových ručiček, jestliže P1 leží vlevo od hrany (úsečky) orientované od P2 k P3. Jestliže Pi = (Pix ,Piy) pak orientace P1, P2, P3 je dána znaménkem determinantu D následující matice: B
B
B
P1x D = P2 x P3 x
6.3
P1 y
B
B
B
B
B
B
B
B
B
1
P2 y 1 P3 y 1
Výpočet Steinerova bodu
Zkonstruování Steinerova bodu pro trojúhelník jehož všechny tři úhly jsou menší než 120° může být provedeno tak jak bylo uvedeno v kapitole 5.2.3. Na výpočet však jednodušší je konstrukce alternativní která je následující (viz. obr. 34). Pro daný trojúhelník ABC sestrojíme na jeho hranách 3 rovnostranné trojúhelníky všechny vně trojúhelníku ABC. Nyní jestliže spojíme vnější vrcholy rovnostranných trojúhelníků s protilehlými vrcholy trojúhelníka ABC dostaneme 3 úsečky a jejich společným průsečíkem je právě Steinerů bod. Pro účely řešení tohoto problému postačí zkonstruování dvou takovýchto přímek a vypočtení jejich vzájemného průsečíku.
Obr.34 Alternativní konstrukce Steinerova bodu
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 55
V implementaci máme při výpočtu k dispozici trojúhelníky které jsou reprezentovány jejich datovou strukturou, která je popsána dále v této práci. Tato struktura udržuje hrany každého trojúhelníka v jednom sledu orientace a samozřejmě udržuje informaci o jeho bodech. K sestrojení rovnostranného trojúhelníka na hraně daného trojúhelníka ABC potřebujeme tedy pouze najít třetí vrchol který není oběma trojúhelníkům společný. Způsob výpočtu tohoto bodu pro hranu daného trojúhelníku bude následující. Nechť AB je hrana ke které hledáme třetí bod rovnostranného trojúhelníku. Nejprve vypočteme souřadnice středu S této hrany. Jak je známo z analytické geometrie výpočet můžeme provést použitím rovnic. Sx =
Ay + B y Ax + B x ,Sy = 2 2
Vzdálenost hledaného bodu od hrany AB je rovna výšce v rovnostranném trojúhelníku proto tedy když máme k dispozici směrový vektor v hrany AB , který je dán rozdílem vx=Bx Sx a vy=By - Sy vypočteme souřadnice hledaného bodu jednoduše z rovnic kde u je směrový vektor kolmý na vektor v a |u| = |v|. B
B
B
B
B
B
B
B
B
B
B
B
Px = S x + u x 3 , Py = S y + u y 3 Tento postup použijeme i pro výpočet druhého rovnostranného trojúhelníka. Po tomto výpočtu obdržíme dva body které spolu s dvěma protilehlými body trojúhelníka ABC tvoří právě Simpsonovy úsečky (hrany). Najít průsečík těchto dvou hran je jednoduchá úloha z analytické geometrie která se sestává z dvou parametrických rovnic přímek (hran). Řešením těchto dvou rovnic obdržíme souřadnice Steinerova bodu v rovině pro daný trojúhelník. Při tomto způsobu výpočtu je opět důležité udržovat hrany všech trojúhelníků v jednom směru orientace aby bylo zaručeno, že konstruovaný rovnostranný trojúhelník bude vně daného trojúhelníka ABC.
6.4
Datové struktury
Datové struktury jsou tvořeny třídami. Při implementaci bylo třeba použít tří základních datových struktur pro reprezentaci geometrických útvarů a také lineárního seznamu. Těmito geometrickými útvary jsou bod, úsečka (hrana), a trojúhelník. Každá ze struktur obsahuje i další obslužné metody které jsou právě rozebrány v této kapitole. Datové struktury jsou až na lineární seznam obsaženy v souboru classes.h 6.4.1 Třída Vertex
Tato třída reprezentuje bod v prostoru. Vlastnosti třídy popisuje následující kód. Máme k dispozici tři souřadnice x, y, z typu double, vlastnost id která identifikuje bod podle přiděleného čísla. Třída umožňuje vytvoření nového bodu použitím konstruktoru který má jako svoje parametry souřadnice a id bodu nebo použitím konstruktoru který má jako svůj parametr ukazatel na již existující bod od kterého jsou parametry pro nový bod okopírovány.
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
class vertex { public: int id; bool mark; bool steiner; double x,y,z; vertex(); vertex(double ax, double ay, double vertex(PVERTEX a_vertex); };
Strana 56
az, int aid);
6.4.2 Třída Edge
Tato třída reprezentuje úsečku (hranu). K dispozici je i několik metod které provádí nezbytné operace jako je převracení a výpočet délky (ceny). Hrana je složena ze dvou bodů A a B. Metoda swap() umožňuje změnu orientace hrany z AB na BA a naopak. Metoda compute_length() vypočte délku hrany tedy Eukleidovskou vzdálenost, která je uložena do vlastnosti length této třídy a je přístupná skrze metodu get_length(). Metody orig() a final() vrací po řadě počáteční a koncový bod orientované hrany. Další důležitou metodou je metoda is_my_point(PVERTEX P), která kontroluje zda zadaný bod jako parametr je bodem této hrany. Pokud ano vrací ukazatel na druhý bod této hrany jinak vrací NULL. Třída má dva konstruktory které umožňují vytvářet novou hranu buďto voláním konstruktoru s dvěma parametry které jsou body nebo jako kopii již stávající hrany voláním konstruktoru spolu s ukazatelem na již existující hranu.
class edge { private: PVERTEX A; PVERTEX B; double length; double compute_length(); public: bool mark; PVERTEX orig() {return A;} PVERTEX final() {return B;} double get_length() {return length;} void swap(); PVERTEX is_my_point(PVERTEX P); edge(PVERTEX vertex_a, PVERTEX vertex_b); edge(PEDGE a_edge); ~edge(); };
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 57
6.4.3 Třída Triangle
Třída Triangle reprezentuje trojúhelník. Tato třída zapouzdřuje tři hrany jenž reprezentují trojúhelník. Máme k dispozici také velikosti úhlů, které jsou přilehlé k jednotlivým vrcholům a také vlastnost obtuse která značí zda je některý z těchto úhlů tupý. Hrany trojúhelníka jsou orientovány v jednom směru tak jak je vidět na obr. 35 a proto můžeme kdykoliv zjistit díky metodám get_A(), get_B() a get_C() všechny tři ukazatele na vrcholy trojúhelníka. Metoda check_angles() vypočte jednotlivé úhly a v případě tupého úhlu nastaví vlastnost obtuse. Další důležitou metodou je get_obtuse_edges(PEDGE &a, PEDGE &b), která vrací dva ukazatele na hrany jenž jsou přilehlé k tupému úhlu. Tyto hrany jsou však ukazatele na nově vzniklé kopie původních hran, které však nesou původní ukazatele na body. Funkce ostatních metod jsou zřejmé z výpisu kódu.
Obr.35 Reprezentace trojúhelníka pomocí třídy Triangle class triangle { private: PEDGE a; PEDGE b; PEDGE c; double alfa,beta,gamma; bool obtuse; void check_angles(); bool acw(PVERTEX A, PVERTEX B, PVERTEX C); public: double get_alfa() {return alfa;} double get_beta() {return beta;} double get_gamma() {return gamma;} triangle(PEDGE a, PEDGE b, PEDGE c); triangle(PVERTEX A, PVERTEX B, PVERTEX C); PVERTEX get_A() {return c->orig();} PVERTEX get_B() {return a->orig();} PVERTEX get_C() {return b->orig();} PEDGE get_a() {return a;} PEDGE get_b() {return b;} PEDGE get_c() {return c;} void get_obtuse_edges(PEDGE &a, PEDGE &b); bool is_obtuse() {return obtuse;} bool mark; };
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 58
6.4.4 Třída List
Jedná se o třídu jenž implementuje dvousměrný lineární seznam. Seznam je implementován tak že je možno ho použít pro ukládání dat libovolného typu. Pro přidávání nové položky slouží metoda add(TDATA data) jenž jako parametr má ukazatel na data které mají být do seznamu přidána. Metody reset(), next() a prev() po řadě provádějí nastavení seznamu do výchozí polohy (seznam ukazuje na první prvek), přechod na další prvek a přechod na předchozí prvek. Metody delete_item() a delete_item(TDATA data), umožňují smazání aktuálního prvku na který je seznam nastaven a nebo vyhledání a smazání prvku který nese ukazatel na data předané jako parametr. Jsou zde také funkce pro smazání seznamu jako celku. Funkce clear_all() smaže jak jednotlivé položky seznamu tak také data na které tyto položky ukazují. Funkce clear_list() maže pouze položky seznamu nikoliv však data na které ukazují. Pokud tedy list nějaká data obsahuje (resp. ukazatele na tyto data) budou nenávratně ztracena. Poslední z této řady je funkce clear_data() která maže pouze data na které ukazují položky seznamu.
6.5
Problém v návrhu algoritmu inkrementálního vkládání
V průběhu řešení této práce vznikl problém s návrhem algoritmu inkrementálního vkládání jehož časová složitost je O(n log(n)). Dle [2] a tak jak je tento algoritmus popsán v kapitole 4 je funkce pro testování legálnosti hrany upravena tak aby byly ošetřeny případy kdy do této funkce vstupuje hrana jenž obsahuje body velkého počátečního trojúhelníka a tedy nepatří do výsledné DT vstupní množiny bodů. Funkce je rozdělena na několik případů podle nichž se o legálnosti hrany rozhodne. Podívejme se na příklad uvedený na obr. 36. Nechť Pi Pj je hrana určená k testování a Pk, Pl jsou body které spolu s hranu tvoří dva B
B
B
B
přilehlé trojúhelníky. Přidání bodu P1 proběhlo bez jakýchkoliv problémů, protože hrany určené k testování jsou v tomto případě P−1 P− 2 , P− 2 P−3 a P−3 P−1 a tedy se jedná o případ 1 kdy indexy u bodů hran jsou negativní a proto rozhodneme, že hrana je legální a tedy nebude prohozena. Problém ovšem nastane hned při přidávání nového bodu Pr . Po přidání tohoto bodu se vytvoří tři nové trojúhelníky a hrany k testování budou tedy P1 P− 2 , P− 2 P−1 a P1 P−1 . Při B
B
B
B
testování hrany P− 2 P−1 by neproběhla výměna hrany jelikož se jedná o stejný případ jako je uveden výše. Dle případu 4 pokud je záporný index z dvojice bodů Pi, Pj menší než záporný index z dvojice bodů Pk a Pl nedojde k prohození jinak ano. Jelikož máme situaci kdy hrana k testování je P1 P−1 a zbylé dva body jsou Pr a P-3 nastane právě tento případ kdy se tato hrana prohodí. Toto prohození je však nepřípustné jelikož by v tomto kroku vytvořil nekonvexní čtyřúhelník P−1 P−3 P1 Pr . B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 59
Obr.36 Problém při vkládání bodu Pr B
B
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
6.6
Strana 60
Úprava heuristického algoritmu pro řešení EStMT
Při implementaci algoritmu pro řešení Eukleidovského Steinerova minimálního stromu byla provedena úprava tohoto algoritmu který je popsán v [3] a také v závěru kapitoly 5.2.3 této práce. Mějme algoritmus ve fázi po provedení kroku č.2 a 3 kdy jsou všechny trojúhelníky jenž mají všechny úhly větší než 120° nahrazeny Steinerovým stromem a poté je na takto vytvořený graf aplikován algoritmus pro určení minimální kostry. Ve výsledném grafu se vyskytují Steinerovy body stupně 1. Tuto situaci znázorňuje obr. 37 kde neterminály jsou zobrazeny červenou barvou a terminály modrou barvou. Tyto body jsou však nežádoucí pro další vylepšování které následuje v kroku 4 a jsou také nežádoucí ve výsledném grafu. Mimoto se v tomto grafu který je vytvořen krokem 3 vyskytují také Steinerovy body stupně 2. Toto spojení je tvořeno hranami trojúhelníka který vystupoval v grafu který byl výsledkem kroku 2 před tím než byl aplikován algoritmus minimální kostry. Tyto body jsou také nežádoucí jelikož se evidentně jedná o takzvanou trojúhelníkovou nerovnost která nám říká, že součet délek dvou stran v trojúhelníku není nikdy kratší než strana třetí. Proto tedy je nutno také odstranit tyto body stupně 2 spolu s incidentními hranami a nahradit je přímou spojnicí. Ovšem po následujícím kroku 4 kdy jsou hrany analyzovány a optimalizovány opět dostaneme pro některé body výše zmíněný případ. Z tohoto důvodu jsem po analýze a testování tohoto algoritmu provedl jeho úpravu jenž ošetřuje výše zmíněné případy. Popis takto upraveného algoritmu je následující (viz. strana 61).
Obr.37 Příklad grafu obdrženého po kroku č.3.
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 61
1. Najdi Delaunayho triangulaci pro danou množinu V. 2. Nahraď všechny trojúhelníky které mají každý úhel menší než 120° EStMT stromy. 3. Urči EMSpT pro graf nalezený v kroku 2. 4. Odstraň všechny Steinerovy body stupně 2 spolu s incidentními hranami. Tyto hrany nahraď přímou spojnicí. 5. Odstraň všechny Steinerovy body stupně 1 dohromady s jejich incidentními hranami. 6. For každou hranu spojující body X a Y do 7. Nalezni hranu YZ která se setkává s XY v nejmenším úhlu. 8. If úhel je menší než 120° then 9. Umísti nový Steinerův bod Sn do trojúhelníku daného hranou XY a YZ . B
B
10. Odstraň hrany XY a YZ 11. Přidej hrany XS n , YS n a ZS n . 12. Odstraň všechny Steinerovy body stupně 2 s incidentními hranami. Tyto hrany nahraď přímou spojnicí. 13. Odstraň všechny Steinerovy body stupně 1 dohromady s jejich incidentními hranami.
6.7
Výsledky získané upraveným algoritmem pro určení EStMT
Po výše zmíněné úpravě jsem provedl měření výsledků pro pět různých velikostí vstupních množin bodů spolu s porovnáním délky minimální kostry. Naměřené hodnoty lze porovnat v tabulce 2 kde n je velikost vstupní množiny bodů, w(EStMT) je cena Steinerova minimálního stromu pro danou množinu bodů, w(EMSpT) je cena minimální kostry pro danou množinu bodů a v posledním sloupci je rozdíl těchto cen vyjádřený v procentech.
[n] 5 10 20 50 100
Tabulka 2 w(EStMT) w(EMSpT) 567 604 1084 1106 1795 1873 2517 2549 3508 3556
Úspora [%] 6,125828 1,98915 4,164442 1,255394 1,349831
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
7.
Strana 63
Závěr
Po prostudování daného tématu jsem dospěl k závěru, že se jedná o zajímavou tématiku ovšem nenalezl jsem v dostupné literatuře žádné další významnější použití Delaunayho triangulace na některou technickou aplikaci. Výsledkem práce je tedy stručný popis jednotlivých algoritmů a některých příkladů použití Delaunayho triangulace a také softwarový projekt řešený jako součást Diplomové práce ve kterém byl implementován algoritmus pro vyhledání Steinerova minimálního stromu jenž používá Delaunayho triangulaci.
VUT v Brně, FSI ústav automatizace a informatiky
Delaunayho triangulace a její aplikace
Strana 65
8.
Literatura
[1]
Šeda M. Algoritmy pro konstrukci Delaunayho triangulace a porovnání jejich časové složitosti. Brno, 2004. 51 s., Diplomová práce na Fakultě strojního inženýrství na ústavu automatizace a informatiky. Vedoucí diplomové práce Doc.RNDr.Ing. Miloš Šeda PhD.
[2]
BERG, de M.; KREVELD, van M.; OVERMARS M.; SCHWARZKOPF O.: Computational geometry: Springer 2000, 367 s. ISBN 3–540–65620–0.
[3]
Šeda M. Solving the Euclidean Steiner Tree Problem Using Delaunay Triangulation. In WSEAS TRANSACTIONS ON COMPUTERS, page 471, Issue 6, Vol.4, June 2005, ISSN 1109–2750.
[4]
Textoris J.: Automatické generování sítí metodou Delaunay. Brno, 1999, 66 s. Diplomová práce na Fakultě strojního inženýrství na ústavu matematiky. Vedoucí diplomové práce RNDr. Rudolf Hlavička CSc.
[5]
TÓTH Z.: Rekonštrucia obrazu pomocou triangulacie, 2004, 49 s., Diplomová práce na fakultě Matematiky, Fyziky a Informatiky Univerzity Komenského v Bratislave. Vedoucí diplomové práce: RNDr. Andrej Ferko PhD.
[6]
HAIGH, Zita K.; SCHEWCHUK, Richard J. Geometric Similarity Metrics for Case– Based Reasoning, Working Notes from the AAAI-94 Workshop (Seattle, Washington), pages 182-187, AAAI Press, August 1994.
[7]
KOHOUT, J. Paralelní Delaunayova triangulace ve 2D a 3D, Plzeň, 2002, 85 s., Diplomová práce na Západočeské univerzitě v Plzni na Fakultě aplikovaných věd Katedře Informatiky a Výpočetní techniky. Vedoucí diplomové práce: Doc. Dr. Ing. Ivana Kolingerová.
[8]
VIRIUS, M.; Od C k C++, Kopp 2000, 227 s., ISBN 80–7232–110–2.
VUT v Brně, FSI ústav automatizace a informatiky