Přechod na menu, Přechod na obsah, Přechod na patičku

Diagramy Voronoia

Úvod

Diagramy Voronoia podávají řešení problému známého pod anglickým názvem „post office problem“. V něm jde o to rozdělit teritorium města na oblasti okolo jednotlivých pošt tak, aby v každé oblasti byla místní pošta ta nejbližší. Jiným příkladem je určení spádových oblastí nemocnic podle vzdálenosti. Matematická formulace problému je následující:

V rovině je dána množina bodů $P=\{p_1,p_2,\ldots,p_n\}$. Diagram Voronoia (zkráceně V-diagram) je rovinné podrozdělení, které má oblasti

$$V(p_i)=\{q\in \mathbb{R}^2;\ \operatorname{dist}(q,p_i)\leq \operatorname{dist}(q,p_j)\ \text{pro všechna}\ j\in P\setminus\{i\} \},$$

kde $\operatorname{dist}(q,p)$ označuje klasickou vzdálenost bodů $q$ a $p$.

Obrázek 9.1
Obrázek 9.1
Příklad diagramu Voronoia.

Množinou bodů v rovině, které mají vzdálenost k bodu $p_i$ menší nebo rovnu vzdálenosti k bodu $p_j\ne p_i$ je polorovina

$$h(p_i,p_j)=\{q\in \mathbb{R}^2;\ \operatorname{dist}(q,p_i)\leq \operatorname{dist}( q,p_j)\}.$$

Oblast Voronia $V(p_i)$ je tedy rovna průniku všech takových polorovin pro $j\ne i$

$$h(p_i,p_j)=\bigcap_{j\in P\setminus\{i\}}h(p_i,p_j).$$

Průnik $n-1$ polorovin umíme podle 5. kapitoly najít v čase $O(n\log n)$. Najít tímto způsobem všechny oblasti by vyžadovalo čas aspoň $O(n^2)$. My si ukážeme algoritmus pro nalezení diagramu Voronoia s časovou náročností $O(n\log n)$. O tom, že taková časová náročnost je reálná, svědčí následující tvrzení:

Věta 9.1

Diagram Voronoia pro $n\geq 3$ bodů má nejvýše $2n-5$ vrcholů a nejvýše $3n-6$ hran.

Důkaz

Označme $m$ počet vrcholů a $h$ počet hran V-diagramu. Přidáním vrcholu $v_{\infty}$ v „nekonečnu“, do kterého vedou všechny hrany V-diagramu, které jsou polopřímky, vytvoříme z V-digramu rovinný graf, který má $n$ oblastí, $m+1$ uzlů a $h$ hran. Viz obrázek.

Obrázek 9.2
Obrázek 9.2
Doplnění V-digramu na planární graf.

Pro něj platí Eulerova věta

$$m+1-h+n=2.$$

Stupeň každého uzlu je aspoň 3, součet všech stupňů je $2h$. Tedy

$$2h\geq 3(m+1).$$

Dosadíme-li do této nerovnosti $h=m+n-1$ z Eulerovy věty, dostaneme

$$2m+2n-2\geq 3(m+1),$$

tedy po úpravě $m\leq 2n-5$. Obdobně po dosazení $m=h-n+1$ je

$$2h\geq 3(h-n+2),$$

což dává $h\leq 3n-6$.

\(\Box\)

Pro množinu $P$ a každý bod $q\in \mathbb{R}^2\setminus P$ definujeme kruh se středem $q$ a poloměrem, který je roven vzdálenosti bodu $q$ od množiny $P$

$$C_P(q)=\{r\in \mathbb{R}^2;\ \operatorname{dist}(q,r)\leq \operatorname{dist}(q,P)\}.$$

Na kružnici ohraničující tento kruh leží vždy aspoň jeden bod množiny $P$.

Je zřejmé, že platí následující tvrzení:

Lemma 9.2

Bod $q$ leží na hraně V-diagramu, právě když na hranici $C_P(q)$ leží aspoň dva body z množiny $P$.

Bod $r$ je vrcholem V-diagramu, právě když na hranici $C_P(r)$ leží aspoň tři body z množiny $P$.

Obrázek 9.3
Obrázek 9.3
Ilustrace lemmatu 9.2.

Metoda zametací přímky

Jeden z možných algoritmů pro konstrukci V-diagramu používá metodu zametací přímky, ale v podobě, která je poněkud sofistikovanější než byly ty, se kterými jsme se prozatím setkali. Zametací přímka $l$ postupuje rovinou shora dolů a nad ní se vytváří V-diagram. Tento diagram nemůže vznikat v celé polorovině $l^+$ nad přímkou $l$, ale pouze v její části, která je sjednocením oblastí nad jistými parabolami. Jak taková oblast vypadá?

Bod $q$ leží v této oblasti, jestliže pro některý bod $p_i\in P\cap l^+$ platí

$$\operatorname{dist}(q,p_i)\leq \operatorname{dist}(q,l).$$

Pak totiž pro všechny body $p_j\in P$ z poloroviny $l^-$ pod přímkou $l$ je

$$\operatorname{dist}(q,p_i)\leq\operatorname{dist}(q,l)\leq\operatorname{dist}(q,p_j).$$

Tedy bod $q$ z této oblasti určitě neleží v oblasti Voronoia $V(p_j)$ pro $p_j\in P\cap l^-$.

Označme $\alpha(p,l)$ parabolu tvořenou body, které mají stejnou vzdálenost od bodu $p$ i přímky $l$. Bod $p$ je ohniskem a přímka $l$ řídící přímkou této paraboly. Nechť $\alpha^+(p,l)$ je množina bodů ležících nad parabolou $\alpha(p,l)$ nebo na ní. Oblast nad zametací přímkou $l$, kde bude V-diagram vytvořen, je sjednocení

$$\bigcup_{p_i\in P\cap l^+}\alpha^+(p_i,l).$$

Hranici této oblasti, tvořenou oblouky parabol, budeme nazývat plážová linie (beach line).

Strom $\mathcal T$ bude v tomto případě vyvážený binární strom, který ve svých listech udržuje pořadí oblouků parabol v plážové linii. Zlomové body plážové linie, tj. body kde se setkávají dva po sobě jdoucí oblouky plážové linie příslušné parabolám s ohnisky $p_1$ a $p_2$, mají stejnou vzdálenost od $p_1$ i $p_2$, neboť

$$\operatorname{dist}(q,p_1)=\operatorname{dist}(q,l)=\operatorname{dist}(q,p_2),$$

a leží tedy na hranách diagramu Voronoia. Tyto zlomové body jsou vnitřními uzly stromu $\mathcal T$. Jejich popis tvaru $\langle p_1\,:\,p_2\rangle$ značí, že se jedná o průsečík parabol s ohnisky $p_1$ a $p_2$ (obě mají řídící přímku $l$) a to ten z průsečíků (obvykle jsou dva), který má oblouk paraboly s ohniskem $p_1$ zleva a oblouk paraboly s ohniskem $p_2$ zprava. Pomocí takového stromu lze pro bod $p$ na zametací přímce $l$ vyhledat oblouk plážové linie, který leží nad bodem $p$. Následující obrázek ukazuje příklad plážové linie a příslušného stromu.

Obrázek 9.4
Obrázek 9.4
Plážová linie
Obrázek 9.5
Obrázek 9.5
Vyvážený binární strom příslušný předchozí plážové linii.

Události

Zásadní význam mají odpovědi na následující dvě otázky:

Na tyto otázky odpovídají následující dvě lemmata. Jejich důkazy jsou technické, proto je nebudeme provádět a raději se omezíme na jejich ilustraci pomocí obrázků.

Lemma 9.3

Nové oblouky v plážové linii vznikají, právě když zametací přímka prochází některým bodem množiny $P$.

Takový bod nazýváme místní událost, anglicky site event. Situaci ilustruje následující obrázek. Zametací přímka $l$ prochází bodem $p\in P$. Nechť $\alpha$ je oblouk plážové linie nad bodem $p$, který je částí paraboly s ohniskem v bodě $p_1\in P$. Po průchodu zametací přímky bodem $p$ se oblouk $\alpha$ rozdělí na dva oblouky $\alpha_1$ a $\alpha_2$ a mezi nimi vznikne v plážové linii nový oblouk $\beta$, který je částí paraboly s ohniskem v bodě $p$.

V průniku dvojice oblouků $\alpha_1$ a $\beta$ a dvojice $\beta$ a $\alpha_2$ leží body hrany diagramu Voronoia, která odděluje oblast $V(p_1)$ od oblasti $V(p)$. Této změně v plážové linii odpovídá následující změna ve stromě $\mathcal T$:

Obrázek 9.7
Obrázek 9.7
Změna ve stromě $\mathcal T$

Vnitřní uzly $\langle p_1\, :\, p_2\rangle$, $\langle p_1\,:\,p\rangle$ a $\langle p\,:\,p_1\rangle$ stromu $\mathcal T$ popisují zlomové body plážové linie pomocí ohnisek příslušných parabol v pořadí zleva doprava tak, jak jsme to popsali výše. Po změně ve stromu $\mathcal T$ je potřeba provést jeho vyvážení.

Nyní uvažujme tři po sobě jdoucí oblouky $\alpha$, $\beta$ a $\gamma$ plážové linie. Nechť $p_1$, $p_2$ a $p_3$ jsou postupně ohniska jim příslušejících parabol. Uvažujme kružnici opsanou těmto třem bodům. Nazvěme $s$ její střed a označme $q$ nejspodnější bod této kružnice, tj. bod s minimální souřadnicí $y$. Jestliže bod $q$ leží pod zametací přímkou $l$, nazýváme jej kruhovou událostí pro oblouk $\beta$ (anglicky circle event).

Při průchodu přímky $l$ tímto bodem má střed $s$ kružnice opsané bodům $p_1$, $p_2$, $p_3$ a $q$ od těchto bodů stejnou vzdálenost jako od přímky $l$, leží tedy na všech třech obloucích $\alpha$, $\beta$ a $\gamma$ plážové linie. Při posunu zametací přímky $l$ pod bod $q$, oblouk $\beta$ mizí z plážové linie. Střed $s$ kružnice opsané je vrcholem diagramu Voronoia, do kterého přicházejí hrany oddělující oblast $V(p_1)$ od $V(p_2)$, oblast $V(p_2)$ od $V(p_3)$ a oblast $V(p_1)$ od $V(p_3)$.

Obrázek 9.9
Obrázek 9.9
Změna ve stromě $\mathcal T $ při průchodu kruhovou událostí
Lemma 9.4

Oblouk paraboly mizí z plážové linie, právě když zametací přímka prochází přes jeho kruhovou událost.

Algoritmus

Začátek algoritmu odpovídá poloze zametací přímky $l$ nad všemi body. Strom $\mathcal T$ je prázdný a fronta událostí $\mathcal Q$ obsahuje všechny body z množiny $P$ (místní události) uspořádané lexikograficky shora dolů a zleva doprava. Při průchodu zametací přímky událostí provádíme v závislosti na tom, zda se jedná o místní nebo kruhovou událost, následující akce:

Tento postup provádíme tak dlouho, až je fronta $\mathcal Q$ prázdná. Strom $\mathcal T$ však prázdný nebude. Jeho aktuální vnitřní uzly určují hrany V-diagramu, které jsou polopřímkami (případně přímkami). To umožní najít pravoúhelník $R$, ve kterém leží všechny body množiny $P$, všechny vrcholy V-diagramu a zlomové body odpovídající těmto zbývajícím vnitřním uzlům stromu $\mathcal T$.

Průběh algoritmu v nejjednodušší možné situaci množiny tří bodů je demonstrován následující animací.

Přesný popis algoritmu je zachycen v následujících pseudokódech.

Algorithm VoronoiDiagram($P$)

Input. A set $P:=\{p_1,\ldots, p_n\}$ of point sites in the plane.

Output. The Voronoi diagram Vor($P$) given inside a bounding box in a doubly-connected edge list $\mathcal D$.

  1. Initialize the event queue $\mathcal Q$ with all site events, initialize an empty status structure $\mathcal T$ and an empty doubly-connected edge list $\mathcal D$.
  2. while $\mathcal Q$ is not empty
  3. do Remove the event with largest $y$-coordinate from $\mathcal Q$.
  4. if the event is a site event, occurring at site $p_i$
  5. then HandleSiteEvent($p_i$)
  6. else HandleCircleEvent($\gamma$), where $\gamma$ is the leaf of $\mathcal T$ representing the arc that will disappear
  7. The internal nodes still present in $\mathcal T$ correspond to the half-infinite edges of the Voronoi diagram. Compute a bounding box that contains all vertices of the Voronoi diagram in its interior, and attach the half-infinite edges to the bounding box by updating the doubly-connected edge list appropriately.
  8. Traverse the half-edges of the doubly-connected edge list to add the cell records and the pointers to and from them.

HandleSiteEvent($p_i$)

  1. If $\mathcal T$ is empty, insert $p_i$ into it (so that $\mathcal T$ consists of a single leaf storing $p_i$) and return. Otherwise, continue with steps 2–5.
  2. Search in $\mathcal T$ for the arc $\alpha$ vertically above $p_i$. If the leaf representing $\alpha$ has a pointer to a circle event in $\mathcal Q$, then this circle event is a false alarm and it must be deleted from $\mathcal Q$.
  3. Replace the leaf of $\mathcal T$ that represents $\alpha$ with a subtree having three leaves. The middle leaf stores the new site $p_i$ and the other two leaves store the site $p_j$ that was originally stored with $\alpha$. Store the tuples $\langle p_j, p_i \rangle$ and $\langle p_i, p_j \rangle$ representing the new breakpoints at the two new internal nodes. Perform rebalancing operations on $\mathcal T$ if necessary.
  4. Create new half-edge records in the Voronoi diagram structure for the edge separating $\mathcal V (p_i)$ and $\mathcal V (p_j)$, which will be traced out by the two new breakpoints.
  5. Check the triple of consecutive arcs where the new arc for $p_i$ is the left arc to see if the breakpoints converge. If so, insert the circle event into $\mathcal Q$ and add pointers between the node in $\mathcal T$ and the node in $\mathcal Q$. Do the same for the triple where the new arc is the right arc.

HandleCircleEvent($\gamma$)

  1. Delete the leaf $\gamma$ that represents the disappearing arc $\alpha$ from $\mathcal T$. Update the tuples representing the breakpoints at the internal nodes. Perform rebalancing operations on $\mathcal T$ if necessary. Delete all circle events involving $\alpha$ from $\mathcal Q$; these can be found using the pointers from the predecessor and the successor of $\gamma$ in $\mathcal T$. (The circle event where $\alpha$ is the middle arc is currently being handled, andhas already been deleted from $\mathcal Q$).
  2. Add the center of the circle causing the event as a vertex record to the doubly-connected edge list $\mathcal D$ storing the Voronoi diagram under construction. Create two half-edge records corresponding to the new breakpoint of the beach line. Set the pointers between them appropriately. Attach the three new records to the half-edge records that end at the vertex.
  3. Check the new triple of consecutive arcs that has the former left neighbor of $\alpha$ as its middle arc to see if the two breakpoints of the triple converge. If so, insert the corresponding circle event into $\mathcal Q$ and set pointers between the new circle event in $\mathcal Q$ and the corresponding leaf of $\mathcal T$. Do the same for the triple where the former right neighbor is the middle arc.

Časová náročnost

Nechť množina $P$, pro kterou provádíme algoritmus, má $n$ bodů. Po průchodu první místní událostí ve frontě bude plážová linie tvořena jedinou parabolou. Při průchodu dalšími místními událostmi přibudou vždy maximálně dva oblouky parabol. Proto celkový počet oblouků, které se někdy vyskytnou v plážové linii nepřevýší $2n-1$. V binárním vyváženém stromě s maximálně $2n-1$ listy trvají operace přidání dvou listů a ubrání jednoho listu a případné vyvážení stromu čas

$$O(\log(2n-1))=O(\log n).$$

Protože každý oblouk v plážové linii jednou vznikne a nejvýše jednou zanikne, bude realizovaných událostí $O(2n-1)=O(n)$. V každé události provedeme konečný shora omezený počet kroků, které trvají konečný čas nebo čas maximálně $O(\log n)$. Tím jsme dokázali:

Věta 4.1

Časová náročnost algoritmu zametací přímky pro vytvoření V-diagramu množiny $n$ bodů je $O(n\log n)$.

Implementace algoritmu pomocí Matlabu je provedena v bakalářské práci Jany Tomšů O algoritmech počítačové geometrie dostupné na

https://is.muni.cz/auth/th/175345/prif_b/