Výpočet Voroného diagramu v maximové metrice

Vstupem algoritmu je seznam bodů v rovině a seznam úseček mezi nimi. Úsečky se mohou protínat jen ve svých koncových bodech a žádné dva body nesmí mít stejnou svislou souřadnici. Souřadnice bodů jsou racionální čísla, jejich jmenovatele považujeme za omezené konstantou.

Výstupem je Voroného diagram v maximové metrice (L_∞) pro množinu ze vstupu. To je seznam úseček a polopřímek, souřadnice vrcholů jsou opět racionální čísla. Velikost výstupu, jak se ukáže později, je omezená lineárně velikostí vstupu.

Obecně uznávané pojmy

Základem je zobecněný Fortunův zametací algoritmus [1], v tomto odstavci ho zkusím popsat hodně stručně. Vstupní množinu bodů zametáme zleva doprava svislou přímkou (sweepline), při tom postupně nastává nějaký sled událostí. Uvažujeme hranici Voroného diagramu mezi sweepline a podmnožinou vstupu, která je od sweepline vlevo; tato hranice je vždy nějaká neomezená křivka, v anglických článcích se označuje jako wavefront. Všechny vrcholy výstupního diagramu pak budou patřit k událostem, kdy se wavefront významným způsobem změní (a to budou, konec konců, všechny uvažované události).

Voroného region daného bodu je množina bodů, které od něj jsou ostře blíž než od kteréhokoli jiného prvku (koncového bodu, izolovaného bodu nebo vnitřku úsečky) na vstupu. Obdobně region úsečky je množina bodů, které jsou ostře blíž od některého vnitřního bodu této úsečky než od všech ostatních prvků vstupu. Hrana Voroného diagramu je průnik uzávěrů dvou různých regionů, pokud je neprázdný. Vrchol Voroného diagramu je průsečík některých hran, pokud je to bod.

V maximové metrice je množina bodů ve vzdálenosti a od zadaného bodu právě obvodem čtverce se středem v tomto bodu a hranou 2·a, přičemž hrany čtverce jsou rovnoběžné s osami souřadnic. Hranu Voroného diagramu mezi dvěma prvky vstupu pak jde přirozeně popsat také jako středy čtverců, jejichž hranice protíná právě zadané dva prvky, a přitom jejichž vnitřek neprotíná žádný prvek vstupu.

Doplňující definice

Uvedené definice teď doplním a přiohnu. Když hranici čtverce protíná víc prvků najednou, budu za ostře nejbližší středu považovat ten z nich, který má průsečík nejblíže levému hornímu nebo pravému dolnímu rohu čtverce. To pomáhá rozseknout problémy především se svislými úsečkami, jinak by diagram mohl ošklivě degenerovat. Výsledek odpovídá tomu, jako bychom měřicí čtverec mírně protáhli podél osy x + y = 0.

Díky předchozí definici má libovolné okolí každého bodu roviny průnik s některým regionem diagramu, takže uzávěr všech regionů je celá rovina. Každý z regionů je navíc souvislý díky tomu, že se úsečky neprotínají.

Programátorskému pohodlí prospěje každý region podrozdělit podle toho, která část čtverce se příslušného prvku dotýká. V případě izolovaného bodu uvažujeme jen dotyky hranou čtverce, které jsou tedy čtyři možné: kdyby byl na vstupu jen jeden bod, jeho Voroného diagram by sestával ze čtyř diagonálních polopřímek. V případě úsečky naopak uvažujeme jen dotyky v rozích čtverce. V souladu s předchozím odstavcem může mít svislá úsečka jen regiony odpovídající dotyku levým horním a pravým dolním bodem. Úsečka je sama hranou svých regionů, a každý z jejích koncových bodů může mít až dva regiony, podle dalších okolností.

Pro shrnutí: výstupní Voroného diagram bude sestávat z regionů odpovídajících dotyku určitého bodu vstupu levou hranou čtverce, resp. pravou, horní, dolní, a z regionů odpovídajících dotyku určité úsečky vstupu levým horním bodem čtverce, resp. levým dolním, pravým dolním, pravým horním. Jednomu bodu ze vstupu můžou náležet žádný, jeden, dva nebo čtyři regiony, podle toho, jaké úsečky v něm končí. Jedné úsečce náleží vždy dva regiony odpovídající protilehlým rohům čtverce.

Z povahy maximové metriky plyne, že hranice dvou regionů (dle nové definice) je vždy podmnožinou polopřímky. Její parametry totiž lze explicitně spočítat.

Konečně, někdy se bude hodit seřadit body na společné svislé ose, aby události nastávaly v jednoznačném pořadí. Volím pořadí zdola nahoru, pohodlná představa je, že je sweepline mírně pootočená proti směru hodinových ručiček.

Wavefront

Wavefront je levá hranice Voroného regionu, který v daný okamžik náleží ke sweepline. V případě maximové metriky je to lomená čára. Každá z jejích hran (tj. úseček, ze kterých je wavefront složená) dělí region náležící sweepline od právě jednoho regionu některého prvku vstupu, kterýžto prvek je aspoň částečně nalevo od sweepline. Ten vztah ale není symetrický: jeden region může být ve wavefront zastupovaný libovolným počtem úseček, případně také žádnou. Region, který je zastoupený aspoň jednou úsečkou, budu nazývat jako aktivní.

Každému bodu wavefront odpovídá jeden měřicí čtverec maximové metriky. O tom čtverci víme, že se dotýká sweepline, takže ho jednoznačně popisuje jeho spodní a horní hranice. Wavefront je topologicky ekvivalentní přímce, body na ní můžu chápat jako uspořádané odspoda nahoru. Navíc pro libovolné dva body na wavefront platí, že ten, který je v pořadí později, má horní i spodní hranici měřicího čtverce neostře větší než odpovídající hranici toho dřívějšího z nich. Nemůžou mít ale obě hranice navzájem stejné, protože jsou to různé body. Díky tomu jde jednoduchými výpočty pořadí bodů na wavefront zjišťovat.

V programu je wavefront reprezentovaná jako seznam regionů v pořadí, jak odpovídají jejím bodům zdola nahoru při daném umístění sweepline. Regiony jako takové uspořádat nejde, protože každý se může v pořadí vyskytovat víckrát. Mezi každými dvěma sousedními regiony se ale wavefront láme v bodu, jehož polohu jde přesně spočítat. To v případě potřeby umožňuje vyhledávat sousední dvojici půlením intervalů (v programu funkce square_top, square_bottom)

Spike

Hranice mezi dvěma sousedními regiony může protínat wavefront. V takovém případě tu hranici (přesněji, dotyčnou její souvislou část) nazvu jako spike. Vychází ze svého levého krajního bodu ve směru, který dva regiony jednoznačně určují. Jedinou výjimkou je hranice mezi dvěma regiony na obou stranách jedné úsečky, tedy ta úsečka samotná; tu jako spike nenazývám, protože by to spíš přidělalo potíže. Každé jiné dvojici sousedních aktivních regionů náleží právě jeden spike. V programu si ho oba příslušné regiony přímo pamatují.

Spike jsou důležité jen proto, že se můžou někde protínat. Průsečík dvou sousedních spike odpovídá v diagramu rohu regionu, který při zametání zmizí z wavefront. To je potřeba ve vhodnou chvíli zpracovat, aby seznam aktivních regionů neustále odrážel jejich skutečné pořadí podél wavefront.

Spike můžou směřovat i doleva, a to za některých okolností až diagonálně (tj. s derivací +-1), podél wavefront. V programu s tím není problém, jen je dobré si tu možnost uvědomit.

Nekonečno

Výsledný Voroného diagram bude obsahovat nějaké polopřímky. Všechny jsou ale diagonální, mají čtyři možné směry; to jde ověřit rozebráním možných sousedících regionů. V programu s nimi proto pracuju jako s obyčejnými úsečkami, které mají jeden koncový bod v některém z rohů plochy. V programu se nicméně pro pohodlnější kreslení smí každý roh vyskytovat v mnoha instancích s různými parametry.

Celý program je zamýšlený pro praktické použití v konečných podmínkách, kdy bude rovina omezená nějakým obdélníkem. Ošetřování zvláštních případů se podle mě zpřehlední přidáním magického regionu Border, který se dá chápat oblast mimo tuto rovinu. Ve wavefront se vyskytuje na jejím horním a dolním konci (představa může být taková, že mu náleží body někde v nekonečnu). Mezi ním a sousedními aktivními regiony jsou obyčejně spike, které vedou z příslušných rohů na levé straně přímo vodorovně. Průsečíky s nimi leží opět v rozích, ale jinak se s nimi pracuje normálně.

Zdá se mi intuitivní, že Border je na začátku běhu programu jediným aktivním regionem v seznamu.

Událost přidání

Aktivní regiony můžou přibýt jedině v okamžiku, kdy sweepline překročí nějaký bod vstupu. Budou to regiony toho nového bodu nebo úseček, které z něj vedou doprava. Jiné regiony zároveň můžou přestat být aktivní: bude to proto, že je nový bod nebo úsečky vedoucí do nového bodu zleva zastíní. Může se stát, že jeden region úsečky se v tom okamžiku rozdělí vejpůl, pokud nový bod zastíní jen prostřední část úsečky. Slovem stínit se rozumí překážet ve vnitřku všech měřicích čtverců, které zastíněnému regionu ve wavefront odpovídaly.

Po přidání bodu se změní nějaká souvislá část wavefront, a zastíněné regiony proto tvoří nějaký souvislý podseznam aktivních regionů. V programu se určí tak, že najdeme pod ním první region, který zastíněný nebude, a obdobně první nezastíněný region nahoře. Může se stát, že to budou dva sousední regiony (pak není zastíněné nic), nebo že dokonce oba nalezené budou týž region: pak je potřeba ho patřičně rozdělit vedví, aby dovnitř šlo vložit regiony odpovídající novému bodu.

Ještě podrobněji k hledání zastíněné části před přidáním regionů nového bodu. Uvážíme horní hranici měřicího čtverce daného dvěma sousedními aktivními regiony a sweepline: buďto je neostře nad nově přidávaným bodem, nebo ostře pod ním. Prvním spodním nezastíněným regionem bude ten nejvyšší aktivní region, který se svým spodním sousedem dává hranici čtverce ostře pod novým bodem. Intuitivně odpovídá tomu, že ukotvíme měřicí čtverec za jeho pravý horní vrchol v novém bodu a zvětšujeme ho, dokud se nedoktne nějakého prvku vstupu: region tohoto prvku hledáme.

Obdobně při hledání prvního horního nezastíněného regionu uvážíme dolní hranici měřícího čtverce pro dva sousední aktivní regiony. Hledáme poslední region shora, který se svým vrchním sousedem dává tuto hranici ostře nad novým bodem.

Kdyby nový bod měl stejnou svislou souřadnici jako nějaký jiný, a ten měl nějaký aktivní region, takhle popsané hledání by dávalo špatné výsledky. Není v něm totiž zabudované zkosení měřicího čtverce. Je mi to líto, ale zatím se mi nepovedlo najít funkční a elegantní řešení toho problému, takže body se stejnou svislou souřadnicí radši zakazuju. Dlužno podotknout, že [2] tenhle problém ani nezmiňuje.

Sousedící zastíněné regiony můžou spolu s levým regionem nového bodu, pokud tento vůbec existuje, dávat vrcholy diagramu. Bod má levý region jedině v případě, že do něj nevede žádná úsečka zleva. Hranice mezi prvním nezastíněným regionem zdola a sousedním zastíněným může dávat vrchol diagramu se spodním regionem nového bodu, pokud tento existuje, nebo se spodním regionem nejnižší úsečky vedoucí z nového bodu doprava, pokud do nového bodu nevede žádná úsečka zleva zespodu. Obdobná pravidla platí pro horní situaci. Přehlednější je tentokrát asi zápis v programovém kódu.

Seznam nově přidávaných regionů považuju za podobně nepřehledný rozbor případů. Každé úsečce vedoucí doprava budou náležet dva nové regiony, a nově přidanému bodu můžou náležet jeden až tři nové regiony podle toho, jakými směry z něj úsečky vedou. Mezi každými dvěma regiony, které budou nově sousedět, je pak ještě potřeba dopočítat spike a k nim příslušné průsečíky.

Událost odebrání

Aktivní region může zmizet také jinde volně v rovině. Jeho horní a spodní soused potom začnou ve wavefront sousedět přímo, takže se v tom místě musí protnout příslušné spike. Bod průsečíku odpovídá středu čtverce, který se dotýká prvků příslušných všem třem zmíněným regionům. Okamžik, kdy je potřeba prostřední z nich odebrat ze seznamu aktivních, nastane, když wavefront protne průsečík spike: v tu chvíli bude sweepline obsahovat pravou hranu zmíněného čtverce.

V programu tedy nezbývá než počítat průsečíky každých dvou sousedních spike a přidávat je se správnou souřadnicí do seznamu událostí. Každý spike si pamatuje právě jeden průsečík, z obou možných ten, který nastane dříve. Když některý aktivní region zmizí, oba jeho spike přestanou dávat smysl a zrovna tak průsečíky, které si tyto spike pamatovaly. Výpočtem může vzniknout průsečík takový, že by si ho nepamatoval žádný spike: takový průsečík nedává smysl už od začátku. Průsečíky, které nedávají smysl, ze seznamu událostí natrvalo vyřadíme.

Průsečíky zpracováváme v zametacím pořadí spolu se vstupními body. Kdyby průsečík nastal ve stejném okamžiku jako nějaký bod vstupu, zpracujeme dřív průsečík, aby později vyhledávání fungovalo správně.

Průsečíky jednoho spike vedle objektu Border a sousedního spike směřujícího někam doprava mají souřadnici nastavenou na nekonečno. To prostě znamená, že budou vždy na konci seznamu událostí. V programu jsem se chtěl používání nekonečna místo čísel co možná vyhnout, takže tyto průsečíky ignoruju a až nakonec si je domyslím ze zbývajících spike.

Při zpracování průsečíku je potřeba dotyčné regiony najít v seznamu aktivních. Nejlepší je si jejich pozici v seznamu přímo pamatovat, což by datové struktuře nemělo dělat potíže. V přiloženém programu místo toho používám opět hledání pomocí horní a spodní hranice měřicího čtverce, což ale může selhat, pokud má v jednom bodě nastat víc průsečíků.

Zbytek zpracování průsečíku je jednoduchý: do diagramu za něj přidáme vrchol a nově sousedním regionům dopočítáme společný spike, který bude z průsečíku vycházet.

Výpočty

Asi hlavní výhodou maximové metriky oproti euklidovské je, že vzdálenosti bodů v racionálních souřadnicích mají racionální hodnoty. Celý program proto pracuje s racionálními čísly, čímž se vyhne všem zaokrouhlovacím chybám. Polohu každého bodu grafu jde vypočítat konstantně velkou soustavou lineárních rovnic. Délka jmenovatele každého bodu výstupu bude díky tomu lineárně omezená nejdelším jmenovatelem na vstupu. Tatáž argumentace platí o všech mezivýsledcích používaných v programu.

To by samo o sobě mohlo vést ke kvadratické složitosti, pokud by nejdelší jmenovatel byl lineárně dlouhý vzhledem k celému vstupu. Byla by škoda tím pošpinit jinak rychlý algoritmus, takže jmenovatele prohlásíme za omezené konstantou. Potom trvá každá početní operace konstantní čas a výstup je lineárně dlouhý vzhledem ke vstupu.

Asi hlavní nevýhodou maximové metriky je, že výpočty vyžadují dlouhé rozbory případů. Nejběsnější funkce v programu se jmenuje square_tuple, kterážto počítá délku strany a polohu měřicího čtverce pro dva dané aktivní regiony a danou polohu sweepline. O něco stručnější, ale hůře pochopitelný, je výpočet směrnice spike (Spike.__init__). Pro zestručnění jsem rozdělil regiony do tří skupin, protože potřebné výpočty se pro členy téže skupiny chovají stejně. Pro pobavení přikládám fotku jednoho z papírů, které jsem popsal při vymýšlení těch skupin.

Datové struktury

Nejdůležitější strukturou je seznam aktivních regionů. Měl by to být ve skutečnosti vyvážený strom. Potřebné operace jsou nalezení prvku zhruba v polovině daného intervalu, nalezení obou sousedů daného prvku a vložení prvku na dané místo, což běžné stromy bez potíží zvládají v logaritmickém čase. V programu se použije strom z modulu blist, pokud je dostupný, jinak ho nahradí obyčejný seznam (na čas výpočtu to nemá velký vliv). Program ale stejnak nepoužívá hledání bisekcí, takže mu to na rychlosti neprospěje.

Seznam událostí nemusí být úplně utříděný, stačí mu dát strukturu haldy. Když je nějaký průsečík prohlášen za nesmyslný, stačí ho rovnou z haldy odebrat (v programu to tak ale není).

Výstupní graf se dá považovat za datovou strukturu sám o sobě. Během zametání do něj nahodile házíme vrcholy diagramu, kdy každý je popsaný aspoň třemi regiony v pevném pořadí. Na požádání je možné vrcholy podle sousedících regionů pospojovat hranami, aby šlo výstup zapsat jako OBJ soubor. Soumístné vrcholy se při tom sloučí, i kdyby mezi nimi třeba hrana nevedla. Vespodu všechny operace vyřizují hashovací tabulky.

Složitost

Zpracování každé události trvá logaritmicky dlouho podle počtu událostí v haldě, protože je potřeba ji vytáhnout, a logaritmicky dlouho podle počtu aktivních regionů za nějaké konstantní množství operací se stromem.

Událostí přidání je právě tolik, co vstupních bodů. Operací odebrání nastane nejvýš tolik, co je v diagramu regionů, a to je lineárně omezené velikostí vstupu. Počet událostí odebrání v haldě (pokud nesmyslné průsečíky mažeme) je menší než počet aktivních regionů, protože každým dvěma sousedům patří nejvýš jeden průsečík. Počet aktivních regionů je nejvýš tolik, co regionů celkem.

Když trvají všechny potřebné početní operace konstantně dlouho, doběhne algoritmus v čase n·log(n).

Ukázkový program

Skript voro.py je psaný pro Python 2.5 až 2.7. Závisí na standardní knihovně heapq, která nabízí běžné operace s haldou, a na nějaké knihovně pro práci s racionálními čísly, kterou případně může nahradit přiložený zdroják mympq.py.

Prvním parametrem skriptu je jméno vstupního souboru, druhým jméno výstupního. Oba jsou ve formátu Wavefront OBJ (shoda jmen náhodná...), souřadnice bodů v nekonečnu jsou na výstupu zarovnány na okraj určitého obdélníka. Skript ze sebe sype hodně ladicích hlášek o všem, co dělá. Když je vstup nekorektní, skript prostě sletí.

Hlavní řídící kód je v posledním odstavci skriptu. Odkomentováním jiné varianty řádku volajícího graph.clean_output je možné zapnout dodatečné zpracování, kdy se na výstupu ponechají jen hrany oddělující nesouvislé části vstupu. Skript má sloužit k řezání plošných spojů, kde je právě tohle potřeba.

Několik testovacích vstupů je přibalených. Další je možné buďto napsat přímo v textovém editoru jako OBJ soubor, nebo nakreslit ve vektorovém editoru jako SVG a převést přibaleným skriptem svg2obj.py. Ten se volá také se dvěma parametry, po řadě název vstupního a výstupního souboru. SVG je ale matoucí tím, že má osu y směrem dolů, konverzní skript souřadnice neopravuje.

Dále je přibalený skript obj2svg.py pro pohodlné prohlížení výsledků. Volá se stejně. Hranice dokumentu jsou nastavené na stejný rámeček, na jaký se zarovnávají body z nekonečna.

Pojmenování objektů a proměnných je v kódu trochu jiné než v tomhle popisu. Region a jemu příslušný prvek vstupu tam nazývám dohromady jako element. Seznam aktivních elementů nazývám jako sweepline. Proměnné element_bottom, element_top jsou po řadě první spodní a vrchní nezastíněný element. Body vstupu se nazývají site, což je v souladu s odkazovanými články.

Obrázky

Reference