Dissertação
dissertacao.jose.borges(1).pdf
Documento PDF (2.3MB)
Documento PDF (2.3MB)
Operador Laplaciano Discreto
via Triangulação de Delaunay Intrı́nseca
José Borges dos Santos Filho
Agosto de 2008
Resumo
O objetivo desta dissertação é apresentar um análogo discreto do operador laplaciano, ou seja, um operador linear definido no conjunto das funções lineares
por partes em uma malha de triângulos que possua o máximo de propriedades
análogas ao operador laplaciano contı́nuo sobre uma superfı́cie. Em particular, mostraremos que se a malha satisfaz ao critério de Delaunay, o laplaciano
obedece a uma versão discreta do princı́pio do máximo, que possui importância
semelhante ao princı́pio do máximo na teoria das funções harmônicas. Apresentamos ainda três aplicações do laplaciano discretizado: a primeira tem como
objetivo obter parametrizações de malhas para efeito de mapeamento de textura;
a segunda consiste na suavização de malhas por meio do processo de difusão;
a terceira e última aplicação visa identificar formas e simetrias de objetos por
meio das curvas de contorno associadas às autofunções do laplaciano.
Abstract
The main goal of this work is to present a discrete analogous of the laplacian
operator, that is, a linear operator on the set of piecewise linear functions over
a triangular mesh that has similar properties to the continuous laplacian over
a surface. Particularly, we will show that if the mesh satisfies a Delaunay
criterion, the laplacian obeys a discrete version of the maximum principle, which
importance in the discrete setting is similar to the importance of the maximum
principle in the theory of harmonic functions. We also present three applications
of the discrete laplacian: the first one has as objective to get parametrizations
of meshes for texture mapping; the second one consists of mesh smoothing
by a diffusion process; the third and last application aims to identify forms
and symmetries of objects by means of the contour curves associated to the
eigenfunctions of the laplacian operator.
Conteúdo
1 Introdução
5
2 Operador Laplaciano
2.1 Operador Laplaciano Contı́nuo . . . . . . . . . . . . . . . . . . .
2.2 Propriedades do Laplaciano . . . . . . . . . . . . . . . . . . . . .
2.3 Discretização do Laplaciano Contı́nuo . . . . . . . . . . . . . . .
2.4 Propriedades do Laplaciano Discreto . . . . . . . . . . . . . . . .
2.5 Laplaciano Discreto Geométrico . . . . . . . . . . . . . . . . . . .
7
7
11
12
13
13
3 Triangulação de Delaunay Intrı́nseca
18
3.1 Triangulação de Delaunay no Plano . . . . . . . . . . . . . . . . . 18
3.2 Algoritmo de Flip de Arestas no Plano . . . . . . . . . . . . . . . 21
3.3 Extensão do Algoritmo de Flip Para Malhas no Espaço R3 . . . . 25
4 Aplicações do Laplaciano Discreto
28
4.1 Mapeamento de Textura . . . . . . . . . . . . . . . . . . . . . . . 28
4.2 Suavização de Superfı́cies . . . . . . . . . . . . . . . . . . . . . . 34
4.3 Identificação de Formas e Simetrias de Malhas . . . . . . . . . . 36
5 Conclusão
39
2
Agradecimentos
Primeiro agradeço a Deus pela oportunidade de desenvolver este trabalho. Agradeço a meus pais, José Borges dos Santos e Therezinha Bispo dos Santos. De
forma especial, agradeço a Vinicı́us Mello, meu orientador, sempre disposto
a auxiliar-me, embora muitas fossem as suas atividades. Destaco e agradeço
as valiosas contribuições do meu colega de mestrado Leonardo Carvalho na
implementação de alguns algoritmos. Agradeço também aos professores, colegas,
funcionários e a tantos outros que de forma direta ou indireta contribuı́ram para
que este trabalho chegasse ao seu final, após exaustivo esforço e dedicação.
3
Para ser grande, sê inteiro: Nada teu exagera ou exclui. Sê todo em cada coisa.
Põe quanto és no mı́nimo que fazes. Assim em cada lago a lua toda brilha,
porque alta vive.
4
Capı́tulo 1
Introdução
É consenso hoje que a Computação Gráfica é uma necessidade em quase todas as áreas do conhecimento, pois através dela é possı́vel visualizar objetos
em fase de projeto ou fora do alcance de nossa percepção visual, bem como
fora da nossa realidade tridimensional. A realidade do dia a dia consiste em
um universo regido pelas leis da fı́sica, e para compreender os fenômenos e
objetos existentes neste mundo fı́sico, inicialmente lançamos mão de modelos
matemáticos abstratos, geralmente contı́nuos, ou seja, que não são encontrados no mundo fı́sico. Em seguida, procuramos uma representação simbólica e
finita desse modelo matemático, ou seja, buscamos uma discretização do modelo matemático que permita sua implementação e visualização por meio de um
computador.
Como exemplo podemos citar o operador laplaciano definido em funções
f : Rn → R duas vezes diferenciáveis, cuja expressão para n = 3 é dada por
∆f =
∂2f
∂2f
∂2f
+
+
.
∂x2
∂y 2
∂z 2
Este operador surge de forma natural no estudo dos processos de difusão térmica,
aparecendo no segundo membro da chamada equação de difusão,
∂u
= η∆u,
∂t
onde u = u(x, y, z, t) representa o campo de temperatura; também aparece em
processos vibratórios na equação da onda,
∆f = λf
onde λ é um coeficiente associado ao meio pelo qual a onda se propaga.
No transcorrer deste trabalho será apresentado um laplaciano contı́nuo intrı́nseco a uma superfı́cie, ou seja, um laplaciano definido em funções cujo domı́nio
é uma superfı́cie, mas nosso principal objetivo será determinar um laplaciano
discretizado que se aproxime ao máximo do laplaciano contı́nuo, ou seja, um
laplaciano cujo domı́nio seja o conjunto das funções definidas em uma superfı́cie
discretizada. É possı́vel obter esta discretização de diversas maneiras e dependendo da escolha pode-se priorizar informações combinatórias ou geométricas de
uma malha. Um das maneiras de se discretizar o laplaciano contı́nuo é através de
5
diferenças finitas, só que este método limita-se ao uso de malhas de retângulos,
limitação que pode acarretar uma série de dificuldades. Apresentamos neste
trabalho um método de discretização que se destaca por priorizar informações
geométricas da malha e por possuir várias propriedades do laplaciano contı́nuo,
entre elas uma versão discreta do princı́pio do máximo, caso a triangulação da
malha satisfaça um critério de Delaunay.
Concluiremos exemplificando três das diversas aplicações do laplaciano discreto: a primeira é o mapeamento de textura, onde inicialmente determinamos
uma malha cuja triangulação é a especial triangulação de Delaunay, na qual
o laplaciano discretizado satisfaz ao princı́pio do máximo, possibilitando esta
propriedade a construção de uma parametrização homeomorfa da superfı́cie discretizada (malha) em uma região poligonal do plano. De posse desta parametrização, faz-se o mapeamento de textura do plano na malha. A segunda aplicação
são as curvas de nı́vel das autofunções do laplaciano, também chamadas de curvas de contorno, permitindo conhecer a sua geometria. Por fim, utilizaremos
o laplaciano para suavização de superfı́cies, isso porque ao aplicá-lo em uma
malha, uma pequena pertubação é dispersada em sua vizinhança, suavizando
assim regiões com altas freqüências, enquanto a forma principal da malha é
apenas levemente deformada.
Este trabalho teve como fonte principal o artigo [1] na parte teórica e os
artigos [3],[13],[4] e [5] na parte de aplicações.
6
Capı́tulo 2
Operador Laplaciano
Neste capı́tulo vamos definir o operador laplaciano em uma superfı́cie, apresentar algumas das suas propriedades e resultados e encontrar um análogo discreto
deste operador, ou seja, determinar um laplaciano definido em uma malha de
triângulos que possua o maior número de propriedades similares as do laplaciano
contı́nuo definido em uma superfı́cie.
2.1
Operador Laplaciano Contı́nuo
O operador laplaciano surge naturalmente do estudo de processos de difusão
térmica, aparecendo no segundo membro da equação do calor, também chamada
de equação de difusão. Existem diversas variações para a equação do calor, a
mais clássica é:
∂u
= η∆u,
(2.1)
∂t
onde u = u(x, y, z, t) representa o campo de temperatura, ∆u é o laplaciano de
u no R3 dado por
∂2u ∂2u ∂2u
∆u =
+ 2 + 2
∂x2
∂y
∂z
e η é o coeficiente de difusão. Essa equação descreve a maneira como o calor se
difunde em um meio homogêneo, isotrópico e que não possua fontes de calor.
Também em processos vibratórios estudados pela fı́sica o laplaciano aparece
no segundo membro da equação da onda,
∂2u
= c∆u,
∂t2
(2.2)
onde c representa coeficiente associado ao meio pelo qual a onda se propaga.
Esta equação modela as vibrações transversais de baixa amplitude de uma membrana S fina fixada em um aro com formato ∂S: se ∂S é um retângulo, estaremos estudando as vibrações de uma membrana retangular; se ∂S é um cı́rculo,
estudamos as vibrações de uma membrana circular (um tambor), e assim por
diante.
A equação da onda e a sua relação com o laplaciano foi estudado de forma
lúdica pelo fı́sico Ernest Chladni, em seu livro “Discoveries Concerning the
Theories of Music” publicado em 1787. Neste livro ele informa as observações
7
obtidas ao pôr areia sobre um prato de metal fino em vibração, notando que
a areia acumulava-se em certas zonas, formando padrões surpreendentemente
complexos (figura 2.1). Este comportamento pode ser explicado pela teoria das
Figura 2.1: Padrões formados pela areia acumulada em algumas regiões de um
prato de metal submetido a várias vibrações (experimentos de Ernest Chladni1700).
ondas estacionárias, que podem ser representadas pelas soluções da equação
∆f = λf,
ou seja, pelas autofunções do laplaciano. Os zeros destas autofunções, denominados de conjuntos nodais, correspondem aos locais onde a areia se acumula no
prato de Chladni.
Veremos nas próximas seções mais detalhes sobre estas autofunções. Por
hora, começaremos com alguns conceitos e resultados de geometria diferencial
que serão utilizados neste trabalho e por isso faz-se necessário conceituá-los.
2.1.1
Geometria Diferencial
Como vimos anteriormente, o laplaciano está definido em um espaço de funções
f ∈ C 2 (Rn ) e quando n = 3 o laplaciano é dado por
∆f =
∂2f
∂2f
∂2f
+
+
.
∂x2
∂y 2
∂z 2
Nesta seção nosso objetivo é encontrar uma expressão para o laplaciano definido
em funções f ∈ C 2 (S) onde S é uma superfı́cie regular, ou seja, queremos um
laplaciano intrı́nseco a uma superfı́cie regular. Mas o que é uma superfı́cie regular? Vejamos a definição segundo [8], nossa principal referência neste assunto:
Definição 1 (Superfı́cies Regulares). Um subconjunto S de R3 diz-se uma superfı́cie regular se, para cada p ∈ S, existirem uma vizinhança aberta V ⊂ R3
de p, um aberto U ⊆ R2 , e uma bijeção ϕ : U −→ V ∩ S com as seguintes
propriedades:
i) ϕ é de classe C ∞ ;
ii) ϕ é um homeomorfismo;
iii) para qualquer q ∈ U a matriz jacobiana Jϕ(q) tem posto dois.
8
Uma função ϕ com tais caracterı́sticas é dita uma parametrização local de
S, ou seja, uma vizinhança parametrizada de p contida em S. Denotaremos
os pontos de U por (u, v), sendo u e v parâmetros locais de S. As derivadas
parciais de ϕ denotaremos por ϕu e ϕv .
Ainda precisamos de mais dois conceitos para definir o laplaciano intrı́nseco
de funções f ∈ C 2 (S), o gradiente e o divergente de f .
Definição 2 (Gradiente). Seja S uma superfı́cie regular e f : S −→ R uma
função diferenciável. O gradiente de f é o campo de vetores tangentes com
→
valores em S, denotado por ∇f , assim definido: Dado p ∈ S e v ∈ Tp S, onde
Tp S é o plano tangente a S em p, temos:
D
E
→
→
dfp ( v ) = ∇f (p), v
→
onde dfp ( v ) é a derivada de f em p.
De acordo com a definição 2 e sabendo que dfp (ϕu ) = fu e dfp (ϕv ) = fv ,
verifica-se que o gradiente em coordenadas é dado por:
∇f =
fu G − fv F
fv E − fu F
ϕu +
ϕv
EG − F 2
EG − F 2
(2.3)
onde E, F e G são os coeficientes da primeira forma fundamental (maiores detalhes consultar [8]). Observe que se S for o plano R2 , E = G = 1 e F = 0,
sabendo que no plano R2 ϕu = (1, 0) e ϕv = (0, 1), o gradiente de f definido no
R2 é, portanto,
∇f = (fu , fv ) .
Vejamos agora qual a definição do divergente de f ∈ C 2 (S), uma vez que,
como já mencionamos anteriormente, é fundamental para o conceito do laplaciano.
Definição 3 (Divergente). Seja w um campo de vetores definido em U ⊂ S, a
→
divergência de w é a aplicação div w: U ⊂ S −→ R tal que, para cada p ∈ U ,
→
→
div w (p) é o traço da aplicação linear D w p ,
→
onde D w p é a matriz da derivada covariante de w no ponto p ∈ S. Explici→
taremos D w p com relação as bases {(1, 0), (0, 1)},{ϕu , ϕv }:
au + aΓ111 + bΓ112 av + aΓ112 + bΓ122
→
(2.4)
D wp=
bu + aΓ211 + bΓ212 bv + aΓ212 + bΓ222
sendo w = aϕu +bϕv e Γkij os sı́mbolos de Christoffel (maiores detalhes consultar
[8]). A partir da definição 3 e de 2.4 concluimos que,
div(w) = (au + aΓ111 + bΓ112 ) + (bv + aΓ212 + bΓ222 ).
(2.5)
Finalmente estamos prontos para definir o laplaciano e obter uma equação
para expressá-lo em uma superfı́cie regular. Vejamos primeiro a definição e logo
após o explicitaremos.
9
Definição 4 (Operador Laplaciano). Dado uma superfı́cie S associaremos a
ela um operador diferencial de segunda ordem, o operador laplaciano, definido
como ∆f = div(∇f ) onde f ∈ C 2 (S), ou seja, f : S → R é uma função real
definida em S duas vezes diferenciável.
A partir de 2.3 e 2.5 o operador laplaciano de uma função f ∈ C 2 (S) onde
S é uma superfı́cie regular em coordenadas é dado por:
∆f = (au + aΓ111 + bΓ112 ) + (bv + aΓ212 + bΓ222 ),
(2.6)
sendo
fv E − fu F
fu G − fv F
e b=
.
2
EG − F
EG − F 2
Observe que no caso particular do plano a = fu e b = fv , consequentemente o
laplaciano de f ∈ C 2 (R2 ) é dado por
a=
∆f = fuu + fvv .
(2.7)
Um resultado importante e que nos interessa neste trabalho é o teorema da
divergência cuja demonstração pode ser vista em [8],
Teorema 1 (Teorema da Divergência). Seja f um campo de vetores tangentes
a uma superfı́cie S, e Ω ⊆ S uma região poligonal. Então
Z
Z
f.nds,
(2.8)
div(f )dA =
Ω
γ
onde γ(s) é a fronteira de Ω e n é o vetor unitário ortogonal a γ ′ (s) e que
aponta para fora de Ω.
A importância deste teorema está no fato dele possibilitar a demonstração
do corolário seguinte também conhecido como primeira identidade de Green.
Corolário 2 (Primeira Identidade de Green). Sejam f, g : Ω ⊂ S → R funções
escalares definidas em uma região Ω ⊂ S que possuem derivadas parciais e
contı́nuas então,
Z
Z
Z
∂f
∇f.∇gdA,
(2.9)
g ds −
∆f gdA =
Ω
γ ∂n
Ω
onde γ é fronteira de Ω ⊂ S e S uma superfı́cie regular.
Demonstração. Substituindo em 2.5 w por g∇f obtemos,
div(g∇f ) = (gu a + gau + gaΓ111 + gbΓ112 ) + (gv b + gbv + gaΓ212 + gbΓ222 ) (2.10)
por 2.3 sabemos que
∇f =
fu G − fv F
fv E − fu F
ϕu +
ϕv
EG − F 2
EG − F 2
∇g =
gu G − gv F
gv E − gu F
ϕu +
ϕv .
EG − F 2
EG − F 2
e
10
Com essas duas equações é fácil mostrar que,
∇f.∇g = agu + bgv ,
(2.11)
onde
fv E − fu F
fu G − fv F
e
b=
.
2
EG − F
EG − F 2
Substituindo agora 2.11 em 2.10 e usando 2.5 obtemos
a=
div(g∇f ) = ∇g∇f + g div(∇f ),
logo
div(g∇f ) = ∇f.∇g + g∆f.
Fazendo v = g∇f em 2.8 obtemos
Z
Z
div(g∇f )dA =
g∇f.nds.
Ω
(2.12)
(2.13)
γ
Por fim, substituindo 2.12 em 2.13 obtemos
Z
Z
Z
∂f
∇f.∇gdA,
∆f gdA =
g ds −
Ω
Ω
γ ∂n
∂f
onde ∇f.n que aparece na integral de linha de 2.13 foi substituı́da por ∂n
, sendo
assim a derivada direcional de f na direção do vetor normal n, também chamada
de derivada normal de f .
A primeira identidade de Green terá fundamental importância nesse trabalho, pois ela possibilitará definir o laplaciano em um conjunto de funções
lineares por partes definidas em uma malha de triângulos, ou seja, um laplaciano discretizado.
2.2
Propriedades do Laplaciano
Na seção anterior definimos o laplaciano em uma superfı́cie,
∆u = div(∇u).
Apresentamos a seguir algumas propriedades desse operador onde o produto
interno utilizado é:
Z
hv, uiC 2 (S) =
uvdA
S
2
sendo u , v ∈ C (S).
1. Elemento neutro
∆u = 0 sempre que u for uma constante. Este resultado é óbvio pela
própria definição do laplaciano.
2. Simetria
h∆u, viC 2 (S) = hu, ∆viC 2 (S) sempre u, v ∈ C 2 (S) e se anulem em ∂S. Esta
propriedade é facilmente verificada usando 2.9.
11
3. Precisão linear
∆u = 0 Sempre que S for parte do plano euclidiano e u = ax + by + cz é a
função linear do plano. De acordo com 2.7, ∆u envolve apenas soma das
derivadas parciais de u de segunda ordem, como u é uma função linear
definida no plano temos ∆u = 0.
4. Suporte Local
Para algum par de pontos p 6= q ∆u(p) independe u(q), ou seja, alterando
o valor da função não afetará o valor local do laplaciano.
5. Princı́pio do máximo
∆u = 0 no interior de S, então u não possui um máximo (ou mı́nimo)
no interior de S, exceto quando u = cte. A verificação desta propriedade
pode ser vista em [12]
2.3
Discretização do Laplaciano Contı́nuo
Uma discretização do laplaciano pode ser obtida de maneiras diferentes, podendo conter mais informações combinatórias ou mais informações geométricas
a depender da discretização escolhida. Se considerado apenas o aspecto combinatório, informações métricas são desprezadas, como o comprimento de arestas
ou ângulos da malha. Veremos mais tarde que esse tipo de discretização do
laplaciano, em geral, não é uma boa escolha.
Definiremos, na próxima seção, o laplaciano discreto para uma malha M de
triângulos, mas antes vejamos qual a definição para uma malha de triângulos.
Definição 5 (Triangulação). Dado um conjunto V = {v1 , v2 , ......, vn } de pontos de uma superfı́cie S, uma triangulação de V é um complexo simplicial M
de ordem dois, ou seja, é um conjunto de triângulos, em que dados quaisquer
dois triângulos, uma das quatro situações ocorre: eles são disjuntos, eles têm
exatamente um vértice em comum, eles têm exatamente uma aresta em comum
ou são iguais.
Vejamos agora a definição do laplaciano em malhas de triângulos.
Definição 6 (Operador Laplaciano Discreto). Considere uma malha triangular
M de uma superfı́cie S, com um conjunto de vértices V , conjunto de arestas E
e um conjunto de faces F . Definimos o operador laplaciano discreto ∆d para
M pela ação linear de funções ui definidas nos vértices, como:
X
(∆d u)i =
wij (ui − uj )
(2.14)
j
onde i e j são ı́ndices associados aos vértices.
O fator determinante na definição 6 é justamente o peso wij . Verificaremos
nas próximas seções que wij está associado a cada aresta da malha M, e que
este valor determinará que aspecto da malha será priorizado. A próxima seção
destaca algumas propriedades do laplaciano discreto ∆d .
12
2.4
Propriedades do Laplaciano Discreto
Cada propriedade a seguir do laplaciano discreto ∆d tem sua motivação centrada
em uma propriedade do laplaciano contı́nuo, essas propriedades são verificadas
a partir das propriedades da matriz cuja as entradas são os pesos wij , sendo que
wij = 0 sempre que i e j não são vértices de uma mesma aresta de M.
1. Simetria
wij = wji . Matriz simétrica com autovalores associados a autovetores
ortogonais. Esta propriedade é verificada devido ao fato de eij e eji representarem a mesma aresta.
2. Suporte local
Mudando o valor da função uj não alterará a ação do laplaciano (∆d u)i ,
se i e j não formam uma aresta na malha. Laplaciano é um conceito local.
3. Precisão Linear
(∆d u)i = 0 em cada vértice, sempre que se tenha uma malha mergulhada
no plano e u é uma função linear deste plano.
4. Princı́pio do Máximo
Caso wij < 0, sempre que i 6= j, o princı́pio do máximo também é satisfeito. Esta propriedade será verificada em seção seguinte.
2.5
Laplaciano Discreto Geométrico
Nosso objetivo nesta seção é obter um laplaciano discretizado que considere
propriedades geométricas de uma malha qualquer M.
Como foi mostrado em 2.9 o operador laplaciano definido sobre uma superfı́cie M satisfaz a primeira identidade de Green:
Z
Z
Z
∂f
h∇f, ∇gi dV
(2.15)
g ds −
∆f gdA =
M
γ ∂n
M
quaisquer que sejam f, g ∈ C ∞ (M), sendo γ a fronteira de M. Suponhamos
M sem bordo e definindo a forma bilinear
Z
h∇f, ∇gi dV,
(2.16)
B(f, g) := −
M
temos
h∆f, gi :=
Z
∆f gdV = B(f, g),
(2.17)
M
logo, o operador ∆ é o operador linear associado à forma bilinear B, com respeito
ao produto interno usual em C ∞ (M).
Nosso objetivo é encontrar um análogo discreto do laplaciano contı́nuo, ou
seja, um operador linear ∆d definido sobre o espaço L(M) de todas as funções
lineares por partes definidas em uma malha de triângulos M no R3 , e que
satisfaça 2.15, pois essa condição garante que −∆d seja simétrico, positivo e
13
semi-definido. Como a forma bilinear B(f, g) faz sentido para quaisquer funções
f, g ∈ L(M) definiremos o laplaciano discreto geométrico como segue:
Definição 7 (Laplaciano Discreto Geométrico). Seja L(M) o conjunto de todas
as funções contı́nuas e lineares por parte definidas em uma malha M, definimos
o laplaciano discreto geométrico como o único operador linear ∆d que satisfaz,
h∆d f, gi := B(f, g)
(2.18)
para todo f, g ∈ L(M).
Tomando a base canônica de L(M), ou seja, as funções βi ∈ C 0 (M) tais
que βi é linear em cada simplexo de M e βi (vj ) = δij , para todo vj ∈ V , toda
função h ∈ L(M) pode ser escrita como
X
h=
ai βi ,
daı́
f=
X
ai βi
,
g=
X
bj βj .
(2.19)
Substituindo as equações de 2.19 no segundo membro da equação 2.18 obtemos,
Z
Z
X
X
−
h∇f, ∇gi = −
ai bj
h∇βi , ∇βj i,
M
M
vi ∈V (M) vj ∈V (M)
fazendo os pesos
wij =
Z
h∇βi , ∇βj i
M
obtemos,
−
Z
M
h∇f, ∇gi = −
X
X
ai bj wij .
(2.20)
vi ∈V (M) vj ∈V
Para continuar desenvolvendo o segundo membro de 2.18 precisamos do
seguinte lema, cuja demonstração será omitida aqui:
Lema 3. Sejam p,q e r pontos do R2 não colineares e f ∈ L(M) uma função
linear por parte, tal que f (p) = 1, f (q) = 0 e f (r) = 0. Então,
1. a direção de ∇f é perpendicular a aresta eqr , com extremidades nos vértices
q e r.
2. |∇f | é igual ao inverso da altura do triângulo com relação a aresta eqr (
ver figura 2.2).
Desenvolveremos wij em 2.20, pois o seu valor será de fundamental importância para o próximo capı́tulo desse trabalho, logo, é necessário aqui determinar sua expressão, dividiremos em dois casos:
1. caso i 6= j
Neste caso wij 6= 0 apenas quando vi e vj são vértices de uma mesma
aresta, ( ver figura 2.3).
Pelo lema 3 temos,
14
p
hprq
∇p
θr
r
θq
q
Figura 2.2: Triângulo de vértices p, q e r, com f (p) = 1, f (q) = 0 e f (r) = 0.
vi
vk
t1
θk
t2
vl
θl
vj
Figura 2.3: triângulos t1 e t2 em que a aresta eij é incidente.
|∇βi | |t1 =
1
1
, |∇βj | |t1 = j
hijk
hik
(2.21)
|∇βi | |t2 =
1
1
, |∇βj | |t2 = j ,
hijl
hil
(2.22)
onde hijk é a altura do triângulo t1 relativa a aresta ejk , além disso,
∠ (∇βi , ∇βj ) |t1 = π − θk , ∠ (∇βi , ∇βj ) |t2 = π − θl ,
logo
Z
M
h∇βi , ∇βj i =
Z
h∇βi , ∇βj i =
t1 ∪t2
Z
h∇βi , ∇βj i +
t1
Z
(2.23)
h∇βi , ∇βj i =
t2
= h∇βi , ∇βj i |t1 Ak + h∇βi , ∇βj i |t2 Al ,
onde Ak e Al são as áreas dos triângulos t1 e t2 respectivamente. Usando
2.21, 2.22 e 2.23, sabendo também que
cos π − θk = − cos θk e
cos π − θl = − cos θl ,
obtemos
Z
1
1
1 1
1 1
h∇βi , ∇βj i = − i j cos θk |vj vk | hijk − i j cos θl |vj vl | hijl
2
2
h
h
h
h
M
jk ik
jl il
15
1
1
h∇βi , ∇βj i = −
2
M
Z
mas
hjik
cos θ
e
|vj vl |
1
= i ,
l
sen θ
hjl
1
|vj vk |
= i
k
sen θ
hjk
daı́
Z
h∇βi , ∇βj i = −
M
1
k
1
2
jk +
hjil
!
l
cos θ |vj vl | ,
cos θk
cos θl
+
k
sen θ
sen θl
.
Concluı́mos então que
Z
1
h∇βi , ∇βj i = − cotg θk + cotg θl ,
2
M
ou seja,
wij = −
2. caso i = j
1
cotg θk + cotg θl .
2
(2.24)
Para este caso temos wij 6= 0 apenas nos triângulos em que vi é um vértice
incidente (ver figura 2.4).
Como consequência temos:
vi
vk
θk
θl
vl
Figura 2.4: [vi , vk , vl ] representa um triângulo qualquer de M em que o vértice
vi ∈ M é incidente, θk e θl ângulos opostos ao vértice vi em cada um desses
triângulos.
wii =
Z
M
h∇βi , ∇βi i =
Z
h∇βi , ∇βi i =
∪[vi ,vk ,vl ]
wii =
X
X
[vi ,vk ,vl ]∈M
Z
h∇βi , ∇βi i
[vi ,vk ,vl ]
2
|∇βi | A[vi ,vk ,vl ] ,
[vi ,vk ,vl ]∈M
onde [vi , vk , vl ] são triângulos de M em que vi é incidente e A[vi ,vk ,vl ] a
área desses triângulos.
16
Usando mais uma vez o lema 3 e observando a figura 2.4 concluı́mos que
wii =
1
X
(hikl )2
[vi ,vk ,vl ]∈M
X
wii =
[vi ,vk ,vl ]∈M
X
wii =
|vk vl | hikl
=
2
×
hikl
hikl
+
tg θk
tg θl
[vi ,vk ,vl ]∈M
X
[vi ,vk ,vl ]∈M
/ 2hikl
1
1
+
k
2tgθ
2 tg θl
|vk vl |
2hikl
.
Finalmente para este caso temos,
wii =
1
2
X
[vi ,vk ,vl ]∈M
cotg θk + cotg θl ,
(2.25)
onde θk e θl são ângulos opostos ao vértice vi .
Alcançamos o objetivo deste capı́tulo apresentando um laplaciano discretizado,
caracterizado pelos pesos wij . Veremos no próximo capı́tulo que se a malha
de triângulos considerada for de Delaunay esses pesos serão todos negativos e
como consequência teremos um laplaciano discretizado obedecendo ao princı́pio
do máximo.
17
Capı́tulo 3
Triangulação de Delaunay
Intrı́nseca
Finalizamos o capı́tulo anterior com as seguintes expressões para os pesos wij :
wij = −
quando i 6= j e,
wii =
1
2
1
cotg θk + cotg θl
2
X
cotg θk + cotg θl
[vi ,vk ,vl ]∈M
(3.1)
(3.2)
quando i = j. Vimos também que se os pesos wij forem todos negativos, o
laplaciano discretizado satisfaz a um princı́pio do máximo. Entretanto, estes
pesos dependem dos ângulos da triangulação. Como um mesmo conjunto de
vértices V = {v1 , v2 , · · · , vn } possui mais de uma triangulação M, não está
claro se o princı́pio do máximo será sempre satisfeito ou não.
Neste capı́tulo mostraremos que, para uma aresta eij ∈ M com extremidades
nos vértices vi , vj ∈ V , wij < 0 se, e somente se, eij satisfaz a um critério
de Delaunay local. Logo, para termos um laplaciano discreto obedecendo ao
princı́pio do máximo basta que todas as arestas da malha M sejam de Delaunay.
Descreveremos também um algoritmo que aplica modificações locais na malha
de modo a obter uma triangulação onde todas as arestas sejam de Delaunay.
Devido as várias propriedades da triangulação de Delaunay e a sua importância no resultado final deste trabalho, separamos este capı́tulo para descrevermos esta ferramenta importantı́ssima em computação gráfica, bem como
em outras áreas. A sua importância está na propriedade de maximizar o menor
de todos os ângulos internos dos seus triângulos.
3.1
Triangulação de Delaunay no Plano
Na definição 5 do capı́tulo anterior, foi apresentado a definição de uma triangulação M de um conjunto de pontos do plano R2 , iniciaremos esta seção
apresentando alguns resultados desta definição em seguida definiremos a triangulação de Delaunay.
18
Em geral M não é única, mas embora não seja única a quantidade de
triângulos e de arestas em qualquer triangulação de V é a mesma, é o que
veremos no teorema a seguir.
Teorema 4. Seja V = {v1 , v2 , ....., vn } um conjunto de pontos do plano, toda
triangulação de V possui exatamente 2n − u − 2 triângulos e 3n − u − 3 arestas,
onde u é o número de pontos de V que estão na fronteira do fecho convexo de
M denotado por Conv(V ) ( maior polı́gono convexo do plano que contém todos
os vértices de V ).
Demonstração. Uma triangulação de V determina uma subdivisão do plano em
T + 1 faces, onde T é o número de triângulos da triangulação. Portanto, pela
fórmula de Euler temos que
n + (T + 1) = a + 2,
(3.3)
onde a é o número de arestas. As u arestas da fronteira de Conv(V ) são comuns
a um triângulo e à face externa. Todas as demais arestas figuram em exatamente
dois triângulos. Logo, somando o número de arestas de todas as faces, processo
no qual cada aresta é contada duas vezes, encontramos
3T + u = 2a,
daı́
3T + u
.
2
Substituindo em 3.3 e resolvendo para T e depois para a, concluı́mos a demonstração com
T = 2n − u − 2
e
a = 3n − u − 3.
a=
É consequência imediata deste teorema que o problema de triângulação é
linear no número de elementos de V .
Embora todas as triangulações de V tenham o mesmo número de triângulos,
em determinadas aplicações é importante que os triângulos da triangulação sejam ”robustos”, ou seja, que não apresentem triângulos com ângulos internos
muito pequenos.
A triangulação de Delaunay é caracterizada por ter cada uma de suas arestas
um cı́rculo circunscrito aos seus vértices e que não possui nenhum outro vértice
em seu interior, são os chamados “cı́rculos vazios”. Esta propriedade será a
base principal na construção do algoritmo que descreveremos a seguir e uma
consequência importante dela é maximizar o menor de todos os ângulos internos
dos triângulos que compõe a triangulação. Antes de apresentar o algoritmo
precisamos de três definições: aresta de Delaunay, triângulo de Delaunay e
triangulação de Delaunay.
Definição 8. Seja V um conjunto de pontos do plano e M uma triangulação
de V , uma aresta eij é dita de Delaunay quando existe um cı́rculo Cij circunscrevendo eij , ou seja, circunscrevendo os vértices vi e vj excluindo todos os
outros vértices de V . Semelhantemente um triângulo tijk ∈ M é dito Delaunay,
quando existe um cı́rculo Cijk circunscrevendo tijk , ou seja, circunscrevendo os
vértices vi ,vj e vk , excluindo todos os outros vértices de V .
19
A seguir apresentamos a definição formal da triangulação de Delaunay.
Definição 9 (Triangulação de Delaunay). Seja V um conjunto de pontos do
plano e M uma triangulação de V , M é denominada Triangulação de Delaunay
se todos os seus triângulos são de Delaunay.
Iniciamos este capı́tulo mensionando o teorema seguinte, vejamos agora a
sua demonstração.
Teorema 5. Seja V um conjunto de pontos do plano e M uma triangulação de
V , todos os triângulos de M são de Delaunay se, e somente se, todas as arestas
de M são de Delaunay.
Demonstração. Primeiro vamos mostrar que se todos as arestas e ∈ M são de
Delaunay então todos os triângulos t ∈ M são de Delaunay.
Suponhamos por absurdo que exista tijk ∈ M que não seja de Delaunay,
então existe um cı́rculo Cijk que circunscreve tijk possuindo um vértice v ∈ V
de M no seu interior (ver figura 3.1).
vi
Cijk
vj
vk
Figura 3.1: Cı́rculo circunscrevendo os vértices vi ,vj e vk
O fato de v está no interior de Cijk implica que uma das arestas de tijk não
é de Delaunay. Sem perda de generalidade suponhamos que v esteja no interior
da região limitada por eij e o cı́rculo Cijk . Afirmamos que qualquer cı́rculo
que circunscreve eij possui no seu interior ou o vértice vk ou o vértice v, como
esse fato contraria a hipótese, concluı́mos que todos os triângulos t ∈ M são de
Delaunay.
Mostremos agora que se todas as arestas e ∈ M são de Delaunay então todos
os triângulos tijk ∈ M são de Delaunay. Seja tijk ∈ M um triângulo qualquer,
como M é de Delaunay, existe um cı́rculo Cijk que circunscreve tijk e que não
possui nenhum outro vértice de M no seu interior, este cı́rculo Cijk mostra que
todas as arestas de t são de Delaunay, e como tijk é um triângulo qualquer,
temos que todas as arestas e ∈ M são de Delaunay.
Vejamos agora através do lema a seguir qual a condição para que uma aresta
e ∈ M seja de Delaunay:
Lema 6. Sejam vi vj e vk vz segmentos do plano que se interceptam em O. Então
para existir um cı́rculo passando por vi e vj e tal que vk e vz sejam exteriores
a esse cı́rculo, é necessário e suficiente que os ângulos do quadrilátero vi vj vk vz
sejam tais que φikz + φjkz > π (equivalentemente φkij + φzij < π).
20
Demonstração. Suponhamos que exista tal cı́rculo (ver figura 3.2). Sejam vk′ e
′
′
vz′ as intersecções do segmento vk vz com o cı́rculo. Então φkij +φzij < φkij +φzij =
π (já que o quadrilátero vi vk′ vj vz′ é inscritı́vel em um cı́rculo). Reciprocamente,
Figura 3.2: Cı́rculo circunscrevendo vi e vj com vk e vz externo.
se φkij + φzij < π, existem pontos vk′ e vz′ tomados sobre o segmento vk vz de
′
′
tal modo que φkij + φzij < φkij + φzij = π. Estes pontos, juntamente com vi e vj
determinam um cı́rculo que exclui vk e vz .
Na próxima seção apresentamos um algoritmo que tem como entrada uma
triangulação qualquer de um conjunto de pontos do plano e retorna a triangulação de Delaunay.
3.2
Algoritmo de Flip de Arestas no Plano
O clássico algoritmo que iremos descrever utiliza uma triangulação qualquer
para chegar à triangulação de Delaunay. Este algoritmo é baseado na mudança
de aresta, segundo o qual a aresta que não obedece ao critério local de Delaunay é trocada, ou seja, para esta aresta deve existir um cı́rculo que contenha
apenas os seus vértices e que não possua nenhum outro vértice no seu interior (são chamados cı́rculos vazios). O algoritmo pára quando todas as arestas
obedecerem a este critério.
Um importante dado de entrada da malha são as combinações de seu grafo
de arestas, assim como o comprimento de cada aresta. Com esses dados a verificação do critério local de Delaunay pode ser feita apenas com o comprimento
das arestas e das combinações locais das mesmas. Um fato importante é que este
procedimento não muda a geometria da malha e não é difı́cil verificar que este
procedimento pára no momento em que todas as arestas da malha obedecem ao
critério local de Delaunay, como faremos a seguir.
O algoritmo que descreveremos a partir de agora tem como entrada uma
triangulação qualquer M0 de um conjunto V de pontos do plano e como saı́da
a triangulação de Delaunay.
Seja o complexo K = (V, E, T ) de vértices, arestas e triângulos, seja também
l uma métrica Euclidiana , tal que, lij = kvi − vj k é o comprimento da aresta
eij ∈ E, a partir destes dados criamos o código seguinte.
Entrada: (K, l)
Saida : (K ′ , l′ ) triangulação de Delaunay
∀e ∈ E Mark(e)
Stack s ← E
While s 6= φ do
Desempilhe eij de s e desmarque eij
If !Delaunay(eij ) then
ekl ← F lip(eij ) e calcule ekl
For all e ∈ {ekj , ejl , eli , eik } do
If !Mark(e) Then
Mark(e) e coloque em s
end If
end For
21
end If
end While
Nosso objetivo a seguir é mostrar que este algoritmo para em um número
finito de etapas.
Teorema 7. O algoritmo de flip de arestas no plano pára em um número finito
de passos.
Demonstração. Seja V um conjunto de pontos do plano e M uma triangulação
de V . Para cada triângulo t ∈ M tome A(t) como a área do cı́rculo que
circunscreve t. Esta área pode ser calculada com os comprimentos a,b e c das
arestas do triângulo t, utilizando a fórmula seguinte,
16πA =
(abc)2
s(s − a)(s − b)(s − c)
onde s = 21 (a + b + c) é o semi-perı́metro de t. Tomemos agora umaPfunção H
definida no conjunto de todas as triangulações de V , sendo, H(M) = t∈M A(t).
Seja M0 , M1 , M2 , M3 , ..... uma sequência de triangulações de V , onde M0 é
uma triangulação qualquer e M1 , M2 , M3 , ..... uma sequência de triangulações
originadas a partir do flip de aresta (algoritmo). A partir do teorema 6.1 de
[10] podemos concluir que H(M0 ), H(M1 ), H(M2 ), ....., é uma sequência estritamente decrescente, (geometricamente é fácil verificar este fato por meio
da figura 3.3). Como o número de triangulações de um conjunto de pontos do
C3
C1
C2
C4
Figura 3.3: C1 e C2 são os cı́rculos circunscritos aos triângulos da malha antes
do flip. C3 e C4 são os cı́rculos circunscritos aos novos triângulos inseridos na
malha. Observe que C3 ∪ C4 ⊂ C1 ∪ C2 e isso mostra que H(Mn ) > H(Mn+1 ).
plano é finita e que H(M0 ), H(M1 ), H(M2 ), ....., é uma sequência estritamente
decrescente, ela atinge um valor mı́nimo em um número finito de flips, o que
mostra que o algoritmo de flip no plano pára em um número finito de etapas
(flips).
A verificação da condição local de Delaunay é fundamental para a implementação do algoritmo, bem como o cálculo do comprimento da nova aresta,
caso uma aresta seja trocada. Tanto o teste de Delaunay como o cálculo
da nova aresta, pode ser feito usando apenas o comprimento das arestas dos
triângulos incidentes da aresta de teste. Vamos verificar esta afirmação utilizando a proposição que segue.
22
Proposição 1. Sejam a,b e c os comprimentos dos lados de um triângulo e φ
o ângulo oposto ao lado de comprimento a, então
s
(a − b + c) (a + b − c)
φ
.
tg =
2
(a + b + c) (−a + b + c)
Demonstração. Da trigonometria sabemos que
tg2
sen2 φ2
φ
1 − cos φ
=
.
=
φ
2
2
1 + cos φ
cos 2
(3.4)
φ
Figura 3.4: Triângulo com lados de comprimentos a, b e c e ângulo φ oposto ao
lado de comprimento a.
Usando a lei dos cossenos
cos φ =
b2 + c2 − a2
2bc
na equação 3.4 obtemos
φ
b2 + c2 − a2
b2 + c2 − a2
= (1 −
)/(1 +
)
2
2bc
2bc
− b2 − 2bc + c2 − a2
2bc − b2 − c2 + a2
=
=
2
2bc + b2 + c2 − a2
(b + c) − a2
tg2
2
=
a2 − (b − c)
(a + c − b) (a + b − c)
=
,
(b + c − a) (b + c + a)
(b + c − a) (b + c + a)
daı́ concluı́mos o teorema com
φ
tg =
2
s
(a + c − b) (a + b − c)
.
(b + c − a) (b + c + a)
(3.5)
Com a equação 3.5 calculamos a cotangente de φ utilizando a seguinte relação
trigonométrica
1 − tg φ2
.
(3.6)
cot gφ =
2tg φ2
23
Utilizando 3.6 e associando a cada aresta eij incidente a dois triângulos tijk e
tijz o peso seguinte,
wij = −
1
cotg φkij + cotg φzij ,
2
(3.7)
será fácil verificar, conforme faremos a seguir, que se cotg φkij + cotg φzij for não
negativo a aresta eij satisfaz o critério local de Delaunay, caso contrário esta
aresta deve ser trocada por ekz e calculado o comprimento desta nova aresta.
Pelo lema 6 eij é de Delaunay se, e somente se φkij + φzij ≤ π
Proposição 2. . Sejam vi , vj , vk e vz vértices coplanares de um quadrilátero,
então φkij + φzij ≤ π se e somente se
cotg φkij + cotg φzij ≥ 0
vi
vk
vz
φzij
φkij
vj
Demonstração. Temos que
cotg φkij + cotg φzij =
=
cos φkij
+
k
sen φij
cos φkij sen φzij + cos φzij sen φkij
sen φkij sen φzij
=
cos φzij
sen φzij
sen (φkij + φzij )
sen φkij sen φzij
,
logo
cotg φkij + cotg φzij =
sen (φkij + φzij )
sen φkij sen φzij
,
(3.8)
como
φkij < π
φzij < π
e
temos
sen φkij sen φzij > 0,
(3.9)
sen (φkij + φzij ) ≥ 0.
(3.10)
se φkij + φzij ≤ π então
De 3.9 e 3.10 concluı́mos que
cotg φkij + cotg φzij =
sen (φkij + φzij )
24
sen φkij sen φzij
≥ 0.
Vamos mostrar agora que se cotg φkij + cotg φzij ≥ 0 então φkij + φzij ≤ π.
Se cotg φkij + cotg φzij ≥ 0 então o numerador e o denominador de 3.8 tem
o mesmo sinal ou o numerador é zero, mas sen φkij sen φzij é sempre positivo
conforme visto em 3.9, logo sen (φkij + φzij ) ≥ 0, concluı́mos então que φkij +φzij ≤
π.
Quando é feito o flip de uma aresta é necessário calcular o comprimento
da nova aresta a ser inclusa na malha. A seguir utizaremos mais uma vez
a tangente do arco metade para determinar uma expressão que permite este
cálculo, envolvendo apenas os comprimentos das arestas dos novos triângulos
incidentes a esta nova aresta. A partir da seguinte relação trigonométrica,
α
φ
Figura 3.5: f é o comprimento da nova aresta que pode ser calculado a partir
dos comprimentos a, b, c, d, e.
tg
φ+α
2
=
tg φ2 + tg α2
1 − tg φ2 tg α2
chegamos a,
(3.11)
b2 + c2 − 2 cos (φ + α)
(3.12)
cos(φ + α) =
1 + tg 2
e concluı́mos que,
f=
p
1 − tg 2
φ+α
2
φ+α
2
Nesta seção apresentamos um algoritmo que retorna a triangulação de Delaunay a partir de uma malha de triângulos no plano R2 , nosso objetivo na
próxima seção é estender este algoritmo para malhas no R3 .
3.3
Extensão do Algoritmo de Flip Para Malhas
no Espaço R3
Na seção anterior apresentamos o algoritmo de flip de aresta para obter a triangulação de Delaunay de um conjunto de pontos V ⊂ R2 , nosso objetivo
a partir de agora é generalizar este algoritmo para uma malha de triângulos
M do R3 , onde o conjunto de pontos singulares V é também o conjunto de
vértices da triangulação. Veremos que de forma idêntica ao algoritmo de flip de
25
aresta no plano, pode-se construir a estrutura de uma triangulação de Delaunay
intrı́nseca a uma superfı́cie S ⊂ R3 a partir de uma malha triangular qualquer
de S, isso porque, ao tomar um par de triângulos adjacentes no R3 sua geometria intrı́nseca pode ser visualizada pelo mapeamento isométrico deles no plano
(veja figura 3.3).
planificacao
resultado
flip
Figura 3.6: Dado um par de triângulos no R3 sua geometria intrı́nseca pode ser
visualizada pelo mapeamento isométrico no plano. Uma aresta que não é de
Delaunay é trocada e inserida uma nova aresta na superfı́cie original.
Portanto, dada uma malha qualquer de triângulos no espaço R3 , uma aresta
que violar o critério local de Delaunay é trocada, e calculado o comprimento
da nova aresta, exatamente como feito para o algoritmo de flip no plano. Esse
processo é repetido até que todas as arestas obedeçam ao critério local de Delaunay.
Sabemos que no caso planar o número de triangulações de um conjunto de
pontos é finito, e esse fato foi fundamental para que o algoritmo de flip no plano
parasse em um número finito de etapas. Já no caso do espaço R3 o número de
triangulações geodésicas de um conjunto de pontos possivelmente é infinito, no
entanto veremos que mesmo com essa complexidade a mais o algoritmo de flip
também pára em um número finito de etapas. Antes de demonstrarmos essa
afirmação, vejamos algumas definições que serão mencionadas.
Já falamos anteriormente em uma superfı́cie plana por partes, vejamos agora
qual é a sua definição.
Definição 10 (Superfı́cie Plana por Partes). Um espaço métrico compacto M
é uma superfı́cie com singularidades cônicas, ou Superfı́cie Plana por Partes
de forma mais simples, se todo ponto admite uma vizinhança isométrica a um
disco do plano euclidiano ou a um cone euclidiano.
Embora a definição anterior seja bem abrangente, neste trabalho nosso interesse está em malhas de triângulos, que são casos particulares de superfı́cies
planas por partes.
Definição 11 (Geodésica). Considere M uma superfı́cie plana por parte, uma
geodésica em M é uma curva α : [0, 1] → M que localmente minimiza distâncias.
A distância geodésica entre dois pontos p, q ∈ M é o comprimento do menor
caminho unindo p a q.
O próximo teorema mostra que o algoritmo de Flip em malhas no espaço
euclidiano R3 também pára em um número finito de etapas, mas antes de apresentá-lo precisamos da proposição seguinte que será de fundamental importância
para a sua demonstração.
Proposição 3. Seja M uma superfı́cie plana por partes, para algum par de pontos p, q ∈ M e um número real L > 0, o número de geodésicas de comprimento
menor que L unindo p a q é finito.
26
Demonstração. Visto que M é uma superfı́cie plana por partes, duas geodésicas
unindo p e q ou são não-homotópicas ou juntas elas formam a fronteira de uma
região contendo pelo menos uma singularidade (veja a prova no corolário 2 de
[9]). Em particular, se α é uma geodésica unindo p e q, podemos encontrar uma
vizinhança U de α formada por arcos simples unindo p e q em M no qual α é
a única geodésica em U . Agora seja Apq (L) uma coleção de todas as geodésicas
α : [0, 1] → M parametrizada de velocidade constante unindo p e q e de comprimento menor ou igual a L. Em particular, a famı́lia Apq (L) é uniformemente
de Lipschitz (com a constante de Lispchitz L) e pelo teorema de Arzela-Ascoli
é um conjunto compacto (da topologia uniforme). Esses argumentos mostram
que cada α é um ponto isolado em Apq (L), portanto Apq (L) é um conjunto
finito.
Como consequência temos o seguinte corolário:
Corolário 8. O conjunto dos comprimentos de todas as geodésicas unindo pares
de pontos em uma superfı́cie plana por partes M é um subconjunto discreto em
R.
Mostraremos agora que o algoritmo de flip de arestas em malhas no espaço
euclidiano R3 também pára em um número finito de passos. Já sabemos que
a geometria intrı́nseca de um par de triângulos adjacentes no espaço euclidiano
R3 pode ser visualizada pelo seu mapeamento isométrico no plano, esse fato
permite construir toda a demonstração de forma idêntica ao teorema 7 da seção
anterior no caso planar.
Teorema 9. O algoritmo de flip de arestas em malhas no espaço euclidiano R3
pára em um número finito de passos.
Demonstração. Seja M uma triangulação geodésica de uma superfı́cie plana por
partes S. Para cada triângulo T de M denotemos por A(T ) a área do cı́rculo
circunscrito a uma cópia isométrica de T no plano R2 . Definamos agora uma
função H no conjunto de todas as triangulações geodésicas por
X
H(M) =
A(T ),
T ∈M
pelo corolário 8 a imagem de H é um subconjunto discreto de R. Seja M, M1 ,
M2 ,... uma sequência de triangulações geradas pelo algoritmo de flip, a partir do
teorema 6.1 de [10] podemos concluir que H(M), H(M1 ), H(M2 ), ....., é uma
sequência estritamente decrescente, como a imagem de H é um subconjunto
discreto de R esta sequência atingirá um mı́nimo em um número finito de passos,
o que mostra que o algoritmo de flip pára em um número finito de etapas.
Com este teorema terminamos a parte teórica do nosso trabalho. No próximo
capı́tulo estudaremos algumas aplicações do laplaciano discreto em computação
gráfica.
27
Capı́tulo 4
Aplicações do Laplaciano
Discreto
Neste capı́tulo apresentaremos três aplicações do laplaciano: a primeira tem
como objetivo obter parametrizações de malhas para efeito de mapeamento de
textura; a segunda consiste na suavização de superfı́cies por meio do processo de
difusão, ou seja, eliminação de ruı́dos ou regiões de altas frequências da malha;
a terceira e última aplicação visa identificar formas e simetrias de objetos por
meio das curvas de contorno associadas às autofunções do laplaciano.
Todas as figuras apresentadas nestas aplicações foram geradas a partir das
implementações que fizemos no Matlab.
4.1
Mapeamento de Textura
Nosso objetivo nesta aplicação é obter parametrizações ϕ : M → Ω ⊂ R2 de
malhas de triângulos M para efeito de mapeamento de textura. Em particular
é necessário que ϕ seja injetiva para ser uma parametrização de M.
Existe uma teoria pronta e bem conhecida para o contexto contı́nuo, o teorema de Radó-Kneser-Choquet:
Teorema 10 (Radó-Kneser-Choquet-RKC). Seja D uma região fechada do R2 .
Suponhamos que uma função ϕ : D → Ω ⊂ R2 é um mapeamento harmônico,
que mapeia ∂D homeomorficamente em ∂Ω, sendo Ω ⊂ R2 uma região convexa
do R2 . Então ϕ é injetiva.
Uma função ϕ é dita harmônica se ∆ϕ = 0. Em particular quando D ⊂ R2
temos ϕ = (ϕ1 (x, y) , ϕ2 (x, y)) e portanto:
∆ϕ1 = 0
∆ϕ2 = 0.
Será que existe uma teoria análoga deste problema no caso discreto? Ou
seja, uma parametrização
ϕ : M → Ω ⊂ R2 ,
28
injetiva tal que,
∆d ϕ1 = 0
∆d ϕ2 = 0.
Sendo ∆d o laplaciano discreto definido no capı́tulo 2 e M uma malha de
triângulos. A resposta para esta pergunta é sim. Essa teoria já foi construı́da
por Michael Floater que demonstrou em [13] um teorema análogo ao de RKC
para o caso discreto garantindo neste teorema que a parametrização seja injetiva.
Para entender o teorema de Floater precisamos da seguinte definição:
Definição 12. Seja f : M → R uma função linear por partes, e suponha
que
P interior de M, existam pesos λij > 0, tais que
P para todo vértice vi no
λ
=
1
e
f
(v
)
=
i
vj ∈Ni λij f (vj ), onde
vj ∈Ni ij
Ni = {vj ; eij é uma aresta de M } .
Denominaremos a função f de função de combinação convexa.
De forma similar denominaremos φ = (f, g) : M → R2 , um mapeamento de
combinação convexa se ambas as componentes f (x, y, z), g(x, y, z) de φ forem
funções de combinação convexa.
Definição 13 (Triangulação fortemente conexa). Uma triangulação M é dita
fortemente conexa se toda aresta interna de M possui pelo menos um vértice
no interior de M.
Teorema 11 (Teorema de Floater). Supondo M uma triangulação fortemente
conexa e φ : M → Ω um mapeamento de combinação convexa que mapeia
homeomorficamente a fronteira ∂M na fronteira ∂Ω de uma região convexa
fechada Ω ⊂ R2 . Então φ é injetiva.
O ponto chave da teoria de Floater é que as funções de combinações convexas
são as análogas discretas das funções harmônicas no caso contı́nuo, visto que estas funções também satisfazem ao princı́pio do máximo, conforme mostraremos
a seguir, tendo esse fato similar importância quanto tem as funções harmônicas
no caso contı́nuo.
No corolário a seguir veremos que se uma função f : M → R, onde M é uma
malha de triângulos, é de combinação convexa, então ela satisfaz ao princı́pio
do máximo.
Corolário 12 (Princı́pio do Máximo Discreto). Seja f : M → R uma função
de combinação convexa, onde M é uma malha de triângulos. Se para algum
vértice vi no interior de M tivermos f (vi ) ≥ f (vj ) qualquer que seja vj ∈ Ni
então f (vi ) = f (vj ) para todo vj ∈ Ni .
Demonstração. Por hipótese f (vi ) > f (vj′ ) para todo vj ∈ Ni , suponhamos por
absurdo que exista vj′ ∈ Ni tal que f (vi ) > f (vj′ ). Então,
X
λij (f (vi ) − f (vj )) > 0.
(4.1)
vj ∈Ni
Por outro lado temos
X
X
X
λij f (vj ) = f (vi ) − f (vi ) = 0,
λij f (vi ) −
λij (f (vi ) − f (vj )) =
vj ∈Ni
vj ∈Ni
vj ∈Ni
29
o que nos leva a uma contradição por 4.1. A contradição surgiu ao considerar a
existência de vj′ ∈ Ni , logo f (vi ) = f (vj ) para todo vj ∈ Ni .
No capı́tulo 2 quando definimos ∆d associamos os pesos wij a cada aresta
eij de uma determinada malha de triângulos M, sendo
wij = −
1
cotg θk + cotg θl , i 6= j,
2
(4.2)
onde θk e θl são ângulos opostos a aresta eij dos triângulos em que eij é uma
aresta incidente e
X
1
wii =
cotg θk + cotg θl
,
i = j,
(4.3)
2
[vi ,vk ,vl ]∈M
onde [vi , vk , vl ] ∈ M são triângulos em que vi é incidente e θk , θl ângulos opostos
ao vértice vi nesses triângulos. Foi visto também no capı́tulo 3 proposição 2 que
cotg θk + cotg θl
é não negativo se, e somente se θk + θl < π, ou seja, se e somente se eij satisfaz
ao critério local de Delaunay.
Aqui associaremos a cada aresta eij os coeficientes
λij =
−wij
,
wii
o que nos leva ao seguinte resultado:
Corolário 13. Seja f : M → R uma função linear por partes tal que ∆d f = 0.
Se M é uma triangulação de Delaunay então f é uma função de combinação
convexa.
Este corolário mostra a importância da triangulação de Delaunay na discretização do laplaciano ∆d que escolhemos, pois, basta que a malha M de
triângulos seja de Delaunay e fortemente conexa para que um mapeamento
linerar por partes φ : M → Ω ⊂ R2 seja injetivo.
Descreveremos a partir de agora o método linear que utilizamos para implementar a parametrização de uma malha M no espaço R3 .
Se φ é um mapeamento linear por partes, ele é completamente determinado
pelos seus pontos φ(vi ) ∈ R2 para vi ∈ V .
Denotaremos por VI o conjunto de vértices interiores de M e VB o conjunto
de vértices da fronteira. O primeiro passo do método é escolher alguns pontos
φ(v) ∈ R2 , para v ∈ VB , tal que a fronteira ∂M de M seja mapeada homeomorficamente num polı́gono φ(∂M) do plano R2 . O segundo passo é escolher
para cada vi ∈ VI um conjunto de valores λij estritamente positivo, tal que
X
λij = 1.
vj ∈Ni
Na implementação que fizemos foram utilizados
λij =
−wij
,
wii
30
Figura 4.1: Parametrização injetiva após flips de arestas
onde wij esta definido em 4.2 e wii em 4.3. Assim os pontos φ(vi ) ∈ R2 são as
únicas soluções do sistema linear de equações
X
λij φ(vj ) vi ∈ VI
(4.4)
φ(vi ) =
vj ∈Ni
Como consequência temos que cada ponto φ(vi ) é uma combinação convexa
de seus vizinhos φ(vj ).
Vamos reescrever 4.4 na seguinte forma,
X
X
φ(vi ) −
λij φ(vj ) =
λij φ(vj )
(4.5)
vj ∈(Ni ∩VI )
vj ∈(Ni ∩VB )
Reescrevendo agora 4.5 na forma de uma equação matricial temos:
Ax = b,
onde x = (φ (vj ))vj ∈VI é um vetor coluna a determinar, b é um vetor coluna
cujo os elementos são formados pelo segundo membro da equação 4.5, e a matriz
A = (aij )vi ,vj ∈VI tem dimensão n × n, com n = |VI |, e elementos
aij =
1,
vi = vj
−λij ,
vj ∈ Ni
0,
outros casos
de posse da matriz A e da matriz b no Matlab é muito fácil determinar a
matriz incógnita x, basta aplicar o comando x = A/b e o Matlab retorna a
matriz x cujas linhas são as coordenadas de cada ponto φ(vi ), onde vi ∈ VI .
As figuras de 4.1 a 4.3 são parametrizações de algumas malhas no R3 a partir
das implementações que fizemos no Matlab.
A figura 4.4 é um exemplo que confirma a necessidade da triangulação da
malha ser de Delaunay. Neste caso a malha original não era de Delaunay e
como consequência obteve-se uma parametrização não injetiva, já após o flip de
aresta a parametrização obtida em 4.5 é injetiva. Como o nosso objetivo final é
aplicar uma textura sobre a malha, esse resultado mostra a necessidade real de
se aplicar o algoritmo de flip de aresta na malha, caso contrário poderá ocorrer
problemas na aplicação da textura.
31
Figura 4.2: Parametrização injetiva após flips de arestas
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figura 4.3: Parametrização injetiva após flips de arestas
Figura 4.4: parametrização não injetiva sem fazer flips de arestas.
Figura 4.5: parametrização injetiva após flips de arestas.
Obtida a parametrização da malha aplicamos uma textura sobre a mesma.
Nas figuras 4.7 e 4.9 a aplicação da textura foi feita por meio de um programa desenvolvido usando a biblioteca OpenGL, precisando apenas informar os vértices
de cada face da triangulação obtida da parametrização via algoritmo descrito
acima e os vértices das suas correspondentes faces na malha. Feito isso, o pro32
grama se encarregou de aplicar a textura.
Figura 4.6: Parametrização injetiva após flips de aresta
Figura 4.7: Aplicação da textura após parametrização
Figura 4.8: Parametrização injetiva do cubo
33
Figura 4.9: Aplicação da textura após parametrização
4.2
Suavização de Superfı́cies
Conforme mencionado na introdução, o laplaciano aparece no segundo membro
da equação de difusão
∂f
= λ∆f.
∂t
Para suavizar uma malha, basta resolver essa equação usando como condição
inicial as funções coordenadas da malha. Com um ligeiro abuso de notação,
representaremos este processo pela equação
∂M
= λ∆(M),
∂t
(4.6)
onde M representa uma malha de triângulos.
Voltando a equação 4.6 e integrando M com relação a t, uma pequena
perturbação será rapidamente propagada em sua vizinhança, essa perturbação
suaviza as regiões de M com altas frequências, enquanto a forma principal da
malha é levemente afetada.
Desenvolvendo 4.6 obtemos:
Mn+1 = Mn + λdt∆(Mn ),
(4.7)
onde Mn = M(t) e Mn+1 = M(t + dt), logo uma sequência de malhas pode
ser construı́da através da seguinte equação:
Mn+1 = (I + λdt∆)Mn .
(4.8)
Com a matriz de ∆d definido no capı́tulo 2 e a equação 4.7 criamos um código
para gerar uma sequência de malhas cada vez mais suaves. A implementação
requer como dados de entrada o comprimento λdt, a matriz de ∆d , a malha M
e o número n de malhas a serem geradas, e como saı́da tem-se uma malha após
n iterações bem mais suavizada.
34
A figura 4.10 é um exemplo de malha com várias rugosidades que foram
suavizadas após uma sequência de iterações a partir do código que implementamos no Matlab:
Figura 4.10: Malha original, malha após 30 iterações, malha após 50 iterações
e malha após 100 iterações de suavização, respectivamente. Para todas λdt =
0, 008
A figura 4.11 mostra que é possı́vel obter resultados similares aos anteriores
aumentando o valor de λdt o que possibilita a dinuição do número de iterações
de suavização e como consequência tem-se um menor tempo de processamento.
Figura 4.11: Esquerda malha original e a direita malha após 80 iterações de
suavização sendo λdt = 0, 01
A figura 4.12 é um exemplo claro de uma malha com regiões de altı́ssimas
frequências, suavizadas após 60 iterações.
Figura 4.12: Esquerda malha com altissı́mas frequências e a direita malha após
60 iterações de suavização, sendo λdt = 0, 008
35
4.3
Identificação de Formas e Simetrias de Malhas
Nesta seção nosso objetivo é determinar curvas de nı́vel das autofunções do laplaciano que permitem identificar formas e simetrias em malhas bem como uma
certa ordenação dos seus vértices. Faremos um hiato para falar dos conjuntos
nodais inicialmente.
4.3.1
Conjuntos Nodais
Estes conjuntos foram estudados pelo fı́sico Ernest Chladni e publicado em
seguida no seu livro “Discoveries Concerning the Theories of Music”. As experiências de Chladni consistiram em colocar areia em um prato de metal submetido a vibrações de diferentes intensidades, ele constatou que para cada vibração a areia se concentrava em determinadas regiões do prato, formando zonas
com padrões complexos. Este comportamento pode ser modelado pelas autofunções f do laplaciano que satisfazem a equação de propagação da onda
∆f = λf.
(4.9)
Muitos matemáticos têm estudado a convergência do espectro do operador
laplaciano, mas nos limitaremos aqui a estudar a geometria das autofunções
deste operador, na verdade vamos estudar osconjuntos nodais, definido como
o conjunto de zeros das autofunções, ou seja, x ∈ R2 ; f (x) = 0 , onde f é uma
autofunção do laplaciano. Este conjunto corresponde as zonas onde a areia se
acumula no prato vibrante de Chladni.
Uma forma de aproximar o operador laplaciano além da que já apresentamos
no capı́tulo 2, é utilizando aproximação linear em cada vértice pelo operador
umbrella (guarda-chuva), sendo este o laplaciano pontual
(∆u)i =
1 X
(ui − uj ) ,
m
(4.10)
vj ∈Ni
onde vj são os vizinhos do vértice vi e m o número desses vizinhos (valência).
1
Observe que os pesos wij nesta aproximação são uniformes, ou seja, wij = m
para todo vj ∈ Ni , e neste caso não é levada em consideração a geometria nem
as combinações existentes na malha, por esta razão apresentamos no Capı́tulo 2
um laplaciano discreto geométrico ∆d que leva em conta a geometria da malha.
Vejamos quais as etapas da implementação no Matlab para gerar as curvas
de nı́vel que no caso desta seção são os conjuntos nodais, ou seja, curvas no nı́vel
zero de funções definidas em malhas planares.
Criamos um código para gerar a matriz de ∆d no Matlab calculamos seus
autovetores que correspondem as autofunções lineares por parte do laplaciano
discretizado, ou seja, com o comando eigs(∆d ) é fornecido pelo Matlab uma
matriz An×m onde n representa o número de vértices da malha, e m o número
de autofunções de ∆d , sendo que cada coluna desta matriz corresponde aos
valores destas autofunções nos vértices da malha, Com esses valores usamos
interpolação linear para calcular os conjuntos nodais de cada autofunção, ou
seja, os valores x ∈ R2 em que cada autofunção se anula.
36
A figura 4.13 foi gerada a partir da implementação que fizemos no Matlab
e representa os conjuntos nodais de várias autofunções de uma malha planar
com fronteira quadrada.
Figura 4.13: Conjuntos nodais de várias autofunções de uma malha planar com
fronteira quadrada
Na próxima seção vamos estudar as curvas de nı́vel das autofunções do laplaciano em malhas no espaço R3 , mas agora não apenas no nı́vel zero. Nosso objetivo é obter para cada autofunção curvas sobre a superfı́cie em vários nı́veis, ou
seja, fixando vários valores para uma autofunção obter os pontos da superfı́cie
onde ela assume esses valores.
4.3.2
Curvas de Contorno
O primeiro autovetor da matriz do laplaciano discretizado ∆d definido no capı́tulo
2 é (1,1,...,1) e está associado ao autovalor 0. O segundo autovetor é denominado autovetor de Fiedler, e tem várias aplicações, como ordenação de vértices
e triângulos em uma malha, essa aplicação pode ser verificada em [11].
Na figura 4.14 apresentamos curvas de nı́veis do autovetor Fiedler em Três
malhas distintas, geradas a partir da implementação que fizemos no Matlab,
estas curvas também são denominadas de curvas de contorno e permitem obter
uma ordem dos vértices e triângulos da malha, como também visualizar simetrias em cada um dos objetos. A figura 4.15 mostra uma diferença nı́tida en-
Figura 4.14: Curvas de contorno do vetor de Fiedler que permitem identificar
formas e simetrias nos objetos, bem como ordenação dos seus vértices.
tre a discretização do laplaciano usando o operador umbrella e operador ∆d .
37
A geometria das curvas de nı́vel diferem bastante, enquanto no operador ∆d
temos curvas bem suaves usando o Laplaciano umbrella as curvas apresentam
ondulações.
Figura 4.15: Curvas de contorno da 7a autofunção do Laplaciano umbrella à
esquerda e à direita do Laplaciano ∆d cujo os pesos são somas de cotangentes
apresentadas no capı́tulo 2.
38
Capı́tulo 5
Conclusão
Nesta dissertação apresentamos um operador laplaciano cujo domı́nio são funções
definidas em uma superfı́cie, no entanto vimos que para obter resultados práticos
é necessário uma versão discreta desse operador, constatamos essa necessidade
na parametrização de malhas de triângulos para efeito de mapeamento de textura, onde foi mostrado que em uma malha que possui todas as suas arestas
obedecendo ao critério local de Delaunay e é fortemente conexa é possı́vel obter
uma parametrização injetiva, evitando problemas como faces sobrepostas na
aplicação da textura; também foi necessário uma discretização do laplaciano
para suavizar malhas por meio do processo de difusão, onde, malhas com altas
frequências foram suavizadas após certo número de iterações; por fim, observamos formas e simetrias em malhas de triângulos ao utilizarmos as curvas de
nı́vel associadas as autofunções do laplaciano discretizado, destacando-se nesta
aplicação o vetor de Fiedler.
Embora na escolha da discretização do laplaciano tenhamos optado por uma
que utilizasse pesos que eram somas de cotangentes, portanto uma discretização
geométrica, verificamos resultados não satisfatórios em algumas aplicações, como
distorções de grandes proporções quando aplicado a textura sobre determinadas
malhas, que poderiam ser minimizadas se escolhida uma parametrização mais
adequada. Na suavização de malhas constatamos um leve encolhimento da
malha, ou seja, uma leve diminuição de valor nas coordenadas da malha, ocasionando assim leve mudança na sua forma geral.
As dificuldades citadas anteriormente servirão de base para estudos futuros,
como encontrar parametrizações que permitam minimizar distorções quando
aplicado a textura sobre malha, assim como buscar alternativas que paralelamente aplicadas a suavização da malha permitam a invariancia na forma geral
da malha. Outra fonte de pesquisas futuras bastante promissora é a ordenação
de vértices em malhas grandes, utilizando-se das curvas de nı́vel dos autovetores
da matriz do laplaciano discretizado, dando destaque para o vetor de Fiedler,
tendo este estudo diversas aplicações.
39
Bibliografia
[1] Bobenko, A. I., and Springborn, B. A., A discrete Laplace-Beltrami
operator for simplicial surfaces, math/0503219, 2005.
[2] Indermitte, C., Liebling, Th. M., Troyanov, M., and Clemençon, H.,
Voronoi Diagrams on Piecewise Flat Surfaces and an Application to
Biological Growth, ScienceDirec-Theoretical Computer Science, 2001.
[3] Fisher M., Springborn B., Bobenko A. I., Schröder P.: An algorithm
for the construction of intrinsic Delaunay triangulations with applications to digital geometry processing, Computing vol. 81, pages 199213, 2007. .
[4] Lévy, B. Laplace-Beltrami Eigenfunctions Towards an algorithm that
“understands” geometry, IEEE-Computer society, 2006.
[5] Desbrun, M., Meyer, M.,Schr, P. Alan H.B., Implicit Fairing of Irregular Meshes using Diffusion and Curvature Flow, SIGGRAPH 99,
pages 317-324, 1999.
[6] Vito, M.C., The Spectrum of the Laplacian in Riemannian Geometry,
math.berkeley.edu,2003.
[7] Wardetzky, M., Mathur, S., Kalberer, F., and Grinspun, E. , Discrete Laplace operators: No free lunch, Eurographics Symposium on
Geometry Processing, 2007 .
[8] Perdigão, M.F., Geometria Diferencial de Curvas e Superfı́cies, Instituto de Matemática Pura e Aplicada (IMPA), 2005.
[9] M. Troyanov, Les surfaces euclidiennes Sa singularité coniques,
L’Enseignement Mathématique, Vol.32 ,1986.
[10] H. Telley, Modélisation et Simulation Bidimensionnelle de la Croissance des Polycristaux, Ecole Polytechnique Federal de Lausanne,
1989.
[11] M. Isenburg and P. Lindstrom. Streaming meshes, In IEEE Visualization, page 30, 2005.
[12] J. Rodney Biezuner, Autovalores do Laplaciano, Notas de Aula, Departamento de Matemática, UFMG.
[13] S. Michael Floater, One-to-one piecewise linear mappings over triangulations, Mathematics of Computation, Vol. 72, pags 685-696, 2003.
40
