CENTER FOR MACHINE PERCEPTION
Filtrace simplexových sítí CZECH TECHNICAL
Diplomová práce
UNIVERSITY
Michal Bušta
[email protected]
Diplomová práce
23. května 2002
Lze získat na ftp://cmp.felk.cvut.cz/pub/users/sara/masters/busta2002.pdf Školitel: Dr. Techn. Radim Šára Tato práce byla podpořena Ministerstvem školství České Republiky grantem MSM 210000012 a Grantovou agenturou České Republiky grantem GACR 102/01/1371. Centrum strojového vnímání, Katedra kybernetiky Fakulta elektrotechnická ČVUT Technická 2, 166 27 Praha 6 fax: (02) 2435 7385, tel: (02) 2435 7637, www: http://cmp.felk.cvut.cz
FILTRACE SIMPLEXOVÝCH SÍTÍ DIPLOMOVÁ PRÁCE Michal Bušta ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZE FAKULTA ELEKTROTECHNICKÁ KATEDRA KYBERNETIKY 2001
Prohlašuji, že jsem předloženou práci vykonal samostatně s použitím literatury, kterou v práci řádně cituji s uvedením plného odkazu na příslušný zdroj. Souhlasím s volným využitím výsledků této práce Fakultou elektrotechniky ČVUT v Praze. V Praze dne 26. ledna 2002 Michal Bušta
Poděkování Rád bych poděkoval vedoucímu práce Dr. Ing. Radimovi Šárovi za vedení práce a cenné rady k této práci. Dále bych rád poděkoval rodině a přátelům, kteří mne podporovali.
Abstract In this thesis, we introduce a complete method for smoothing of polyhedral surfaces. We introduce a new algorithm which minimizes integrated gradient of discrete curvature. The algorithm is based on representation surface by simplex mesh. This method of smoothing is scale, rotation and translation invariant and does not change genus and shape of smoothed surface. The algorithm avoids surface shrinking, and has a simple implementation. Value of discrete mean curvature on simplex mesh is given by the inverse radius of the sphere circumscribing a simplex. We show on synthetic and real data it gives good results. Keywords: Surface reconstruction, surface smoothing, surface filtration, simplex mesh, geometric modeling.
Abstrakt V této práci je předložen kompletní postup pro vyhlazování 3D povrchu. Je zde představen algoritmus filtrace založený na minimalizaci integrálu gradientu diskrétní střední křivosti pomocí reprezentace povrchu simplexovou sítí. Použitá reprezentace zaručuje obecnost použití na povrchy libovolného topologického typu a snadnou implementaci. Velikost diskrétní stření křivosti definované na simplexové síti odpovídá inverznímu poloměru opsané koule nad vrcholem simplexu. Navržený algoritmus je invariantní vůči otočení, změně měřítka, translaci a posunutí souřadného systému. Algoritmus nesmršťuje filtrovaný povrch, zachovává jeho topologický typ a tvar. Na syntetických a reálných datech bylo prokázáno, že algoritmus dává dobré výsledky. Klíčová slova: Rekonstrukce povrchů, vyhlazování povrchů, filtrace povrchů, simplexová síť, geometrické modelování.
Obsah 1 Úvod
1
2 Vyhlazování a filtrace 3D povrchu
3
2.1
Terminologie používaná v literatuře . . . . . . . . . . . . . . . . . . .
4
2.2
Formulace problému optimalizace povrchu . . . . . . . . . . . . . . .
5
2.3
Obecné požadavky na filtr . . . . . . . . . . . . . . . . . . . . . . . .
6
3 Přehled souvisejících prací
7
3.1
Laplaceovské metody vyhlazování . . . . . . . . . . . . . . . . . . . .
8
3.2
Taubinovo vyhlazování . . . . . . . . . . . . . . . . . . . . . . . . . .
9
3.3
Bilaplaceovské vyhlazování . . . . . . . . . . . . . . . . . . . . . . . . 10
3.4
Vyhlazování pomocí diskrétní střední křivosti . . . . . . . . . . . . . 10
3.5
Modifikované vyhlazování pomocí diskrétní střední křivosti . . . . . . 11
4 Simplexová síť
13
4.1
k-Simplexová síť . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.2
Reprezentace povrchu simplexovou sítí . . . . . . . . . . . . . . . . . 15
4.3
Simplexový úhel a střední křivost . . . . . . . . . . . . . . . . . . . . 16
4.4
Metrické parametry simplexové sítě . . . . . . . . . . . . . . . . . . . 18
5 Převod mezi triangulací a simplexovou sítí
21
5.1
Převod triangulace na simplexovou síť
5.2
Převod simplexové sítě na triangulaci . . . . . . . . . . . . . . . . . . 22
6 Vyhlazování simplexové sítě
. . . . . . . . . . . . . . . . . 22
23
6.1
Vyhlazování simplexové sítě . . . . . . . . . . . . . . . . . . . . . . . 23
6.2
Algoritmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
7 Implementace
30
7.1
Reprezentace dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
7.2
Rychlost algoritmu . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 i
8 Experiment 8.1 Smršťování objektu . . 8.2 Zachování tvaru . . . . 8.3 Zachování genusu . . . 8.4 Závislost na vzorkování
. . . a
. . . . . . . . . . . . . . . . . . . . . . . . filtrace hran
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
36 36 40 44 45
9 Výsledky vyhlazování na reálných datech 47 9.1 Použitá data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 9.2 Výsledky vyhlazování . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 10 Závěr
50
Přílohy: CD se zdrovými kódy programu a návodem.
ii
Kapitola 1 Úvod V posledních letech se rekonstrukce 3D objektů stává mocným nástrojem reverzního inženýrství. Model objektu je navzorkován a převeden do vhodné reprezentace pro další zpracování. Navzorkovaná data jsou získávána různými technikami, například pomocí rangefinderů nebo stereografie. Obvyklý postup je ten, že se získají snímky objektu z různých pohledů a následně se sloučí. Většina chyb způsobených měřením je odstraněna klasickými metodami teorie zpracování signálů, a to na úrovni dat ze senzorů. Při slučování dat vznikají další chyby (např. chybná kalibrace kamery). Na sloučená data už nelze použít konvenční filtry pro digitální zpracování vzhledem k tomu, že neexistuje obecná teorie, jak reprezentovat 3D data. Neupravená data nejsou vhodná zejména pro vizualizaci předmětů. Data je tedy nutné před vizualizací upravit, resp. vyhladit. Existuje několik přístupů, jak tuto problematiku řešit. Základní metodou, která již v současné době nemá takový význam, je parametrizace povrchu funkcí z = f (x, y). Nevýhodou této metody je, že tato parametrizace výrazně omezuje třídu objektů, které lze takto vyjádřit. Metoda je vhodná například pro lokální modelování zemského povrchu. Další často používanou metodou je neparametrická reprezentace dat povrchem. Vyhlazování je u tohoto postupu dáno volbou aproximace. Tento přístup stále vyžaduje časově náročnou interakci uživatele s počítačem, vzhledem k ne příliš dobrým výsledkům při zachování hran a navíc je tato úloha výpočetně náročná. V současné době mají asi největší význam metody vyhlazování, které pracují s reprezentací bodů v síti. Metody, které byly vyvinuty pracují nejčastěji s reprezentací dat v triangulované síti a jejich přehled je v sekci ”Související práce”. Tyto metody pracují s diskrétním povrchem (povrch je reprezentován body a jejich sousedností). Aby bylo možné diskrétní povrch efektivně filtrovat, je nutné zavést parametry, které budou lokálně charakterizovat povrch. Pro úlohu filtrace jsou klíčovými pojmy normálový vektor a křivost. Tyto charakteristiky jsou důležité, protože jejich variace ovlivňují negativním způsobem výsledný ”vzhled” vizualizovaného povrchu. Pro popis křivosti 1
diskrétního povrchu se používají zejména střední a gaussovská křivost[11]. Pomocí těchto charakteristik, problém filtrace může být vyjádřen jako minimalizace variací křivosti. Minimalizováním těchto variací je určena žádaná hodnota dané veličiny v bodě diskrétního povrchu. Dalším problémem je tedy volba ”ideální” pozice bodu, ve které bude mít daný parametr žádanou hodnotu. Určení ideální pozice vrcholu (a obtížnost jejího určení) je dána reprezentací povrchu. Tato práce řeší problém tak, že jako reprezentace diskrétního povrchu byla zvolena simplexová síť[11]. Díky této reprezentaci je možné efektivně kontrolovat pozici vrcholu pomocí jeho souvislých vrcholů a tím se úloha značně zjednoduší. Práce navazuje na výzkumy v této oblasti. Je zde předložen algoritmus filtrace založený na základní charakteristice navzorkovaných povrchů - diskrétní střední křivosti a to s použitím reprezentace dat pomocí simplexové sítě.
2
Kapitola 2 Vyhlazování a filtrace 3D povrchu Pokud je na počítači rekonstruován předmět, pak jsou nejprve naměřena data, která jsou zatížená chybou. Z naměřených dat je v počítači sestavena reprezentace objektu. Reprezentací objektu může být například triangulace, simplexová síť aj. Tato reprezentace je zatížena chybou obsaženou v datech.
Apriorní model
Objekt
Merení
Data
Filtrace
Parametry
Vyhlazování
Reprezentace
Obrázek 2.1: Model rekonstruovaného povrchu a jeho reprezentace 3
Reprezentovaný objekt (povrch) je před dalším zpracováním vhodné filtrovat (vyhladit), aby dobře reprezentoval původní objekt. Situace je schématicky znázorněna na obrázku 2.1. Pokud se upravuje vzhled reprezentovaného povrchu, jsou rozlišovány dva základní pojmy: vyhlazování a filtrace. Filtrací je míněna optimalizace reprezentace tak, aby odpovídala jak naměřeným datům, tak apriornímu modelu. Apriorním modelem může být například po částech spojitý hladký povrch. Vyhlazováním je míněna optimalizace reprezentace tak, aby vyhovovala apriornímu modelu a aby splňovala dodatečná heuristická kritéria jako například “zachování objemu” uzavřeného povrchu a podobně. Rozdíl mezi vyhlazováním a filtrací je schematicky znázorněn na obrázku 2.1. Do bloku označeného “Filtrace” vstupují naměřená data a aktuální reprezentace, zatímco do bloku “Vyhlazování” vstupuje pouze aktuální reprezentace povrchu. Chybu v naměřených datech (určení pozice vrcholu) je nutné odstranit zejména proto, že velikost této chyby se nemění s hustotou vzorkování rekonstruovaného předmětu. (Co je sotva patrné při řídkém vzorkování, je velice rušivé při hustém vzorkování). Při vysoce kvalitní rekonstrukci povrchu je tedy nutné před vizualizací povrch filtrovat nebo vyhlazovat. Filtrace 3D povrchu je velice specifický problém a kapitola 2 je obecným úvodem do této problematiky.
2.1
Terminologie používaná v literatuře
V této sekci je stručný úvod do terminologie užívané v literatuře o filtraci a vyhlazování povrchů. Přehled terminologie je uveden zejména proto, že terminologie v této problémové oblasti není zcela sjednocena. • Fairing: V anglické literatuře pojem “fairing surfaces” (ze slova “fair” - pěkný) nemá český ekvivalent. Tímto pojmem se rozumí proces vyhlazování (nebo také odstraňování vysokých frekvencí povrchu). Ve výukových materiálech se uvádí, že “fairing” je nadstavba dolnofrekvenční propusti známé z teorie zpracování signálů. Tento termín se používá zejména v literatuře, která se zabývá geometrickým modelováním předmětů. • Smoothing: Českým ekvivalentem je “vyhlazování”. Tento termín má podobný význam jako termín “fairing”. Používá se zejména v literatuře, která se zabývá rekonstrukcí předmětů. Na rozdíl od termínu “fairing” se také používá tak, že zastřešuje pojmy “vyhlazování” a “filtrace”. V této práci se používají české pojmy “filtrace”, “vyhlazování” a “optimalizace”. Bude-li v dalším textu hovořeno o “optimalizaci”, pak jsou tím míněny oba pojmy 4
filtrace a vyhlazování.
2.2
Formulace problému optimalizace povrchu
Při optimalizaci diskrétního povrchu je prvním krokem definice kritéria ”kvality” povrchu. Tímto kritériem může být například střední křivost, gaussovská křivost nebo změny uvedených křivostí. Kritérium kvality odpovídá “apriornímu modelu povrchu” (Obr. 2.1). Filtrací a vyhlazováním se pak snažíme toto kritérium minimalizovat. Optimalizaci povrchu lze vyjádřit jako minimalizaci funkcionálu J:
Z |kQ | + |kD |dΩ,
J=
(2.1)
Ω
kde: Ω je povrch, kQ člen funkcionálu hodnotící lokální kvalitu povrchu, kD člen funkcionálu hodnotící lokální vzdálenost od naměřených dat, | · | je absolutní hodnota. Kritérium je tedy složeno ze dvou částí, první udává lokální kvalitu povrchu a druhá udává lokální “vzdálenost” reprezentace od naměřených dat. Pokud je v rovnici (2.1) člen kD , pak se jedná o filtraci povrchu, jinak hovoříme o vyhlazování. Je důležité si uvědomit, že z tohoto důvodu nelze proces vyhlazování provádět až do úplné konvergence (tj. nalezení extrému). Například při volbě kritéria kvality tak, aby minimalizovalo změny hlavních křivostí povrchu, by povrchy s genusem nula vždy konvergovaly ke kouli. Jelikož jsou filtry pro diskrétní povrchy většinou iterační, zanedbává se část kD funkcionálu (2.1). Tato část se zanedbává vzhledem k tomu, že měření vzdálenosti mezi dvěma diskrétními povrchy je časově náročná operace a její vyhodnocování by značně zpomalilo algoritmus. Kritérium 2.1 je definováno pro spojité povrchy. Pro diskrétní reprezentaci modelu je nutné kritérium 2.1 diskretizovat.
5
2.3
Obecné požadavky na filtr
Při vyhlazování povrchu určeného pro vizualizaci jsou nejrušivějším jevem variace křivosti. Drobné odchylky od ideálního povrchu nejsou pro pozorovatele tak rušivé jako drobné variace křivosti. Cílem filtrace povrchu je tedy omezit variace střední křivosti (a tím také tzv. normálového šumu) při zachování reprezentace povrchu. Na filtr se kladou následující požadavky: • Zachování tvaru a topologického typu (genu) předmětu. Toto je jeden ze základních požadavků na filtraci. Pokud by jakýkoliv filtr tento požadavek nerespektoval, výsledný model by nereprezentoval rekonstruovaný předmět. • Invariantnost vůči měřítku a souřadnému systému. Je důležité, aby algoritmus dával ty samé výsledky jak při natočení objektu do jiného úhlu, nebo při změně měřítka. • Obecnost. Obecnost použití pro jakýkoliv typ (zejména topologický) objektu úzce souvisí s použitou reprezentací. Velká část vyvinutých filtrů pro filtraci klade omezení na to, jaká data se mohou použít, a to značně omezuje možnost použití. Je proto vhodné volit reprezentace povrchu tak, aby pokryla co největší třídu objektů. • Zachování velikosti. Při filtraci 3D dat, zejména u uzavřených povrchů, může dojít k tak zvanému ”smršťování”. Je proto žádoucí tento efekt omezit, popřípadě vyloučit. Cílem práce bylo vytvořit vyhlazovací filtr, který bude vyhovovat výše uvedeným požadavkům. Tato práce se tedy zabývá vyhlazováním 3D povrchů.
6
Kapitola 3 Přehled souvisejících prací V této kapitole je uveden přehled základních metod pro vyhlazování povrchů reprezentovaných pomocí sítí. Všechny zde uvedené metody byly navrženy a implementovány s reprezentací triangulované sítě. Vyhlazováním povrchu, který je reprezentován pomocí simplexové sítě se zabývá práce[5], kde je představena metoda rekonstrukce a vyhlazování povrchu založená na adaptaci simplexové sítě. Při rekonstrukci touto metodou se nejprve okolo povrchu vygeneruje simplexová síť a ta se poté adaptuje tak, aby co nejlépe reprezentovala daný povrch. Tato práce se zabývá vyhlazováním již rekonstruovaného povrchu. Vzhledem k tomu, že všechny uvedené metody pracují s diskrétním povrchem reprezentovaným triangulací, bude v následujících sekcích použito toto značení: Pi je pozice i-tého vrcholu triangulace. PN j(i) je pozice j-tého sousedního vrcholu i. Sousedními vrcholy jsou míněny vrcholy spojené hranou.
7
3.1
Laplaceovské metody vyhlazování
Jednou ze základních metod pro vyhlazování 3D povrchu je Laplaceovské vyhlazování [17]. Laplaceovské vyhlazování je iterační algoritmus, který v každé iteraci vyrovnává vzájemné pozice sousedních vrcholů. Nová pozice vrcholu se vypočítá jako:
1 Pi (k + 1) = n
n X
PN j(i) (k),
(3.1)
j=1
kde: Pi (k + 1) je nová pozice vrcholu vypočítaná v k-té iteraci vyhlazovacího algoritmu, Qi (k) je pozice i-tého sousedního vrcholu, n je počet sousedních vrcholů. V každé iteraci je tedy pozice vrcholu nahrazena průměrnou hodnotou pozic sousedních vrcholů. Předností tohoto algoritmu je hlavně jednoduchá implementace, přičemž dává poměrně dobré výsledky. Při tomto vyhlazování jsou rozlišovány dvě základní modifikace. První je sekvenční verze, kdy je stará pozice okamžitě nahrazena novou. Druhá verze provede náhradu starých pozic až poté, co jsou vypočítány nové pozice vrcholů v celé síti. Vlastnosti: • Omezuje vysoké frekvence. • Při jiném než konstantním vzorkování povrchu nepřirozeně deformuje povrch[11]. • Zmenšuje povrch, obvykle konverguje do jednoho bodu. Hlavní nevýhodou této metody vyhlazování povrchu je nepřirozená deformace povrchu při nekonstantním vzorkování[11] a konvergence vyhlazovaného povrchu do jednoho bodu.
8
3.2
Taubinovo vyhlazování
Taubinovo vyhlazování[15] dává při vyhlazování povrchu velice dobré výsledky. Hlavní rozdíl oproti laplaceovskému vyhlazování je posilování nízkých frekvencí povrchu. Taubin použil analogii z teorie zpracování signálů. Pokud triangulovaný povrch představuje signál, lze zapsat proces vyhlazování jako:
x0 = f (K)x,
(3.2)
∆x = −Kx,
(3.3)
kde: f (K) je přenosová funkce filtru, x je vstupní signál (triangulace), x0 je výstupní signál, K je diskrétní aproximace Laplaceova operátoru. Přenosová funkce byla Taubinem zvolena jako: f (K) = (1 − λK)(1 − µK),
(3.4)
λ > 0,
(3.5)
µ < −λ,
(3.6)
kde: λ je pozitivní reálná konstanta, µ je negativní reálná konstanta. Pokud je přenosová funkce funkce filtru zvolena jako: f (K) = (I − λK), pak jde právě o laplaceovské vyhlazování, neboť rovnice 3.2 může být pak upravena do tvaru x0 = x + λ∆x. Taubinovo vyhlazování lze pak interpretovat tak, že po jednom kroku laplaceovského vyhlazovaní je aplikován krok zpět, pro zabránění efektu zmenšování povrchu.
9
Vlastnosti: • Potlačuje vysoké frekvence. • Posiluje nízké frekvence. • Potlačuje efekt zmenšování. • Ve vyhlazeném povrchu se objevují vlny o nízkých frekvencích[11]. Ačkoliv tento algoritmus dosahuje velice dobrých výsledků, vytváří na povrchu filtrovaného objektu vlny o nízkých frekvencích, a to značně zhoršuje výsledný ”vizuální dojem”.
3.3
Bilaplaceovské vyhlazování
Bilaplaceovské vyhlazování povrchu[11] je modifikací Taubinova vyhlazování. Při volbě konstant Taubinova vyhlazování λ = −µ se vzorec pro výpočet nové pozice vrcholu redukuje na:
f (K) = (1 − λK)(1 + λK).
(3.7)
Vlastnosti: • Omezuje vysoké frekvence. • Mírně deformuje původní povrch [11].
3.4
Vyhlazování pomocí diskrétní střední křivosti
Vyhlazování založené na diskrétní střední křivosti [7] je jednou z neúčinnějších metod pro vyhlazování 3D povrchu, a to zejména pro data určená pro vizualizaci. Algoritmus má následující matematický zápis:
Pi (k + 1) = Pi (k) + λH(Pi (k))~n(Pi (k)), kde: H(Pi ) je diskrétní střední křivost v daném bodu, ~n(Pi ) je normálový vektor v bodě Pi , 10
(3.8)
Pi (k + 1) je nový bod, který nahradí starý bod triangulace Pi (k), λ je malá kladná reánlá konstanta. Vlastnosti: • Nesmšťuje povrch. • Velice dobře vyhlazuje povrch. • Produkuje nerovnoměrné rozložení vrcholů v triangulaci [11], [7]. Dalším problémem této metody je určení nomálového vektoru ve vrcholu triangulace (obvykle se průměrují normály trojúhelníků, které daný vrchol triangulace tvoří).
3.5
Modifikované vyhlazování pomocí diskrétní střední křivosti
Toto vyhlazování odstraňuje efekt nerovnoměrného rozložení vrcholů filtrovaného povrchu u metody zmíněné v sekci 3.4. Vektor medián definovaný v [11] je dán rovnicí:
1X PN j(i) − Pi ~ = nX m 1 k PN j(i) − Pi k n
(3.9)
Modifikované vyhlazování pomocí diskrétní střední křivosti je dáno rovnicemi[11]:
Pi (k + 1) = Pi (k) + λF(Pi (k)),
(3.10)
kde:
F=
~ |H|m cos(θ)
pokud cos(θ) > ε,
~ |H|m 2H~n − cos(θ) 0
kde: 11
pokud cos(θ) < −ε, pokud | cos(θ)| ≤ ε,
(3.11)
Obrázek 3.1: Výsledky filtrace triangulovaných povrchů: a) původní koule b) koule s přidaným šumem c) filtrace laplaceovským vyhlazováním d) Taubinovo vyhlazování e) vyhlazování pomocí střední křivosti ε je malá reálná konstanta, ~ je vektor medián, m ~n je normálový vektor, θ je úhel mezi normálovým vektorem a mediánem, H je diskrétní střední křivost. Metoda je založena následujícím principu. Pokud je normálový vektor ~n téměř ~ pak se hodnota vrcholu nemění. Jinak jsou ortogonální s vektorem mediánu m, vrcholy v každé iteraci posouvány ve směru mediánu tak, že rychlost posunu v normálové složce je dána diskrétní střední křivosti.
Uvedené metody tvoří pouze základní přehled metod, které se v praxi používají. Většina z uvedených metod má své modifikace, které potlačují nevhodné vlastnosti dané metody, nicméně základní charakteristiky těchto metod zůstávají stejné. Obrázek 3.1 ukazuje výsledky diskutovaných metod (Obrázek 3.1 je převzat z [11]).
12
Kapitola 4 Simplexová síť Jako reprezentace dat byla zvolena simplexová síť [5]. Simplexová síť byla vybrána vzhledem k jejím unikátním vlastnostem, zejména díky možnosti jednoduchým způsobem vyjádřit diskrétní střední křivost a normálovou složku ve vrcholu simplexové sítě. Dalším přínosem reprezentace povrchu pomocí simplexové sítě je možnost definovat pozici vrcholu simplexové sítě pomocí jeho sousedních vrcholů. Tato vlastnost umožňuje efektivně modifikovat lokální geometrii simplexové sítě. Topologická dualita simplexové sítě a triangulace umožňuje převést triangulovaný povrch na povrch reprezentovaný pomocí simplexové sítě. Převeditelnost triangulace na simplexovou síť zvyšuje použitelnost metody filtrace, neboť velká část systémů pro rekonstrukci pracuje s triangulací. V této kapitole je stručný popis simplexové sítě a jejích vlastností, které tuto reprezentaci činí vhodnou pro úlohu filtrace nebo deformace.
4.1
k-Simplexová síť
k-Simplexová síť (k SM) M je diskrétní datová struktura tvořená množinou vrcholů V (M ) a relací sousednosti N (M ). Každý vrchol simplexové sítě má konstantní počet sousedních vrcholů. Vrchol k-simplexové sítě má vždy k + 1 sousedních vrcholů. V této práci je použita 1-simplexová (1SM) a 2-simplexová síť (2SM). Množina vrcholů V (M ) je tvořena vrcholy vi :
V (M ) = {vi , i = 1, 2, .., n}. Kde: vi je název (index) daného vrcholu, n počet vrcholů simplexové sítě. 13
(4.1)
Vlastnosti vrcholu simplexové sítě: Každému vrcholu vi ∈ V (M ) je přiřazen vektor souřadnic (poloha). Poloha vrcholu simplexové sítě bude v dalším textu značena symbolem P (vi ), nebo zkráceně Pi . Relace sousednosti N (M ) udává, které vrcholy vi ∈ V (M ) jsou “sousední vrcholy” (vrcholy spojené hranou):
N (M ) = {N (vi ), i = 1, 2, ..., n},
(4.2)
kde: n je počet vrcholů simplexové sítě, N (vi ) je množina sousedních vrcholů vi .
Vlastnosti N (vi ): Relace sousednosti N (vi ) je množina sousedních vrcholů vrcholu vi . Pro k-simplexovou síť má N (vi ) vždy k + 1 prvků. Množina N (vi ) má dále následující vlastnosti:
vj ∈ N (vi ) právě když vi ∈ N (vj ),
(4.3)
vi ∈ / N (vi ),
(4.4)
jestliže vk , vl ∈ N (vi ), potom P (vk ) 6= P (vl ).
(4.5)
Relace sousednosti tedy každému vrcholu kSM přiřazuje množinu k+1 sousedních vrcholů. Na relaci sousednosti klademe dvě omezení. První omezení udává, že žádný vrchol simplexové sítě není sám sobě sousedem, druhé omezení klade podmínku na různost polohy sousedních vrcholů a vylučuje duplicity. Pro jednodušší zápis vzorců bude v dalším textu využito následujícího zápisu:
PN 1(i) = P (N1 (vi )), kde: N1 (vi ) je první soused vrcholu vi . Označení PN 1(i) je tedy souřadnice prvního sousedního vrcholu vrcholu vi .
14
(4.6)
4.2
Reprezentace povrchu simplexovou sítí
V této práci je diskrétní povrch reprezentován pomocí 2-simplexové sítě a 1-simplexové sítě. 1-simplexová síť je polygon a je vhodná pro reprezentaci křivek, 2-simplexová síť je vhodná pro reprezentaci povrchu. V této práci je povrch reprezentován pomocí 2-simplexové sítě a jeho hranice pomocí 1-simplexové sítě. Hranice spojitého povrchu odpovídá intuitivní představě bodů na samém okraji povrchu.
Obrázek 4.1: Reprezentace povrchu pomocí simplexové sítě. Červeně je znázorněna hranice simplexové sítě reprezentována pomocí 1-simplexové sítě, modře povrch reprezentovaný pomocí 2-simplexové sítě. Pro názornost je celá situace znázorněna na obrázku 4.1. Červeně je znázorněna hranice povrchu, která je reprezentována pomocí 1-simplexové sítě, modře pak povrch reprezentovaný pomocí 2-simplexové sítě.
15
4.3
Simplexový úhel a střední křivost
V této sekci budou zavedeny parametry, které popisují lokální geometrii simplexové sítě. Jsou jimi simplexový úhel a diskrétní střední křivost [5].
Definice 1 Nechť M je 2-simplexová síť, vi je vrchol simplexové sítě a (PN1 (i) , PN2 (i) , PN3 (i) ) jsou souřadnice jeho tří sousedních vrcholů. Tři sousední vrcholy vrcholu vi pak tvoří normálovou rovinu simplexu s normálovým vektorem ~ni :
~ni =
PN1 (i) × PN2 (i) + PN2 (i) × PN3 (i) + PN3 (i) × PN1 (i) , kPN1 (i) × PN2 (i) + PN2 (i) × PN3 (i) + PN3 (i) × PN1 (i) k
(4.7)
kde: × je vektorový součin, k · k je velikost vektoru. Definice 2 Nechť vi je vrchol 2SM, S1 je opsaný kruh nad třemi sousedícími vrcholy (PN1 (i) , PN2 (i) , PN3 (i) ) vrcholu vi s poloměrem ri a středem Ci . Nechť S2 je opsaná koule nad vrcholy (Pi , PN1 (i) , PN2 (i) , PN3 (i) ) s poloměrem Ri a středem v bodě Oi . Pak simplexový úhel ve vrcholu bodě Pi , ϕi ∈ h −π, πi je definován jako:
ri sign(P P sin(ϕi ) = R ni ), i N1 (i) · ~ i kO C k cos(ϕi ) = Ri i sign(Oi Ci · ~ni ). i
(4.8)
Lokální geometrie je znázorněna na obrázku 4.2. Simplexový úhel je nezávislý na pozici sousedů vrcholu Pi na kruhu S1 a na pozici vrcholu Pi na polokouli S2 . Poloměr a střed opsané koule je vypočítán řešením determinantu [2]:
x2 + y 2 + z 2 x21 + y12 + z12 x22 + y22 + z22 x23 + y32 + z32 x24 + y42 + z42
x x1 x2 x3 x4 16
y y1 y2 y3 y4
z z1 z2 z3 z4
1 1 1 1 1
= 0,
(4.9)
Pi
Pi
ϕi
PN3
S1
ri
ri
PN1
PN2 Ri
Ri S2
S2
Obrázek 4.2: Simplexový úhel Střední křivost je charakteristika, která popisuje zakřivení povrchu. Nad simplexovou sítí je definována diskrétní střední křivost, definována pomocí simplexového úhlu jako:
Hi =
sin(ϕi ) ri ,
(4.10)
kde: Hi je diskrétní střední křivost v bodě Pi , ri je poloměr kruhu opsaného nad sousedícími vrcholy Pi , ϕi je simplexový úhel v bodě Pi . Absolutní hodnota střední křivosti může být také vypočítána jako převrácená hodnota poloměru opsané koule nad vrcholem a jeho sousedními vrcholy.
17
4.4
Metrické parametry simplexové sítě
Aby bylo možné dobře popsat lokální geometrii simplexové sítě, jsou zavedeny metrické parametry [5] simplexové sítě i = (1i , 2i , 3i ), které popisují relativní polohu vrcholu vzhledem k jeho sousedům. Metrické parametry jsou definovány pomocí ortogonální projekce vrcholu simplexové sítě do jeho normálové roviny. Nechť Fi je ortogonální projekce bodu Pi do jeho normálové roviny. Situace je znázorněna na obrázku 4.3. Pak bod Fi lze zapsat jako lineární kombinaci pozic jeho sousedních vrcholů:
Fi = 1i PN1 (i) + 2i PN2 (i) + 3i PN3 (i) ,
(4.11)
1i + 2i + 3i = 1.
(4.12)
Skaláry 1i , 2i , 3i se nazývají metrické parametry. Při volbě metrických parametrů podle rovnice (4.12) je volba metrických parametrů jedinečná.
Pi . PN1
ni PN3
Fi PN2
Obrázek 4.3: Metrické parametry simplexové sítě Z obrázku 4.3 je patrné, že pokud vyjádříme pozici vrcholu Pi relativně vzhledem k jeho sousedním vrcholům, pak lze pozici vrcholu rozdělit do normálové a tečné složky. Tečná složka pozice vrcholu Fi (4.11) leží v rovině dané sousedícími vrcholy a normálová složka je daná vzdáleností bodu Pi od této roviny.
18
Pi
ϕi
θ1 θ 2 Li ri
di Ri
Obrázek 4.4: Výpočet normálové složky pozice vrcholu
19
Velikost normálové složky je odvozena pomocí obrázku 4.4: Pro úhly θ1 , θ2 platí vztahy:
− di , tan(θ1 ) = ri L i r + i tan(θ2 ) = L di . i
(4.13) (4.14)
Vzhledem k tomu, že dále platí ϕi = π − θ1 − θ2 , rovnice (4.13), (4.13) upravíme:
− tan(ϕi ) = tan(θ1 + θ2 ) =
tan(θ1 ) + tan(θ2 ) . 1 − θ1 θ2
(4.15)
Úpravou rovnic (4.13), (4.14) a (4.15) je dán vztah pro výpočet normálové složky pozice vrcholu:
Li =
2 2 q(ri − di )tan(ϕi ) , ri + a ri2 + (ri2 − d2i ) tan2 (ϕi )
|ϕi | < |ϕi | >
π 2
π 2
(4.16)
: a = 1,
(4.17)
: a = −1.
(4.18)
Nevýhodu tohoto vzorce uvedeném v práci [5] je to, že obsahuje numerickou nestabilitu pro ϕi rovno π2 a proto je tento případ nutné vyhodnocovat zvlášť. Po zavedení metrických parametrů může být pozice vrcholu simplexové sítě vyjádřena pomocí pozice jeho sousedů, metrických parametrů a střední křivosti jako:
Fi = 1i PN1 (i) + 2i PN2 (i) + 3i PN3 (i) , Pi = Fi + Li ~ni . Kde: ri je poloměr opsaného kruhu nad vrcholy (PN1 (i) , PN2 (i) , PN3 (i) ), di je vzdálenost mezi Fi a středem opsaného kruhu, Li je velikost normálové složky definovaná 4.16, ~ni je normálový vektor vrcholu vi . 20
(4.19)
Kapitola 5 Převod mezi triangulací a simplexovou sítí Pokud je povrch, který bude filtrován reprezentován pomocí triangulace, je prvním krokem filtrace pomocí navrženého algoritmu převod triangulovaného povrchu na povrch reprezentovaný pomocí simplexové sítě. Převod dat je založen na duálním grafu simplexové sítě a triangulace, jak je znázorněno na obrázku 5.1.
Obrázek 5.1: Duální graf simplexové sítě a triangulace
21
5.1
Převod triangulace na simplexovou síť
Algoritmus převodu se dá rozdělit do třech kroků: 1. V prvním kroku se pro každý trojúhelník triangulace vygeneruje vrchol simplexové sítě, který je umístěn do středu trojúhelníku. (Těžiště trojúhelníku bylo zvoleno, aby implementace byla rychlá a převedený povrch dobře reprezentoval původní.[5]) 2. V druhém kroku se spojí vrcholy simplexové sítě tak, že pokud dva trojúhelníky sdílí hranu, pak jsou odpovídající vrcholy simplexové sítě spojeny. Pokud trojúhelník nějakou hranu nesdílí, pak je vygenerován nový vrchol simplexové sítě, který se umístí do středu dané hrany. 3. Ve třetím kroku se spojí sousedící vrcholy vygenerované v kroku 2. Už převodem triangulace na simplexovou síť se povrch mírně vyhladí. Převod dat také omezuje použití navržené metody filtrace na povrchy s dostatečnou hustotou dat.
5.2
Převod simplexové sítě na triangulaci
Převod simplexové sítě na triangulaci je úloha výpočetně náročnější. Základem algoritmu je nalezení všech stěn simplexové sítě. Stěnou simplexové sítě je míněna elementární smyčka v simplexové síti [5] (ekvivalentem českého slova “stěna” je anglické slovo “face”, termín se používá v topologických grafech). V práci byly implementovány dva způsoby převodu triangulace na simplexovou síť. První způsob převodu spočívá v tom, že do středu každé stěny je umístěn vrchol (vrchol triangulace) a poté jsou vrcholy spojeny hranou s vrcholy vygenerovanými v sousedních stěnách. Druhý způsob převodu spočívá v tom, že jednotlivé stěny simplexové sítě jsou triangulovány. Vrcholy simplexové sítě jsou pak i vrcholy triangulace. Oba způsoby převodu simplexové sítě jsou vhodné pro zaoblené povrchy. Obecně je problém převodu simplexové sítě na triangulaci další optimalizační problém, který přesahuje rámec této práce. Domnívám se, že správným způsobem převodu simplexové sítě je způsob podobný prvnímu s tím rozdílem, že poloha nově vygenerovaného vrcholu neleží ve středu stěny. Jeho poloha musí být taková, aby respektovala zakřivení povrchu v jeho okolí.
22
Kapitola 6 Vyhlazování simplexové sítě Činnost navrženého algoritmu spočívá v minimalizaci variací střední křivosti v bodech simplexové sítě. V této kapitole je vysvětlen princip, popis algoritmu filtrace a popis jeho vlastností.
6.1
Vyhlazování simplexové sítě
Předpokládejme, že povrch Ω je C 3 spojitý (parametrizace povrchu má spojité parciální derivace do řádu 3). Jako kritérium pro vyhlazování simplexové sítě bylo zvoleno kritérium:
Z EH (S) =
(∇H)2 dΩ,
(6.1)
Ω
kde: Ω je povrch, H je střední křivost (zde definovaná na spojitém povrchu), ∇ je derivace skalární funkce definovaná jako: ∂H ∂u , ∆H = ∂H ∂v
kde: u, v je parametrizace spojitého povrchu [13]. 23
(6.2)
Toto kritérium minimalizuje integrál gradientu diskrétní střední křivosti. Kritérium bylo zvoleno proto, že právě změny střední křivosti jsou jedním z nejrušivějších faktorů ovlivňující výsledný vzhled povrchu. Pro extrém funkcionálu (6.1) je nutná podmínka:
∆H = 0,
(6.3)
∆H = ∇∇H,
(6.4)
kde:
je Laplaceův operátor. Rovnici (6.4) lze řešit pomocí difuzní rovnice, která je dána vztahem:
∂H = λ∆H, ∂t
(6.5)
kde: t je čas, λ je difuzní konstanta. Difuzní rovnice konverguje k řešení rovnice (6.4) pro volbu difuzní konstanty 0 < λ < 1. Protože byla zvolena diskrétní reprezentace povrchu, je nutné rovnici (6.5) diskretizovat. Za předpokladu kruhového okolí vrcholu vi simplexové sítě N (vi ) lze aproximovat Laplacián v difuzní rovnici (6.5) jeho diskrétní verzí:
∆Hi = Hi −
1 |N (vi )|
24
X j∈N (vi )
Hj
(6.6)
Vzhledem k tomu, že na počátku filtrace není splněna podmínka kruhového okolí simplexové sítě (simplexová síť není uniformní), rovnice (6.6) je dále upravena o konstanty, které aproximují kruhové okolí:
∆Hi = Hi − X1
X (wj ) j∈N (vi )
(wj Hj ),
(6.7)
j∈N (vi )
wj =
1 , kPi PN j(i) k
(6.8)
kde: wj jsou váhové konstanty aproximující kruhové okolí, k · k je velikost vektoru. Diskretizací difuzní rovnice 6.5 a použitím vztahu 6.7 je dán vztah pro výpočet nové hodnoty diskrétní střední křivosti ve vrcholu simplexové sítě:
Hi (k + 1) = Hi (k) − λ∆Hi ,
(6.9)
kde: Hi (k + 1) je nová hodnota diskrétní střední křivosti, ∆Hi je diskrétní aproximace daná rovnicí (6.7). Porovnáním s postupem uvedeným v sekci 3.1 je patrné, že metoda laplaceovR ského vyhlazování je ve skutečnosti minimalizací funkcionálu J = min Ω ∇zdΩ při parametrizaci z = f (x, y). Další částí vyhlazování je zvýšení uniformity vzorkování simplexové sítě. Uniformita vzorkování simplexové sítě je zvyšována filtrací tečné složky pozice vrcholu simplexové sítě (4.12). Cílem filtrace tečné složky je tedy to, aby jednotlivé sousední simplexy měly stejnou vzdálenost. Filtrace tečné složky je detailně popsána v následující sekci 6.2.
25
6.2
Algoritmus
Algoritmus využívá parametrické reprezentace simplexové sítě. Vrchol simplexové sítě je tedy vyjádřen jako: Pi = 1i PN 1(i) + 2i PN 2(i) + 3i PN 2(i) + Li Ni . Vlastní filtr se skládá ze dvou částí F1 a F2 . První filtr F1 je difuzní filtr normálové složky pozice vrcholu a druhý filtr vyrovnává vzorkování simplexové sítě. Difuzní filtr F1 filtruje střední křivost a je dán rovnicí (6.9). Filtr F1 vypočítá novou hodnotu diskrétní střední křivosti v daném vrcholu simplexové sítě. Hodnotu tečné složky pozice vrcholu vypočítáme pomocí vztahu (4.16). Pro názornost je situace znázorněna na obrázku 6.1.
Pi Pi
S2
Fi S2
Obrázek 6.1: Alternativní pozice normálové složky vrcholu v závislosti na diskrétní střední křivosti (poloměru opsané koule nad simplexem).
26
Při filtraci tečné složky filtrem F2 je postupováno následujícím způsobem: 1. Je vypočítána ideální pozice tečné složky (4.11) vrcholu simplexu Fˆi . 2. Je vypočítána nová pozice vrcholu simplexu Fi (k + 1): Fi (k + 1) = Fi (k) + λ2 (Fi (k)Fˆi ),
(6.10)
kde: Fi (k + 1) je nově vypočítaná pozice tečné složky vrcholu, λ2 je rychlostní konstanta (v intervalu (0, 1 )). 3. Z nově vypočítané pozice vrcholu jsou spočítány nové hodnoty metrických parametrů dle rovnice (4.11). Volba ideální pozice tečné složky vrcholu Fˆi : Cílem filtrace tečné složky je zlepšení uniformity vzorkování simplexové sítě. Ideálním případem proto je, pokud průmět vrcholu simplexu Fi (obrázek 4.3) leží ve středu opsané kružnice nad sousedními simplexy. Experimentálně bylo prokázáno, že tento bod nelze použít, protože touto volbou se simplexová síť značně deformuje. Při volbě Fˆi do těžiště trojúhelníku daného sousedními vrcholy je při filtraci dosahováno dobrých výsledků, nicméně tato volba způsobuje smršťovaní “stěn” simplexové sítě, která jsou tvořena malým počtem vrcholů. Stěnou simplexové sítě je míněna elementární smyčka vrcholů simplexové sítě [5] (v [5] je stěna simplexové sítě označena výrazem “face”). Z uvedených důvodů byl zvolen kompromisní bod a to střed vepsané kružnice do trojúhelníka, daného sousedními vrcholy vrcholu simplexové sítě. Při této volbě dochází k menšímu smršťování stěn simplexové sítě. Tomuto jevu lze zamezit úpravou triangulované sítě na vstupu. Nyní, když jsou definovány filtry pro tečnou a normálovou složku pozice vrcholu F1 a F2 , lze v každém kroku iteračního algoritmu vypočítat novou pozici vrcholu simplexové sítě tak, že pomocí nových metrických parametrů (vypočítané filtrem F2 ) je vypočítána nová pozice projekce vrcholu simplexové sítě Fi pomocí rovnice 4.11. Hodnotu nové tečné složky pozice vrcholu vypočítáme pomocí vztahu (4.16). Výsledná pozice vrcholu je pak dána rovnicí 4.19. Optimalizace tečné a normálové složky jsou navzájem nezávislé.
27
Pro názornost je schematické vyjádření celého algoritmu znázorněno na obrázku 6.2. Jednotlivé bloky ve schématu znamenají: SM - simplexová síť, SMC - kopie simplexové sítě SM, F1 - filtr F1 definovaný v rovnici 6.9, F2 - filtr F2 definovaný v rovnici 6.10. Algoritmus nejprve inicializuje simplexovou síť. Inicializací se rozumí vypočítání metrických parametrů simplexové sítě a střední křivosti ve všech bodech simplexové sítě. Poté se vytvoří její kopie. V každém kroku iterace se pak pomocí parametrů simplexové sítě vypočítají nové parametry pro její kopii. Na základě nových metrických parametrů a informace o pozici sousedních vrcholů se vypočítají nové pozice vrcholů. Z obrázku 6.2 je také patrné, že paměťová náročnost algoritmu je 2M , kde M je velikost paměti jedné simplexové sítě. Jako podmínka zastavení algoritmu je omezení maximálního počtu iterací a minimální změna diskrétní střední křivosti. Algoritmus tedy vyhlazuje dokud není dosaženo daného počtu iterací, nebo dokud je změna diskrétní střední křivosti větší, než zadaná hodnota.
28
Obrázek 6.2: Schematické znázornění algoritmu filtrace
29
Kapitola 7 Implementace Program pro filtraci 3D povrchů byl implementován jako modul pro vývojové prostředí Matlab. Jednotlivé moduly jsou napsány v programovacím jazyku C++. V následujících sekcích jsou popsány použité abstraktní datové struktury použité při implementaci a implementované funkce.
7.1
Reprezentace dat
• Pomocné datové struktury: Pomocné datové struktury použité při implementaci slouží ke snadnějšímu přístupu k jednotlivým elementům povrchu. Jednotlivé elementy povrchu jsou u triangulace trojúhelníky a vrcholy triangulace, u simplexové sítě pak vrcholy simplexové sítě a simplexy. Pole ukazatelů. Abstraktní datová struktura POLE je dynamické pole ukazatelů s přímým přístupem. K jednotlivým prvkům pole se přistupuje pomocí indexu. Obousměrný seznam DLIST - obousměrný seznam je dynamická datová struktura, která mezi svými členy obsahuje kromě datové části i vazebné členy, které ukazují na předchozí a následující prvky seznamu. Obousměrný seznam navíc obsahuje ukazatel na první prvek seznamu. Přístup k jednotlivým datovým prvkům je sekvenční a to omezuje jeho použití na seznamy s malým počtem prvků. • Reprezentace vrcholů: Vrchol triangulace i simplexové sítě je reprezentován pomocí abstraktní třídy Vertex. Tato třída obsahuje pozici vrcholu, identifikační číslo daného vrcholu, dynamický list trojúhelníků, které daný vrchol tvoří a ukazatel na simplex, jehož je daný vrchol středem. 30
• Reprezentace trojúhelníků triangulace: Trojúhelník triangulace je reprezentován pomocí abstraktní třídy Triangle. Každý trojúhelník má své jedinečné identifikační číslo, ukazatele na vrcholy, které daný trojúhelník tvoří, ukazatele na sousední trojúhelníky a ukazatele na vrcholy simplexové sítě, které reprezentují daný trojúhelník. • Reprezentace simplexů: Jeden simplex simplexové sítě je reprezentován pomocí abstraktní třídy SM. Každý simplex obsahuje ukazatel na vrchol daného simplexu, ukazatele na sousední simplexy, normálový vektor simplexu, ukazatel na opsanou kouli nad simplexem, ukazatel na opsaný kruh simplexu, ortogonální projekci středu simplexu do roviny simplexu a metrické parametry. • Reprezentace triangulace: Triangulace je definována množinou vrcholů a množinou trojúhelníků, které tvoří triangulaci. V programu je triangulace reprezentována dvěma dynamickými poli. Pole vrcholů je dynamické pole třídy Vertex a pole triangulace je reprezentováno dynamickým polem třídy Triangle. Pro přístup k triangulaci slouží třída TriVert. • Reprezentace simplexové sítě: Simplexová síť je definována pomocí vrcholů simplexové sítě a simplexů. V programu je simplexová síť reprezentovaná pomocí dynamického pole vrcholů a dvou dynamických polí, které reprezentují simplexovou síť. Jedno pole reprezentuje vnitřní simplexy simplexové sítě (povrch bez hranice) a druhé hranici. Pro “přístup” k simplexové sítí se v programu používá třída SMVert: Třída SMVert obsahuje ukazatel na statické pole vrcholů, ukazatel na dvě statické pole simplexů (jedno reprezentující hranici povrchu a druhé reprezentující povrch) a počet simplexů, respektive vrcholů tvořících povrch.
7.2
Rychlost algoritmu
V této sekci je provedena experimentální analýza rychlosti implementovaného algoritmu. Vzhledem k podstatě problému je očekávána lineární složitost [10]. Cílem experimentu bylo ověření toho, že skutečná doba běhu implementovaného algoritmu má “lineární charakter”. Při experimentu byly filtrovány triangulované povrchy různých velikostí. Velikostí triangulovaného povrchu je zde míněno počet trojúhelníků, které povrch tvoří. Byla měřena doba převodu triangulované sítě na simplexovou síť, doba filtrace pro různý počet iterací a doba převodu simplexové sítě na triangulaci. Výsledky měření 31
0.3 data linear fit 0.25
t [s]
0.2
0.15
0.1
0.05
0
0
0.5
1
1.5 Velikost instance
2
2.5 4
x 10
Obrázek 7.1: Závislost doby převodu triangulace na simplexovou síť na počtu trojúhelníků triangulace. Modře jsou znázorněna data, červeně lineární interpolace dat. byly pro názornost zpracovány ve formě grafů. Algoritmus byl testován na počítači s procesorem AMD THUNDERBIRD 900 MHz. • Převod triangulace na simplexovou síť. Při tomto experimentu byla měřena doba převodu triangulace na simplexovou síť. Výsledky měření jsou zpracovány v grafu 7.1. Z grafu je patrné, že závislost má lineární charakter. Drobné odchylky od lineární aproximace závislosti (v grafu znázorněno červeně) jsou dány tím, že velikost instance triangulace je zde měřena počtem trojúhelníků dané triangulace a není uvažován počet vrcholů, které danou triangulaci tvoří. • Filtrace simplexové sítě. V tomto experimentu byla měřena doba filtrace simplexové sítě v závislosti na počtu iterací filtrace a na velikosti simplexové sítě. Velikost simplexové sítě je dána počtem vrcholů. Výsledky jsou zpracovány ve dvou grafech. V grafu 7.2 je závislost doby filtrace na počtu iterací, v grafu 7.3 je závislost doby filtrace na velikosti simplexové sítě. Z obou grafů (7.2, 7.3) je patrné, že obě závislosti jsou lineární. 32
180
160 I: 320 I: 1280 I: 5120 I: 20480 I: 7680 I: 24580
140
120
t [s]
100
80
60
40
20
0 50
100
150
200
250
300
350
400
450
500
k
Obrázek 7.2: Závislost doby filtrace simplexové sítě na počtu iterací pro simplexové sítě různé velikosti. Legenda udává velikost simplexové sítě, kde I je počet simplexů, k je počet iterací.
33
35
data linear fit
30
25
t [s]
20
15
10
5
0
−5
0
0.5
1
1.5 Velikost instance
2
2.5 4
x 10
Obrázek 7.3: Závislost doby filtrace simplexové sítě na velikosti simplexové sítě. • Převod simplexové sítě na triangulaci. V tomto experimentu byla měřena doba převodu simplexové sítě na triangulaci v závislosti na velikosti simplexové sítě. Výsledky jsou zpracovány v grafu 7.4. Z grafu 7.4 je patrný lineární charakter závislosti. To, že závislost není zcela lineární, je dáno tím, že tento převod je závislý na navzorkování simplexové sítě. Závislost je patrná z algoritmu převodu uvedeném v sekci 5.2.
34
1.5
data linear fit
t [s]
1
0.5
0
−0.5
0
0.5
1
1.5 Velikost instance
2
2.5 4
x 10
Obrázek 7.4: Závislost doby převodu simplexové sítě na triangulaci na velikosti simplexové sítě. Modře jsou znázorněna data, červene lineárni interpolace.
35
Kapitola 8 Experiment Pro ověření funkčnosti filtru byla použita série experimentů, která demonstruje vlastnosti algoritmu. V experimentech je ověřeno, že při vyhlazování objektu nedochází ke smršťování, předmět zachovává tvar a genus. V experimentech byl použit implementovaný algoritmus vyhlazování s konstantami λ1 = 0, 2 a λ2 = 0, 01.
8.1
Smršťování objektu
Pro ověření toho, že filtrační algoritmus nesmršťuje filtrovaný objekt, byl sestaven následující jednoduchý experiment. Byla vygenerována koule a na ní byl testován algoritmus filtrace. Při generování koule se postupovalo následujícím způsobem. Jako základ byl zvolen rovnostranný čtyřstěn se středem v počátku souřadného systému, který byl iteračně dělen. V jedné iteraci se každý trojúhelník rozdělil na čtyři rovnostranné trojúhelníky. Nově vzniklé vrcholy se vynásobily konstantou tak, aby jejich metrická vzdálenost od počátku souřadného systému byla rovna poloměru koule. Hodnota poloměru byla zvolena 2. Poté se takto vygenerovaná triangulovaná koule převedla do simplexové sítě. Pro pro číselné vyjádření odchylky vrcholů od jejich ideální polohy bylo použito následujících charakteristik: Střední hodnota:
1 x¯ = n
n X i=1
kde: x¯ je střední hodnota měřené veličiny, 36
xi ,
(8.1)
xi je i-tá naměřená hodnota měřené veličiny. Střední kvadratická chyba aritmetického průměru:
s¯ =
v n uX u (¯ x − xi )2 u t i=1
n(n − 1)
,
(8.2)
kde: s¯ je střední kvadratická chyba aritmetického průměru, x¯ je střední hodnota měřené veličiny. V tomto experimentu je měřena střední hodnota filtrované koule, a kvadratická chyba této hodnoty. Výsledky toho experimentu jsou uvedeny v tabulce 8.1 a na obrázcích 8.1 a 8.2.
Obrázek 8.1: Efekt smršťování - modře původní koule, žlutě koule po 500 iteracích Z obrázku 8.1 je patrné, že se původní koule s koulí po 500 iteracích navzájem překrývají. Pro zobrazování 3D objektů byl použit volně šiřitelný software pro Unixové systémy “Geomview”. Na obrázku 8.1 je vidět, jak tento software řeší situaci, kdy se zobrazí dva téměř stejné objekty různé barvy do stejné scény. V grafu 8.2 je možné odečíst, o kolik se zmenšila hodnota poloměru koule převodem z triangulace do simplexové sítě. Původní poloměr triangulované koule byl 2. 37
2
1.999
R
1.998
1.997
1.996
1.995
0
500
1000
1500
2000
2500 k
3000
3500
4000
4500
5000
Obrázek 8.2: Závislost střední hodnoty poloměru koule na počtu iterací. R je střední hodnota poloměru koule, k je počet iterací. Z grafu 8.2 je patrné, že efekt smršťování je potlačen na minimum. Smršťování je dáno zaokrouhlováním ve výpočtech a ne příliš robustní implementací některých geometrických úloh (například výpočet opsané koule nad simplexem).
38
n 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2200 2400 2500 2600 2800 3000 3200 3400 3500 3600 3800 4000 4200 4400 4500 4600 4800 5000
¯ R 1.998082 1.998019 1.997957 1.997896 1.997834 1.997772 1.997710 1.997649 1.997587 1.997525 1.997463 1.997402 1.997340 1.997278 1.997216 1.997155 1.997093 1.997031 1.996970 1.996908 1.996846 1.996723 1.996599 1.996538 1.996476 1.996353 1.996229 1.996106 1.995983 1.995921 1.995860 1.995737 1.995613 1.995490 1.995367 1.995306 1.995244 1.995121 1.994998
s¯ 3.791463e-06 3.741337e-06 3.897696e-06 4.071892e-06 4.252752e-06 4.436904e-06 4.622923e-06 4.810068e-06 4.997908e-06 5.186176e-06 5.374695e-06 5.563344e-06 5.752041e-06 5.940723e-06 6.129344e-06 6.317871e-06 6.506274e-06 6.694534e-06 6.882633e-06 7.070558e-06 7.258295e-06 7.633176e-06 8.007217e-06 8.193909e-06 8.380377e-06 8.752629e-06 9.123951e-06 9.494331e-06 9.863760e-06 1.004812e-05 1.023223e-05 1.059975e-05 1.096630e-05 1.133190e-05 1.169655e-05 1.187852e-05 1.206025e-05 1.242300e-05 1.278482e-05
¯ ref − R ¯ R 0 -6.252901e-05 -1.243135e-04 -1.860569e-04 -2.477902e-04 -3.095258e-04 -3.712667e-04 -4.330122e-04 -4.947605e-04 -5.565097e-04 -6.182582e-04 -6.800045e-04 -7.417477e-04 -8.034870e-04 -8.652217e-04 -9.269514e-04 -9.886757e-04 -1.050394e-03 -1.112107e-03 -1.173814e-03 -1.235515e-03 -1.358898e-03 -1.482257e-03 -1.543926e-03 -1.605590e-03 -1.728898e-03 -1.852181e-03 -1.975439e-03 -2.098672e-03 -2.160280e-03 -2.221881e-03 -2.345065e-03 -2.468224e-03 -2.591360e-03 -2.714471e-03 -2.776018e-03 -2.837559e-03 -2.960622e-03 -3.083663e-03
Tabulka 8.1: Závislost poloměru nezašuměné filtrované koule na počtu iterací
39
8.2
Zachování tvaru
Tento experiment byl sestaven tak, aby výsledek prokázal schopnost vyhladit povrch objektu a zachovat přitom jeho tvar. Algoritmus byl vyzkoušen na kouli s jednou polokoulí hladkou a druhou zatíženou normálovým náhodným šumem. Koule byla generována stejně jako v předchozím experimentu s tím rozdílem, že při poslední iteraci generování koule byl vrchol násoben hodnotou: g(τ − 0.5), pokud hodnota jeho x-ové souřadnice byla záporná. Hodnota τ je náhodná hodnota rovnoměrného rozdělení N(1, 0) (v Matlabu příkaz rand) a g je velikost přidávaného šumu. Poloměr koule byl zvolen 2 a konstanta g udávající velikost šumu byla zvolena 0,1.
N=0
N = 10
N = 100
N = 2000
Obrázek 8.3: Filtrace koule zatížené normálovým šumem - N je počet iterací Z obrázků 8.3 a 8.4 je patrné, že filtrací koule, která je zašuměná normálovým šumem vznikne opět koule, při dostatečném počtu iterací algoritmu dokonce koule téměř shodná s původní nezašuměnou koulí. Na filtrované kouli navíc nejsou známky deformace. 40
Obrázek 8.4: Rozdíl před filtrací a po filtraci koule s polovinou zatíženou normálovým šumem. Žlutě nefiltrovaná koule, modře koule po 1000 iteracích 2
R
1.998
1.996
1.994
0
500
1000
1500
2000
2500 k
3000
3500
4000
4500
5000
Obrázek 8.5: Závislost střední hodnoty poloměru zašuměné koule na počtu iterací. R je střední hodnota poloměru koule, k je počet iterací. 41
−4
7
x 10
6
s
5
4
3
2
1
0
500
1000
1500
2000
2500 k
3000
3500
4000
4500
5000
Obrázek 8.6: Závislost střední kvadratické chyby střední hodnoty poloměru zašuměné koule na počtu iterací. Na obrázku 8.4 jsou vidět dvě koule. Zašuměná nefiltrovaná koule je zobrazena žlutě. Filtrovaná koule je zobrazena modře. Aby byly vidět rozdíly mezi oběma koulemi, filtrovaná koule je zobrazována transparentně (koule je částečně průhledná). V grafu 8.5 je vidět, podobně jako v předchozím experimentu, že se střední hodnota poloměru zmenšuje. Zmenšování je však v řádu 10−3 a lze ho zanedbat. V grafu 8.6 je vidět, že střední kvadratická chyba s rostoucím počtem iterací klesá. Je také vidět, že velké lokální změny střední křivosti jsou vyhlazeny velice rychle.
42
n 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2200 2400 2500 2600 2800 3000 3200 3400 3500 3600 3800 3900 4000 4200 4400 4500 4600 4800 5000
¯ R 2.000582 1.999324 1.998953 1.998731 1.998576 1.998456 1.998354 1.998265 1.998182 1.998104 1.998030 1.997958 1.997888 1.997819 1.997752 1.997685 1.997619 1.997554 1.997489 1.997425 1.997361 1.997235 1.997109 1.997047 1.996984 1.996860 1.996737 1.996614 1.996492 1.996431 1.996370 1.996248 1.996188 1.996127 1.996006 1.995886 1.995825 1.995765 1.995645 1.995525
s¯ 6.500232e-04 3.365996e-04 2.844681e-04 2.569514e-04 2.391457e-04 2.264285e-04 2.167527e-04 2.090486e-04 2.026987e-04 1.973234e-04 1.926765e-04 1.885919e-04 1.849526e-04 1.816743e-04 1.786940e-04 1.759640e-04 1.734470e-04 1.711134e-04 1.689394e-04 1.669054e-04 1.649953e-04 1.614944e-04 1.583507e-04 1.568922e-04 1.555006e-04 1.528959e-04 1.504991e-04 1.482806e-04 1.462166e-04 1.452363e-04 1.442877e-04 1.424779e-04 1.416136e-04 1.407742e-04 1.391652e-04 1.376418e-04 1.369095e-04 1.361957e-04 1.348203e-04 1.335094e-04
¯ ref − R ¯ R 0 -1.257763e-03 -1.629161e-03 -1.851017e-03 -2.005754e-03 -2.126237e-03 -2.227559e-03 -2.317448e-03 -2.400029e-03 -2.477753e-03 -2.552089e-03 -2.623983e-03 -2.694044e-03 -2.762691e-03 -2.830218e-03 -2.896839e-03 -2.962712e-03 -3.027960e-03 -3.092674e-03 -3.156930e-03 -3.220786e-03 -3.347486e-03 -3.473073e-03 -3.535517e-03 -3.597758e-03 -3.721698e-03 -3.845010e-03 -3.967788e-03 -4.090104e-03 -4.151108e-03 -4.212017e-03 -4.333576e-03 -4.394235e-03 -4.454821e-03 -4.575784e-03 -4.696494e-03 -4.756763e-03 -4.816977e-03 -4.937251e-03 -5.057335e-03
Tabulka 8.2: Závislost poloměru zašuměné filtrované koule na počtu iterací
43
8.3
Zachování genusu
Expreriment byl sestaven stejně jako v předchozí sekci s tím rozdílem, že vyhlazovaným objektem byl torus. Výsledek filtrace je vidět na obrázku 8.7.
N=0
N = 100
N = 500
N = 1000
Obrázek 8.7: Filtrace torusu zatížený normálovým šumem - N je počet iterací Je vidět, že při filtraci torusu bylo dosaženo podobných výsledků jako u filtrované koule, na vyfiltrovaném torusu jsou patrné drobné deformace, které by byly odstraněny při dostatečném počtu iterací. Podstatné je, že nejsou patrné podstatné deformace ani známky ohýbání torusu.
44
8.4
Závislost na vzorkování a filtrace hran
Pro tento experiment byla vygenerována kostka se zkosenou hranou, zatížená normálovým šumem. Vygenerovaná kostka má nerovnoměrné vzorkování povrchu. Cílem toho experimentu je zjištění závislosti vyhlazovacího algoritmu na vzorkování simplexové sítě a jak algoritmus vyhladí “hranatý” předmět. Pro zobrazování byl sestaven modul, který umožňuje barevně zobrazovat změny křivosti ve vrcholech simplexové sítě. Pro zobrazení těchto změn byla použita barevná paleta “hot” (z anglického slova horký - barvy přecházejí postupně z černé přes červenou do bílé). Simplexová síť je tedy obarvena tak, že v místech s malou lokální změnou střední křivosti je barva černá a v místech, kde je lokální změna střední křivosti velká, je simplexová síť bílá. (Bílá barva je pouze v místě maximální hodnoty).
N=0
N = 4000
N=0
N = 4000
Obrázek 8.8: Filtrace kostky s nerovnoměrným vzokováním - N je počet iterací 45
Výsledky filtrace jsou znázorněny na obrázku 8.8. Na horních obrázcích je patrné různé zaoblení jednotlivých rohů filtrované kostky. Toto různé zaoblení je dáno vzorkováním původní triangulace respektive tím, jakým způsobem byly spojovány jednotlivé vrcholy do trojúhelníků. Tuto závislost lze odstranit tak, že vstupní triangulaci před filtrací upravíme. Úpravou je zde míněno rozbití velkých trojúhelníků na malé tak, aby všechny trojúhelníky měly zhruba stejný objem a síť je poté retriangulována. Tento krok předzpracování je poměrně jednoduché implementovat, ale je mimo rozsah této práce.
46
Kapitola 9 Výsledky vyhlazování na reálných datech Aby bylo možné plně ukázat výsledky vyhlazování, a schopnost algoritmu pracovat s reálnými daty, byl otestován na datech získávaných pomocí stereografie[3]. Data byla nejprve převedena z triangulace do simplexové sítě a poté vyhlazována.
9.1
Použitá data
Jako data byly použity výsledky z projektu stereo vidění[3] katedry kybernetiky ČVUT [4]. Byly vyhlazovány dvě lidské hlavy, každá získaná jinou metodou rekonstrukce. První hlava byla označena pracovním názvem ”BRUNT” a druhá ”KRUSTA”. Hlava BRUNT byla složena z 5392 vrcholů a 10101 trojúhelníků. Vrcholy jsou rozloženy téměř rovnoměrně po celém povrchu hlavy, a hlava netvoří uzavřený povrch. Hlava KRUSTA byla složena z vrcholů 12451 a 24580 trojúhelníků. Vrcholy jsou rozloženy poměrně rovnoměrně po celém povrchu, ale v některých místech jsou v triangulaci velké trojúhelníky. Tato hlava také netvoří uzavřený povrch.
9.2
Výsledky vyhlazování
Hlavy BRUNT a KRUSTA byly nejprve převedeny do simplexové sítě algoritmem popsaným v sekci 5.1. Poté byly obě hlavy vyhlazovány s parametry algoritmu λ1 = 0.2 a λ2 = 0.01 50 iteracemi. Výsledky vyhlazování jsou na obrázcích 9.1 a 9.2.
47
Obrázek 9.1: Vyhlazování hlavy BRUNT 1) triangulovaná hlava 2) hlava převedená do simplexové sítě 3) hlava po 50 iteracích Z obrázků je patrné, že při vyhlazování obou hlav bylo dosaženo dobrých výsledků (filtrované hlavy se zdají být hladší). Na obrázku 9.2 je vidět (v okolí nosu hlavy), že v místech, kde jsou mezi simplexy větší vzdálenosti (v původní triangulaci byl v tomto místě trojúhelník s větším obsahem než ve zbytku triangulace) se povrch mírně “vyboulil”. Při vyhlazování navrženou metodou je proto vhodné věnovat pozornost předzpracování triangulované sítě před převodem na simplexovou síť.
48
Obrázek 9.2: Vyhlazování hlavy KRUSTA 1) triangulovaná hlava 2) hlava převedená do simplexové sítě 3) hlava po 50 iteracích
49
Kapitola 10 Závěr V práci byl navržen a implementován algoritmus pro vyhlazování diskrétních povrchů. Algoritmus má dobré vyhlazovací schopnosti, které jsou dány hodnotícím kritériem kvality povrchu (6.1), založeným na diskrétní střední křivosti. Díky použité reprezentaci povrchu pomocí simplexové sítě je možné efektivně vyjádřit diskrétní střední křivost (4.10) a normálový vektor (4.7). Velikost diskrétní střední křivosti definované na simplexové síti odpovídá poloměru opsané koule nad simplexem a tím dává názornou představu o lokálním zakřivení povrchu. Navržený algoritmus, jak bylo experimentálně ověřeno v sekci 8.1, smršťuje povrch minimálním způsobem a lze tento nežádoucí efekt vyhlazování zanedbat. V experimentu v sekci 8.2 byly ověřeny vyhlazovací schopnosti algoritmu. Dále bylo ověřeno, že implementovaný algoritmus zachovává tvar vyhlazovaného povrchu a nedeformuje ho. Algoritmus velice rychle vyhlazuje velké lokální změny křivosti povrchu. V sekci 8.3 bylo experimentálně ověřeno, že algoritmus zachovává genus vyhlazovaného povrchu a nedeformuje ho. Výsledky vyhlazování pomocí navrženého algoritmu jsou mírně závislé na simplexové síti (na její lokální topologii). Pro odstranění této závislosti je možné před vyhlazováním předzpracovat vstupní data (triangulaci nebo simplexovou síť). Navržený algoritmus filtrace byl také otestován na reálných datech čímž byla ověřena jeho použitelnost v praxi. ”Slabé místo” navrženého algoritmu je převod dat ze simplexové sítě na povrch reprezentovaný pomocí triangulace. Problém tohoto převodu je diskutován v sekci 5.2. Problém toho převodu je v této práci zmíněn pouze okrajově a přesahuje rámec této práce. Pro dosažení lepších výsledků (zejména pro po částech rovné povrchy) by bylo třeba se tímto problémem zabývat detailněji. Pro dosažení dobrých výsledků vyhlazování povrchu touto metodou je vhodné upravit vstupní data (diskrétní povrch) tak, aby tento povrch měl rovnoměrné vzor50
kování. Tato metoda filtrace 3D povrchů je použitelná na většinu typů povrchů. Díky tomu, že odstraňuje variace střední křivosti, je vhodná zejména při rekonstrukci geometrických dat ve velkém rozlišení, která jsou určena pro vizualizaci.
51
Literatura [1] Amenta, Bern, and Eppstein. Optimal point placement for mesh smoothing. In SODA: ACM-SIAM Symposium on Discrete Algorithms (A Conference on Theoretical and Experimental Analysis of Discrete Algorithms), pages 528–537, 1997. [2] W. H. Beyer. CRC Standard Mathematical Tables, p. 227. 28th ed. Boca Raton, FL: CRC Press, 3-rd edition, 1987. [3] Czech Technical University Prague Centre for Machine Perception. Stereoscopic Vision. World Wide Web, http://cmp.felk.cvut.cz/demos/Stereo.html, 1999. [4] Czech Technical University Prague Centre for Machine Perception. World Wide Web, http://cmp.felk.cvut.cz/, 2001. [5] Hervé Delingette. Simplex meshes: a general representation for 3D shape reconstruction. In Conf. on Computer Vision and Pattern Recognition (CVPR ’94), pages 856–859, June 1994. [6] Herve Delingette. General object reconstruction based on simplex meshes. Technical Report RR-3111, INRIA Sophia-Antipolis, 1997. [7] Mathieu Desbrun, Mark Meyer, Peter Schröder, and Alan H. Barr. Implicit fairing of irregular meshes using diffusion and curvature flow. Computer Graphics, 33(Annual Conference Series):317–324, 1999. [8] S. Karbacher and G. ausler. A new approach for modeling and smoothing of scattered 3d data. In R. N. Ellson and J. H. Nurre, editors, Three-Dimensional Image Capture and Applications, 3313 of SPIE Proceedings, pages 168–177, Bellingham, Washington, 1998. The International Society for Optical Engineering., 1998. [9] D. T. Lee. Computational geometry. ACM Computing Surveys, 28(1):27–31, 1996. 52
[10] Vladimír Mařík, Štěpánková, Jiří Lažanský, and kolektiv. Umělá inteligence 3, pages 262–310. ACADEMIA Praha, 2001. [11] Yutaka Ohtake, Alexander G. Belyaev, and Ilia A. Bogaevski. Polyhedral surface smoothing with simultaneous mesh regularization. In Geometric Modeling and Processing (GMP), pages 229–237, 2000. [12] Stanislav Racek. Objektově orientované programování v C++. Nakladatelství KOPP, 1995. [13] Martin Raussen. Elementary Differential Geometry: ves and Surfaces, Lecture Notes. World Wide http://www.math.auc.dk/˜raussen/mrteachgg.html, 1999.
CurWeb,
[14] Karel Rektorys. Přehled užité matematiky. SNTL, 3-rd edition, 1973. [15] Gabriel Taubin. A signal processing approach to fair surface design. Computer Graphics, 29(Annual Conference Series):351–358, 1995. [16] Jan Paseka (Masarykova univerzita). Geometrické modelování. World Wide Web, http://www.math.muni.cz/pub/math/people/Paseka/lectures, 1998. [17] J. Vollmer, R. Mencl, , and H. Môller. Improved laplacian smoothing of noisy surface meshes. In P. Brunet and R. Scopigno, editors, Computer Graphics Forum (Eurographics ’99), volume 18(3), pages 131–138. The Eurographics Association and Blackwell Publishers, 1999.
53
