Universidade da Beira Interior Departamento de Inform ...agomes/PhDthesis-FranciscoMorgado.pdf · ii. Tese realizada sob a ... 2004, Lecture Notes in Computer Science, Vol.3039, Springer-Verlag,

Universidade da Beira InteriorDepartamento de Informatica

Poligonizacao de Curvas eSuperfıcies Implıcitas

Nao-Homogeneas com PreservacaoTopologica

Jose Francisco Monteiro Morgado

Tese submetida para obtencao do grau de Doutor

em Engenharia Informatica

Covilha - Janeiro 2005

ii

Tese realizada sob a orientacao do

Professor Doutor Abel Joao Padrao Gomes

Professor Auxiliar do Departamento de Informatica daUniversidade da Beira Interior

iii

A Paula e aos meus filhos Francisco e Diogo

v

Agradecimentos

Ao meu orientador, Doutor Abel Joao Padrao Gomes, pela confianca, dedicacao eorientacao que muito contribuiram para a elaboracao desta tese.

Ao meu colega Frutuoso Silva, pelo debate de ideias que contribuiu para melhoraro trabalho aqui apresentado.

Ao PRODEP III (Accao 5.3 - Formacao Avancada de Docentes do Ensino Su-perior), ao IT (Instituto de Telecomunicacoes), a ESTV-ISPV (Escola Superior deTecnologia de Viseu - Instituto Superior Politecnico de Viseu) e ao Departamento deInformatica da UBI, pelos apoios concedidos, sem os quais este trabalho nao teria sidopossıvel.

Aos colegas do Departamento de Informatica da ESTV que, durante os trabalhosde doutoramento, mostraram o seu interesse e apoio atraves dos seus comentarios eboa disposicao.

A minha famılia, Paula, Francisco e Diogo, por toda a paciencia e compreensaodemonstradas nas minhas ausencias. Um muito Obrigado.

vii

Resumo

Ha sistemas matematicos (por exemplo, Maple, MathLab e Mathematica) capazesde representar curvas e superfıcies implıcitas definidas por funcoes polinomiais emR2 e em R3, respectivamente. Contudo, estes sistemas tem dificuldades notorias narepresentacao de curvas e superfıcies mais gerais definidas por funcoes reais, e taopouco preservam a sua forma topologica. Por outras palavras, estes sistemas, assimcomo a generalidade dos algoritmos actuais, nao representam correctamente curvase superfıcie implıcitas, em particular as suas singularidades (por exemplo, apices eauto-interseccoes) e componentes sem variacao de sinal, pontos isolados inclusive.

Alem dos pontos isolados e componentes sem variacao de sinal de curvas esuperfıcies, os algoritmos existentes na literatura nao geram correctamente superfıciesimplıcitas nao-homogeneas (por exemplo, as superfıcies de Steiner e de Cartan). Porexemplo, o eixo negativo dos ZZ em R3 faz parte da superfıcie de Cartan, mas naoha sistema ou algoritmo conhecido que o represente. De facto, todo o esforco deinvestigacao internacional em torno das curvas e superfıcies implıcitas tem incididoem curvas e superfıcies dimensionalmente homogeneas, independentemente de seremvariedades topologicas (i.e. manifolds) ou nao.

Do ponto de vista da topologia, uma componente sem variacao de sinal de umacurva implıcita (respectivamente, superfıcie) e o subespaco conexo maximal da curva(respectivamente, superfıcie) descrito por uma funcao real que tem sempre o mesmosinal, positivo ou negativo, numa pequena vizinhanca de cada um dos seus pontos,a nao ser nos proprios pontos da curva (respectivamente, superfıcie) onde a funcaoe nula. Estas componentes sem variacao de sinal nao podem ser detectadas pelosmetodos numericos tradicionais porque todos eles se baseiam no Teorema do ValorIntermedio, ou seja, todos assumem que a funcao muda de sinal de um lado para ooutro da curva.

Por conseguinte, um novo metodo numerico foi especialmente desenvolvido paradeterminar quer componentes sem variacao de sinal, quer componentes com variacaode sinal. Designou-se este metodo por metodo generalizado da falsa posicao (GFP).Basicamente, o processo de amostragem de cada componente de uma curva e com-binado com a particao binaria (BSP) dum subespaco rectangular Ω contido em R2,designado por espaco ambiente. O algoritmo BSP subdivide recursivamente o espacoambiente em duas partes atraves de um segmento de bisseccao. Depois, e determinadocada ponto de interseccao entre o segmento de bisseccao e a curva implıcita usando o

ix

x

metodo GFP. O algoritmo BSP foi projectado para representar nao so curvas implıcitasem R2, mas tambem curvas resultantes da interseccao de uma superfıcie implıcita comum plano em R3. Marcante e tambem o facto de que o algoritmo GFP poder ser usadopara calcular pontos isolados de curvas e superfıcie implıcitas atraves duma sequenciade mınimos absolutos que tendem para zero.

Abstract

Mathematical systems (e.g. Maple, MathLab, and Mathematica) can easily plotimplicit curves and surfaces defined implicitly by polynomial functions in R2 andin the space R3, respectively. However, none of these systems preserve the topologicalshape of general implicit curves and surfaces. In other words, they cannot correctlydisplay implicit curves and surfaces, in particular their singularities (e.g. cusps andself-intersections) and sign-invariant components, including isolated points.

In addition to isolated points and sign-invariant components of curves and sur-faces, current algorithms found in the literature cannot correctly generate and displayinhomogeneous implicit surfaces (e.g. Steiner surface and Cartan surface). For ex-ample, the negative z-axis belongs to the Cartan surface, but in practice no currentsystem or algorithm is capable of plotting it. In fact, the research focus in this field hasbeen on homogeneously 2-dimensional surfaces, regardless whether they are manifoldsor not.

Topologically speaking, a sign-invariant component of an implicit curve (respec-tively, surface) is a maximal connected subspace of the curve (respectively, surface)such that its describing real function evaluates always either positive or negative ina small neighborhood of each point of it, unless at each curve or surface point atwhich the function vanishes. These sign-invariant components cannot be detected bytraditional numerical methods because all them are based on the Intermediate ValueTheorem, i.e. they all assume that the function changes sign from one side to otherside of the curve.

As a consequence, a new numerical method was specially designed for samplingboth sign-invariant components and sign-variant components, called generalized falseposition (GFP) method. The sampling process of a curve is interspersed with the spacebinary partitioning (BSP) of its ambient space in R2. It recursively subdivides thespace in two parts by a bisection segment at each step. For each bisection segment, itcomputes the intersection with the implicit curve by using the GFP method. The BSPalgorithm can be used to plot any implicit curve in R2, as well as the any implicit curvethat results from intersecting an implicit surface and a plane in R3. Remarkably, theGFP method is also used to compute isolated points of implicit curves and surfaces.

xi

Publicacoes

Durante os trabalhos de doutoramento foram efectuadas as seguintes publicacoes:

- Jose Morgado and Abel Gomes. Rendering implicit curves with isolated pointsthrough binary space partitioning. Virtual - The Portuguese Journal of Com-puter Graphics, 2004.

- Francisco Morgado and Abel Gomes. A derivative-free tracking algorithm forimplicit curves with singularities. Proceedings of the 3rd International Work-shop on Computer Graphics and Geometric Modeling (CGGM’2004), also inM. Bubak, D. Albada, P. Sloot, and J. Dongarra (eds.), Proceedings of theInternational Conference on Computational Conference (ICCS’2004), Krakow,Poland, June 6-9, 2004, Lecture Notes in Computer Science, Vol.3039, Springer-Verlag, 2004.

- Francisco Morgado e Abel Gomes. Representacao BSP de curvas implıcitas 2D.Actas do 12o Encontro Portugues de Computacao Grafica (EPCG’2003), Porto,Portugal, Outubro 8-10, 2003.

- Amandio Marques and Francisco Morgado. Application of a genetic algorithmto a scheduling assignment problem soft computing and complex systems. CentroInternational de Matematica, Coimbra, Portugal, 2003.

- Francisco Morgado and Abel Gomes. A non-uniform binary space partition algo-rithm for 2D implicit curves. Proceedings of the 2nd International Workshop onComputer Graphics and Geometric Modeling (CGGM’2003), also in V. Kumar,M. Gavrilova, C. Tan, and P.L’Ecuyer (eds.), Proceedings of the InternationalConference on Computational Science and Its Applications (ICCSA’2003), Mon-treal, Canada, May 18-21, 2003, Lectures Notes in Computer Science, 3(2669),Springer-Verlag, 2003.

- Francisco Morgado and Pedro Almeida. Interactive Environment for TeachingDescritive Geometry. M. Nistal, M. Iglesias and L. Rifon (eds.), Computers andEducation Towards a Lifelong Learning Society, Kluwer Academic Publishers,pages 251–262, 2003.

- Francisco Morgado e Pedro Almeida. Ensino da Geometria Descritiva Assistidopelo Computador. M. Nistal, M. Iglesias and L. Rifon (eds.), Proceedings of the

xiii

xiv

4th International Symposium of Educative Computer Science (also 6o CongresoIberoAmericano de Informatica Educativa (IE’2002)), Vigo, Spain, November21-23, 2002.

- Francisco Morgado and Abel Gomes. A new 2D implicit curve-tracking algo-rithm. K. Wojciechowski (ed.), Proceedings of the International Conference onComputer Vision and Graphics (ICCVG’2002), Zakopane, Poland, September25-29, 2002.

- Francisco Morgado and Abel Gomes. Fast representation of implicit curvesthrough space subdivision. P. Santos, J. Jorge and A. Marcos (eds.), Proceedingsof the 1st Ibero-American Symposium in Computer Graphics (SIACG’2002),Guimaraes, Portugal, 2002.

Estrutura da Tese

A presente tese esta organizada em cinco capıtulos, nomeadamente:

- No Capıtulo 1 introduz-se os conceitos fundamentais acerca das funcoes reaisde variaveis reais. Estas funcoes sao usadas para descrever a geometria deobjectos geometricos, em particular curvas e superfıcies implıcitas. Faz-se aindaum levantamento do estado-da-arte dos algoritmos utilizados para amostrar epoligonizar estes objectos geometricos.

- No Capıtulo 2 introduz-se um novo metodo numerico e suas variantes paradeterminar zeros e extremos de funcoes reais de variavel real. O calculo deum zero ou de um extremo e feito atraves da mesma formula de interpolacao.Como se vera, isto e particularmente importante na amostragem de componentesde curvas e superfıcies que nao sao detectadas pelos algoritmos actuais.

- O Capıtulo 3 descreve um algoritmo de representacao de curvas implıcitas nao-homogeneas em que a amostragem e feita por subdivisao BSP.

- No Capıtulo 4 descreve-se um algoritmo de representacao de superfıcies implıcitasnao-homogeneas em que a amostragem e feita por subdivisao regular do espacoambiente em celulas cubicas.

- Por fim no Capıtulo 5 apresentam-se algumas conclusoes do trabalho realizado,assim como as suas contribuicoes para o domınio da computacao grafica. A tesetermina com a indicacao de alguns pontos deixados em aberto e que poderao serassunto de trabalho futuro.

O objectivo principal desta tese foi o de resolver o problema de representacaode curvas e superfıcies dimensionalmente degeneradas, assim como o problema depreservacao da forma topologica (furos, componentes e singularidades) destes objectosgeometricos durante a sua representacao grafica.

xv

Indice

Agradecimentos vii

Resumo ix

Abstract xi

Publicacoes xiii

Estrutura da Tese xv

Lista de Figuras xxi

Lista de Tabelas xxiii

Lista de Sımbolos e Siglas xxv

1 Curvas e Superfıcies Implıcitas: Estado-da-Arte 1

1.1 Conceitos basicos associados a funcoes reais . . . . . . . . . . . . . . . 2

1.1.1 Imagem de uma funcao . . . . . . . . . . . . . . . . . . . . . . . 2

1.1.2 Grafico de uma funcao . . . . . . . . . . . . . . . . . . . . . . . 2

1.1.3 Conjunto de nıvel de uma funcao . . . . . . . . . . . . . . . . . 3

1.2 Representacao de curvas implıcitas . . . . . . . . . . . . . . . . . . . . 3

1.2.1 Algoritmos baseados na continuacao da curva . . . . . . . . . . 4

1.2.2 Algoritmos baseados na conversao da representacao . . . . . . . 5

1.2.3 Metodos de subdivisao do espaco . . . . . . . . . . . . . . . . . 5

xvii

xviii INDICE

1.3 Representacao de superfıcies implıcitas . . . . . . . . . . . . . . . . . . 11

1.3.1 Definicao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.3.2 Algoritmos de representacao de superfıcies implıcitas . . . . . . 11

1.4 Problemas numericos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.5 Sumario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2 Computacao de zeros e extremos de funcoes 15

2.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.2 Metodo generalizado da falsa posicao . . . . . . . . . . . . . . . . . . . 16

2.2.1 Determinacao dos zeros . . . . . . . . . . . . . . . . . . . . . . . 17

2.2.2 Determinacao dos extremos . . . . . . . . . . . . . . . . . . . . 17

2.2.3 Analise da convergencia . . . . . . . . . . . . . . . . . . . . . . 18

2.2.4 Implementacao do metodo GFP . . . . . . . . . . . . . . . . . . 19

2.3 Resolucao das descontinuidades . . . . . . . . . . . . . . . . . . . . . . 21

2.4 Condicoes de paragem . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.5 Metodo TGFP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.5.1 Convergencia do metodo TGFP . . . . . . . . . . . . . . . . . . 26

2.6 Metodo NGFP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.7 Resultados experimentais . . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.7.1 Resultados experimentais dos metodos GFP, TGFP e NGFP . . 29

2.7.2 Comparacao com outros metodos . . . . . . . . . . . . . . . . . 33

2.8 Sumario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3 Curvas Implıcitas Nao-Homogeneas 35

3.1 Subdivisao do domınio . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.1.1 Estrutura de dados: arvore BSP . . . . . . . . . . . . . . . . . . 36

3.1.2 Linha de bisseccao . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.1.3 Determinacao dos subespacos Ω− e Ω+ . . . . . . . . . . . . . . 39

3.2 Algoritmo BSP e condicoes de paragem . . . . . . . . . . . . . . . . . . 41

INDICE xix

3.3 Determinacao dos pontos da curva . . . . . . . . . . . . . . . . . . . . 42

3.3.1 Pontos com variacao de sinal . . . . . . . . . . . . . . . . . . . . 44

3.3.2 Pontos sem variacao de sinal . . . . . . . . . . . . . . . . . . . . 47

3.3.3 Pontos: algoritmo geral . . . . . . . . . . . . . . . . . . . . . . . 50

3.4 Determinacao de pequenas componentes . . . . . . . . . . . . . . . . . 53

3.4.1 Pontos isolados . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3.4.2 Teste do ponto isolado . . . . . . . . . . . . . . . . . . . . . . . 56

3.4.3 Sobre-particao dos subespacos terminais . . . . . . . . . . . . . 57

3.5 Determinacao de pontos de auto-interseccao . . . . . . . . . . . . . . . 60

3.6 Visualizacao da Curva . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

3.7 Comparacao e resultados experimentais . . . . . . . . . . . . . . . . . . 63

3.7.1 Comparacao geral . . . . . . . . . . . . . . . . . . . . . . . . . . 63

3.7.2 Resultados experimentais . . . . . . . . . . . . . . . . . . . . . . 68

3.8 Sumario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

4 Superfıcies Implıcitas Nao-Homogeneas 75

4.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

4.1.1 Algoritmos baseados na subdivisao do espaco . . . . . . . . . . 75

4.1.2 Algoritmos de visualizacao directa da superfıcie . . . . . . . . . 77

4.2 Subdivisao nao-adaptativa: estrutura de dados . . . . . . . . . . . . . . 77

4.3 O algoritmo ImplicitSurface . . . . . . . . . . . . . . . . . . . . . . . 79

4.4 Poligonizacao das celulas . . . . . . . . . . . . . . . . . . . . . . . . . . 80

4.4.1 Celulas do tipo 0: pontos isolados . . . . . . . . . . . . . . . . . 82

4.4.2 Celulas do tipo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 84

4.4.3 Celulas do tipo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 86

4.4.4 Celulas do tipo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.4.5 Celulas do tipo 4 . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.4.6 Celulas do tipo 5 . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.5 Triangulacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

xx INDICE

4.6 Resultados experimentais . . . . . . . . . . . . . . . . . . . . . . . . . . 90

4.6.1 Apreciacao geral . . . . . . . . . . . . . . . . . . . . . . . . . . 90

4.6.2 Desempenho do algoritmo . . . . . . . . . . . . . . . . . . . . . 92

4.7 Sumario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

5 Conclusoes, Contribuicoes e Trabalho Futuro 95

5.1 Conclusoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

5.2 Contribuicoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

5.3 Trabalho futuro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

Bibliografia 99

Lista de Figuras

1.1 Subconjuntos definidos por uma funcao f . . . . . . . . . . . . . . . . . 2

1.2 Algoritmos baseados na continuacao da curva . . . . . . . . . . . . . . 4

1.3 Subdivisao quadtree. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.4 Subdivisao BSP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.5 Subdivisao usando aritmetica de intervalo. . . . . . . . . . . . . . . . . 9

1.6 (a) Funcao com variacao de sinal (b)(c) Funcao com potencial positivoe negativo (d) Funcao com e sem variacao de sinal. . . . . . . . . . . . 13

2.1 Ilustracao do metodo generalizado da falsa posicao. . . . . . . . . . . . 16

2.2 Assimptota vertical . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.3 Metodo TGFP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.4 Metodo NGFP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.5 Graficos de funcoes reais com zeros e extremos. . . . . . . . . . . . . . 32

3.1 Subdivisao e arvore BSP. . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.2 Determinacao da linha BSP para tres subespacos. . . . . . . . . . . . . 39

3.3 Classificacao de pontos em relacao a linha BSP. . . . . . . . . . . . . . 41

3.4 Curvas com componentes de tipos distintos. . . . . . . . . . . . . . . . 46

3.5 Curvas sem variacao de sinal. . . . . . . . . . . . . . . . . . . . . . . . 49

3.6 Possibilidade de haver mais zeros em XiXi+1. . . . . . . . . . . . . . . 52

3.7 Existencia de outro zero no interior de XiXi+1 quando f(Xi).f(Xi+1) = 0. 52

3.8 Determinacao de um ponto isolado num subespaco Ω. . . . . . . . . . . 55

xxi

xxii LISTA DE FIGURAS

3.9 Interseccao da curva com o quadrado inscrito na vizinhanca circular deP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.10 Tres casos onde e necessario particionar o subespaco. . . . . . . . . . . 58

3.11 Curvas algebricas: (esquerda) processo de subdivisao BSP e (a direita)poligonizacao. Regioes verdes e roxas e onde f > 0 e f < 0, respecti-vamente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

3.12 Curvas racionais: (esquerda) processo de subdivisao BSP e (a direita)poligonizacao. Regioes verdes e roxas e onde f > 0 e f < 0, respecti-vamente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

3.13 Curvas transcendentes. . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

3.14 Esquemas de subdivisao de curvas racionais sem assimptotas: (esquerda)BSP nao-uniforme, (centro) BSP uniforme, (direita) quadtree uniformecom AI. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

3.15 Comparacao dos tres metodos: BSP, BSU, Quadtree com AI. . . . . . . 73

4.1 Criacao da celula C(i, j, k). . . . . . . . . . . . . . . . . . . . . . . . . . 80

4.2 Varios tipos de celulas. . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

4.3 Processo de determinacao do ponto isolado. . . . . . . . . . . . . . . . 82

4.4 Celulas do tipo 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

4.5 Determinacao dos pontos de um segmento 1-dimensional de uma su-perfıcie. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

4.6 Ponto de uniao entre um fragmento 2-dimensional e um segmento 1-dimensional de uma superfıcie numa celula do tipo 4. . . . . . . . . . . 88

4.7 Triangulacao de um fragmento 2-dimensional de uma superfıcie. . . . . 89

4.8 Superfıcies algebricas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

4.9 Superfıcies algebricas e racionais. . . . . . . . . . . . . . . . . . . . . . 91

Lista de Tabelas

2.1 f(x) = x2 + 5 em]− 1, 8[, mınimo local . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.2 f(x) = sin2(x) + 1 em ]1, 1.8[, maximo local . . . . . . . . . . . . . . . 30

2.3 f(x) = x3

6− x2

2+ x + 3 em ]− 3.0, 2.0[, zero secante . . . . . . . . . . . 31

2.4 f(x) = tg(x)tg(x) − 103 em ]1.3, 1.4[, zero secante . . . . . . . . . . . . . 31

2.5 f(x) = x ex − 10 em ]− 10.0, 10.0[ zero secante . . . . . . . . . . . . . . 31

2.6 f(x) = |x|1/3 em]− 0.5, 1.5[, zero extremo . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.7 Comparacao entre metodos numericos em termos do numero Nf deestimativas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.1 y3 − x2 + 2xy = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

3.2 y3 − x2 + 2xy − x = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

3.3 (x− 1y2+1

)(y − 1x2+1

) = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . 71

3.4 0.004+0.11x−0.177y−0.174x2+0.224xy−0.303y2−0.168x3+0.327x2y−0.087xy2−0.013y3+0.235x4−0.667x3y+0.745x2y2−0.029xy3+0.072y4 =0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

4.1 Tempos de processamento e de visualizacao de superfıcies (em segundos). 92

xxiii

Lista de Sımbolos e Siglas

f - Funcao de Rn em R

CDf - Contradomınio de f

Cf - Conjunto que contem o contradomınio de f

Ω - Subespaco de R2 ou R3

Ω+ - Semi-espaco de Ω a direita da linha de bisseccao

Ω− - Semi-espaco de Ω a esquerda da linha de bisseccao

Ω - Subespaco de R2 ou R3

d(x) ≈ f(x+dx)−f(x)dx

, com dx ≈ 0

CSG - Constructive Solid Geometry

CAD - Computed Aided Design

BSP - Binary Space Partition

GFP - Generalized False Position

TGFP - Tangent-extended Generalized False Position

NGFP - Newton-extended Generalized False Position

TCP/IP - Transmission Control Protocol / Internet Protocol

∇f - Gradiente de f

xxv

Capıtulo 1

Curvas e Superfıcies Implıcitas:Estado-da-Arte

Um dos problemas classicos em computacao grafica e a representacao de curvas esuperfıcies definidas implicitamente por funcoes reais. A importancia destes objectosgeometricos advem de varios factos, entre os quais:

• Permitem uma descricao concisa de objectos geometricos, que e util pela econo-mia de memoria que representa e pela transmissao rapida atraves do canalTCP/IP;

• Sao manipulaveis atraves de operacoes de blending, operacoes de deformacao e,em geral, utilizaveis em animacao [22, 39, 45].

• Facilidade na aplicacao de algoritmos de ray casting e ray tracing [23, 59] nageracao de imagens realısticas de cenas, incluindo a representacao de objectostranslucidos, efeitos de sombras,etc.

• Podem ser usados na definicao de primitivas de modeladores geometricos. Nummodelador CSG (Constructive Solid Geometry), cada primitiva geometrica edefinida pela combinacao booleana de semi-espacos, sendo que cada semi-espacoe, normalmente, descrito por uma funcao real de tres variaveis reais [42, 43]. Porexemplo, um cilindro e o resultado da interseccao de tres semi-espacos definidospor funcoes implıcitas.

Por forma a compreender a natureza e as propriedades das curvas e superfıciesimplıcitas, este capıtulo comeca por abordar os conceitos matematicos que lhes estaoassociadas. Segue-se uma descricao sumaria dos principais algoritmos utilizados narepresentacao destas curvas e superfıcies. A ideia e, pois, mostrar as vantagens elimitacoes de tais algoritmos por forma a contextualizar o trabalho de investigacaoque conduziu a elaboracao e escrita dos Capıtulos 2, 3 e 4 da presente tese.

1

2 CAPITULO 1. CURVAS E SUPERFICIES IMPLICITAS: ESTADO-DA-ARTE

m

n

Image f f

x f(x)

c f -1(c)

Figura 1.1: Subconjuntos definidos por uma funcao f .

1.1 Conceitos basicos associados a funcoes reais

As funcoes reais sao utilizadas para representar os mais diversos objectos geometricosem computacao grafica. A Figura 1.1 ilustra os tres tipos principais de objectos asso-ciados a uma funcao: imagem, grafico e conjunto de nıvel. As definicoes matematicasdestes objectos podem encontrar-se em varios livros. Na presente tese adoptou-se asdefinicoes usadas por Baxandall [8].

1.1.1 Imagem de uma funcao

Definicao 1.1 Seja Ω um subconjunto de Rm. A imagem da funcao f : Ω ⊂ Rm → Rn

e o subconjunto de Rn dado por Image f = y ∈ Rn | y = f(x),∀x ∈ Ω.

Por exemplo, seja f : R → R2 a funcao definida por f(t) = (x(t), y(t)) =(cos t, sin t), t ∈ R. A imagem de f e a circunferencia x2 + y2 = 1 de centro (0, 0) eraio 1 em R2, uma vez que x(t)2 + y(t)2 = cos2 t + sin2 t = 1.

Os objectos geometricos gerados atraves de imagens de funcoes sao conhecidosem computacao grafica por curvas e superfıcies parametricas. As curvas e superfıciesde Bezier [12] e as B-splines [16] sao exemplos bem conhecidos deste tipo de objectos.

1.1.2 Grafico de uma funcao

Definicao 1.2 Seja Ω um subconjunto de Rm. O grafico da funcao f : Ω ⊂ Rm → Rn

e o subconjunto do espaco produto Rm+n = Rm × Rn definido por Graph f = (x,y) |x ∈ Ω e y = f(x)

Por exemplo, o grafico da funcao f(x) = x2 e o conjunto de todos os pontos(x, f(x)) = (x, x2), ou seja uma parabola. Repare-se que o grafico de f e representadono espaco produto R1+1 = R1 × R1.

1.2. REPRESENTACAO DE CURVAS IMPLICITAS 3

Objectos geometricos gerados por graficos de funcoes sao casos particulares deconjuntos de nıvel. De facto, prova-se que todo o grafico duma funcao e um conjuntode nıvel [8]. Como se vera adiante, curvas e superfıcies implıcitas sao conjuntos denıvel de funcoes reais.

1.1.3 Conjunto de nıvel de uma funcao

Definicao 1.3 Seja Ω um subconjunto de Rm, c = (c1, ..., cn) um ponto de Rn ef : Ω ⊂ Rm → Rn. O conjunto de nıvel de f , denotado por f−1(c) e entao definidopor f−1(c) = x ∈ Ω | f(x) = c

Na definicao anterior, a funcao f e denominada em geral por campo escalar oupotencial da funcao, dependendo do contexto em que se da a utilizacao de f . Nestatese, optou-se por usar a designacao de potencial da funcao. Assim, determinar f−1(c)equivale a determinar os pontos x ∈ Ω onde o potencial e c. Note-se que f−1 naodeve ser confundido com a funcao inversa de f , pois podera acontecer que f nao tenhafuncao inversa.

Daqui em diante, assumir-se-a que c = 0, logo determinar f−1(0) —denominadopor conjunto de nıvel 0 ou nucleo Nf de f— equivale a determinar os pontos de Ωonde o potencial e nulo, i.e. Nf = x ∈ Ω : f(x) = 0. Se n = 1 e m = 2, entao Nf

representa uma curva ou um ponto. Por exemplo, o nucleo Nf de f(x, y) = x2 +y2−1e uma circunferencia de centro (0, 0) e raio 1. Se n = 1 e m = 3, entao Nf representauma superfıcie, uma curva ou um ponto. Por exemplo, o nucleo Nf de f(x, y, z) =x2 + y2 + z2, e o ponto (0, 0, 0).

E importante referir que a definicao de Nf nao permite determinar directamenteos pontos de um objecto, como no caso das funcoes parametricas, onde para cada valordo(s) parametro(s) se obtem o ponto do objecto atraves do calculo de f . Por outrolado a determinacao de Nf envolve a resolucao (numerica) da equacao f(x) = 0, ou

seja pretende-se determinar os pontos cujos transformados por f e nulo. E por estefacto que as curvas e superfıcies parametricas sao mais faceis de visualizar do que ascurvas e superfıcies implıcitas.

Em resumo, pode dizer-se que um conjunto de nıvel e um subconjunto do domınioda funcao, a imagem e o contradomınio da funcao, e o grafico e um subconjunto doespaco produto.

1.2 Representacao de curvas implıcitas

Em computacao grafica, ha tres famılias principais de algoritmos para representarcurvas implıcitas: algoritmos baseados na continuacao da curva, algoritmos baseadosna conversao da representacao e algoritmos baseados na subdivisao do espaco.


(b)

0

P

1

2

3 4 5

6

7

Q

(a) (c)

f(P)

Q P

X

Q P

Figura 1.2: Algoritmos baseados na continuacao da curva

1.2.1 Algoritmos baseados na continuacao da curva

A ideia base destes algoritmos consiste em determinar o proximo ponto Q da curvaa custa do ponto corrente P . A determinacao do proximo ponto pode ser feitono espaco da imagem pixelizada ou no espaco da curva (domınio da funcao). Naliteratura, e possıvel encontrar varios algoritmos para calcular o proximo ponto dacurva, nomeadamente:

• Procura no espaco da imagem. Em vez de procurar o proximo ponto, procura-seo pixel no espaco discreto da imagem que mais perto esta do proximo ponto Q[13]. Basicamente, calcula-se o valor da funcao nos 8 pixeis vizinhos do pixelactual P , sendo Q o pixel vizinho mais proximo da curva (Figura 1.2(a));

• Procura no espaco do objecto. Ha duas formas conhecidas de calcular o proximoponto Q no espaco do objecto. Talvez a mais conhecida seja a que usa o correctorde Newton [4], como se ilustra na Figura 1.2(b). A partir do ponto P , calcula-seo ponto X atraves de um vector ortogonal ao gradiente de f em P . De seguida,aplica-se o corrector de Newton a X de modo a obter o ponto Q da curva.Embora rapido, este algoritmo usa a derivada da funcao, o que dificulta a suaaplicacao a curvas definidas por funcoes nao-diferenciaveis [34].

Outra forma de calcular o proximo ponto Q e atraves da interseccao da curvaC com uma circunferencia centrada no ponto corrente P [36]. O ponto Q podeser calculado atraves de um metodo numerico desprovido de derivadas (e.g. ometodo da falsa posicao) ao longo de um arco de circunferencia (Figura 1.2(c)).

Os algoritmos de continuacao sao atractivos, pois concentram o esforco de procurae computacao do proximo ponto da curva onde existe maior probabilidade de eleexistir. Isto significa que estes algoritmos se adaptam a geometria local da curva. Noentanto, os algoritmos que fazem a procura no espaco da curva adaptam-se melhora geometria local da curva, pois sao calculados mais pontos onde a curvatura e maiselevada e menos pontos onde a curva tem curvatura mais reduzida.


A desvantagem destes algoritmos e que precisam do conhecimento previo dumponto (inicial) em cada componente da curva, o que acarreta um pre-processamentoao algoritmo. Alem disso, a determinacao dos pontos isolados —que sao componentessingulares— duma curva e uma tarefa para o qual nao existe algoritmo numericoconhecido.

1.2.2 Algoritmos baseados na conversao da representacao

A ideia base destes algoritmos e a de converter a representacao implıcita duma curvanuma representacao parametrica, sendo depois aplicado um algoritmo de representacaode curvas parametricas [56, 6, 27]. Por exemplo, uma representacao parametrica dacurva x3 + y3 = 3axy, denominada de folium de Descartes, e a seguinte: x(t) = 3at

1+t3,

y(t) = 3at2

1+t3. Uma representacao parametrica de uma circunferencia e x(t) = cos(t),

y(t) = sin(t), t ∈ [0, 2π], mas ha outras possıveis como, por exemplo, a seguinte:x(t) = cos(2t), y(t) = sin(2t), t ∈ [0, π]. Isto significa que a representacao parametricade uma curva nao e unica [8].

Contudo, nem sempre e possıvel determinar uma parametrizacao global da curva[15], a nao ser que a curva seja regular, i.e. ∇f 6= 0 em todos os seus pontos. Estefacto inviabiliza a construcao de um algoritmo generico para a representacao de curvasimplıcitas atraves da formulacao parametrica. E por isto que Abhyankar e Bajajsugerem algoritmos de conversao especıficos para cada tipo de curva [1, 3, 2].

Recorde-se, no entanto, que a conversao duma representacao parametrica narepresentacao implıcita duma curva e sempre possıvel, sendo este processo conhecidopor implicitacao, a qual se baseia na teoria de eliminacao [28]. Contudo, este processode eliminacao gera funcoes implıcitas de grau elevado. E por isso que varios autorestem proposto metodos especıficos de implicitacao para diferentes classes de funcoes,em particular para funcoes parametricas racionais [31, 46].

1.2.3 Metodos de subdivisao do espaco

O princıpio destes algoritmos e o de dividir para conquistar. Sao metodos que sub-dividem o espaco inicial em subespacos mais pequenos ate que um dado criterio deparagem se verifique [9, 55]. Os metodos de subdivisao exigem recursos significativosdo computador, em particular um maior tempo de processamento e mais memoria[48]. No entanto, com o desempenho crescente dos processadores modernos e daevolucao das arquitecturas dos computadores, os metodos baseados na subdivisaotem-se tornado cada vez mais atractivos.

No ambito da representacao de curvas implıcitas em R2, ha duas formas desubdivisao do espaco 2-dimensional: quadtrees [48] e BSP trees [37]. Como se mostranas Figuras 1.3 e 1.4 respectivamente, a principal diferenca entre elas e que, na


(a) Região inicial e a árvore

(b) Primeira partição

(c) Segunda partição

1 2

3 4

2 1 3 4

2 1 3 4

1

. . .

Figura 1.3: Subdivisao quadtree.

subdivisao quadtree o espaco e dividido sempre em 4 partes, podendo ser iguais ounao, ao passo que na subdivisao BSP o espaco e dividido em duas partes.

Duas das tecnicas mais conhecidas de controlo da subdivisao quadtree sao:aritmetica intervalar [48] e aritmetica afim [19, 32].

Aritmetica Intervalar (AI). A aritmetica intervalar foi introduzida por Moore[35] no contexto de investigacao matematica, tendo sido recentemente aplicada narepresentacao de curvas implıcitas algebricas [48]. A aritmetica intervalar consistenum conjunto de operacoes simples sobre intervalos, nomeadamente:

(1) Adicao: [a, b] + [c, d] = [a + c, b + d]

(2) Subtraccao: [a, b]− [c, d] = [a− d, b− c]

(3) Produto: [a, b]× [c, d] = [min(ac, ad, bc, bd), max(ac, ad, bc, bd)]

(4) Divisao: [a, b]/[c, d] = [a, b]× [1/d, 1/c]; considera-se que 0 /∈ [c, d]

A questao que se coloca e saber entre que valores varia a funcao f para qualquerponto (x, y) ∈ Ω, o que e equivalente a determinar o contradomınio CDf de f . Por


l1 +

-

l2 l1

l3

l1

+ -




l1

l3 l2

Figura 1.4: Subdivisao BSP.


exemplo, considere-se a curva definida pela funcao f(x, y) = xy + y no espaco Ω =(x, y) ∈ R2 : −1 ≤ x ≤ 1 ∧ 2 ≤ y ≤ 3. Se 0 /∈ CDf , entao nao existe qualquersegmento da curva em Ω.

Embora nao se conheca nenhum algoritmo que permita determinar com exac-tidao o contradomınio de f , a aritmetica intervalar permite determinar um conjuntoCf que contem todos os valores de f , i.e Cf ⊃ CDf . O processo de determinacao deCf e feito do seguinte modo:

(x, y) ∈ Ω

x ∈ [−1, 1] e y ∈ [2, 3]

x× y ∈ [−3, 3]

x× y + y ∈ [−1, 6]

f(x, y) ∈ [−1, 6]

Recorde-se que uma curva implıcita e definida por Nf = (x, y) ∈ Ω : f(x, y) =0. Tendo em conta que 0 ∈ [−1, 6], entao existe grande probabilidade de a curvaestar parcial ou totalmente contida em Ω. Assim sendo, subdivide-se Ω em 4 novosespacos (Ω1,Ω2,Ω3,Ω4), sendo de seguida repetido o processo para cada um dos novosespacos criados. O processo termina quando algum criterio de paragem for verificado.

Note-se que se a funcao anterior for apresentada na forma f(x, y) = y(x + 1)obtem-se:

(x, y) ∈ Ω

x ∈ [−1, 1] e y ∈ [2, 3]

x + 1 ∈ [0, 2]

y × (x + 1) ∈ [0, 6]

f(x, y) ∈ [0, 6]

Logo em Ω podera existir algum segmento da curva pois 0 ∈ [0, 6]. Isto e, paradiferentes decomposicoes de f , pode-se obter diferentes intervalos Cf onde esta contidoo seu contradomınio CDf .

A analise intervalar e, pois, um metodo que permite determinar um conjuntoCf que contem sempre o contradomınio CDf da funcao que se esta a representar.O facto de se determinar um conjunto Cf que pode nao ser totalmente coincidentecom o contradomınio da funcao leva a que, por vezes, sejam efectuadas subdivisoesdesnecessarias em zonas de Ω em que a curva nao passa, como acontece no cantoinferior esquerdo da Figura 1.5(a). Estas subdivisoes redundantes sao tanto maisquanto maior for a diferenca entre Cf e CDf , fenomeno este que e mais acentuado emfuncoes descritas por expressoes mais extensas. Pode dizer-se que o ideal seria ter ummetodo que permitisse determinar CDf com exactidao, pois assim saber-se-ia logo seera necessario efectuar mais subdivisoes no subespaco Ω ou nao.


(a) 0.004 + 0.11x− 0.177y − 0.174x2 (b) y − 1x2 = 0

+0.224xy − 0.303y2 − 0.168x3 + 0.327x2y−0.087xy2 − 0.013y3 + 0.235x4 − 0.667x3y

+0.745x2y2 − 0.029xy3 + 0.072y4 = 0

Figura 1.5: Subdivisao usando aritmetica de intervalo.

Aritmetica Afim (AA). A aritmetica afim e uma tecnica complementar a AI quepermite minimizar a diferenca entre Cf e CDf [19]. De facto, quanto mais pequeno forCf , maior probabilidade existe em que 0 /∈ Cf , logo e mais provavel que este subespaconao seja mais subdividido.

Enquanto que na AI uma variavel assume valores de um intervalo —como eusual—, na AA uma variavel e ainda representada da seguinte forma:

x = x0 + x1ε1 + x2ε2 + ... + xnεn

onde xi sao coeficientes reais e εi sao variaveis que assumem valores no intervalo [−1, 1].Cada εi (noise symbol) representa uma percentagem (negativa ou positiva) de xi paraa variacao total de x.

Como exemplo, seja x ∈ [a, b] em AI. A forma correspondente em AA e dadapor x=x0 + x1ε1, onde

x0 =b + a

2e x1 =

b− a

2

Em geral, se x=x0 + x1ε1 + x2ε2 + ... + xnεn em AA, entao o correspondenteintervalo para a variavel x em AI e o seguinte:

[x0 − ε, x0 + ε] com ε =n∑

i=1

|xi|

Sejam, x e y duas variaveis em AA dado por:

x = x0 + x1ε1 + ... + xnεn


y = y0 + y1ε1 + ... + ynεn

e α ∈ R um escalar real. Entao, as operacoes basicas em AA sao:

(1) Adicao/Subtraccao: x± y=(x0 ± y0) + (x1 ± y1)ε1 + ... + (xn ± yn)εn

(2) Multiplicacao por escalar: αx=(αx0) + (αx1)ε1 + ... + (αxn)εn

(3) Adicao/Subtraccao com escalar: x± α=(x0 ± α) + x1ε1 + ... + xnεn

A operacao de multiplicacao entre x e y e dada por z = xy, onde

z0 = x0y0

zi = x0yi + y0xi, (i = 1, ..., n)

zk = uv

com u =∑n

i=1 |xi|, v =∑n

i=1 |yi| e z=z0 + z1ε1 + ...+ znεn + zkεk. O termo extra zkεk,representa o erro introduzido pela aproximacao.

Como exemplo, considere-se a funcao f(x) = x(10− x) com x ∈ [4, 6]. Assim:

x = 5 + 1ε1

10− x = 5− 1ε1

x× (10− x) = 25 + 0ε1 − 1ε2

donde se conclui que, por utilizacao da AA, o intervalo ao qual f pertence e [25 −1, 25+1] = [24, 26]. Se se usar as operacoes definidas em AI, o intervalo para a funcaof sera [16, 36]. Note-se, no entanto, que o contradomınio de f e o intervalo [24, 25].

Este exemplo mostra que a AA permite encontrar um intervalo que se aproximamais do contradomınio da funcao do que o metodo AI, e assim reduzir o numero desubdivisoes desnecessarias em relacao ao metodo AI. Embora os calculos efectuadospela AA tendam a ser mais dispendiosos em termos de tempo do que na AI, a suamaior precisao pode fazer com que o algoritmo tenha convergencia mais rapida paraa curva.

Uma comparacao mais detalhada entre os varios metodos de representacao decurvas implıcitas baseados na aritmetica intervalar foi levada a cabo por Martin [32].Alguns destes metodos baseiam-se nas operacoes da AI e da AA, usando diferentes de-composicoes da funcao. Estas decomposicoes podem ser efectuadas de varias maneiras:potencias basicas, metodo de Horner, polinomios de Bernstein, metodo de Taubin,metodo de Rivlin e, ainda, metodo de Gopalsamy.

No entanto estes metodos nao sao validos para qualquer tipo de funcao. Porexemplo, considere-se a curva f(x, y) = y− 1

x2 da Figura 1.5(b) que tem uma assimp-tota vertical em x = 0. Como se mostra, os metodos AI e AA efectuam sucessivas

1.3. REPRESENTACAO DE SUPERFICIES IMPLICITAS 11

subdivisoes na vizinhanca tubular da assimptota. Isto deve-se ao facto de que aesquerda e a direita da assimptota a funcao tem valores de sinais contrarios, o queimplica que a assimptota e interpretada como pertencente a curva, i.e. f(x, y) = 0,o que nao e verdade. Uma assimptota vertical ocorre quando temos uma divisao por0, ou quando o valor da funcao numa pequena vizinhanca da assimptota tende para±∞.

1.3 Representacao de superfıcies implıcitas

1.3.1 Definicao

Segundo Bloomenthal [10], uma superfıcie implıcita e uma entidade geometrica 2-dimensional imersa no espaco tridimensional. Dentre as superfıcies implıcitas maisconhecidas encontram-se as quadricas, que sao um caso particular das superfıciesalgebricas representadas por funcoes polinomiais de grau 2 em R3. Superfıcies esfericas,cilindros, elipsoides sao quadricas usadas num grande numero de modeladores geome-tricos e sistemas CAD, ate porque existe um grande numero de algoritmos especial-izados na interseccao entre quadricas [47].

As super-quadricas sao uma generalizacao das quadricas, pois os expoentes dosmonomios podem ser diferentes de 2 [7]. De uma forma geral uma super-quadrica tema forma:

f(x, y, z) = (αx)n + (βy)n + (γz)n = K

Mas ha tambem superfıcies implıcitas definidas por funcoes reais mais complexas como,por exemplo, as definidas por funcoes exponenciais, logarıtmicas, trigonometricas, etc.

1.3.2 Algoritmos de representacao de superfıcies implıcitas

De um modo geral, ha tres metodos para representar superfıcies implıcitas, nomeada-mente: emissao de raios (ray casting), distribuicao de partıculas e marcha de cubos(marching cubes).

• Emissao de raios. Este metodo de representacao introduzido por Appel [5] ebaseado na emissao de raios que tentam intersectar o objecto. Deste modo,pretende-se determinar os pontos de interseccao entre os raios e os objectos.Nao e muito utilizada na visualizacao em tempo real, uma vez que a maioria dohardware grafico esta optimizado para o processamento de superfıcies triangu-lares [59, 23].

• Distribuicao de partıculas. Esta representacao baseia-se na distribuicao de par-tıculas sobre a superfıcie do objecto. Essas partıculas em geral sao triangulos


ou cırculos, orientados segundo a normal a superfıcie do objecto. O processode distribuicao das partıculas na superfıcie segue um modelo baseado em fısica,onde forcas de atraccao sao impostas para manter as partıculas sobre a superfıcie,enquanto que forcas de repulsao sao usadas para manter as partıculas com umdeterminado afastamento umas das outras, garantindo assim uma distribuicaouniforme. As forcas de atraccao em geral sao calculadas em funcao do gradientee do sinal de f [18, 58]. Embora esta solucao seja rapida nao produz resultadosmuito satisfatorios porque nem sempre as partıculas aproximam adequadamentea superfıcie, sendo alem disso necessario cuidar da conectividade das partıculaspor forma a bem formar a superfıcie.

• Marcha de cubos. E um dos metodos mais conhecidos na poligonizacao desuperfıcies implıcitas [30]. Este metodo e baseado na subdivisao do espaco emcelulas cubicas. A estrutura de dados correspondente e conhecida por octree.Cada celula e classificada em relacao a superfıcie que se esta a representaratraves do calculo do valor de f nos vertices da celula [49]. Se em dois verticesde uma aresta da celula, f tem sinais opostos entao considera-se que existealgum fragmento da superfıcie entre estes vertices. Procede-se entao ao calculodo(s) ponto(s) da interseccao da referida aresta com a superfıcie. O factodo metodo marching cubes subdividir o espaco em pequenas celulas cubicasindependentemente de elas conterem algum fragmento da superfıcie, torna oalgoritmo mais lento e exigente em termos de espaco de armazenamento.

• Subdivisao adaptativa. Desbrun [17] e Velho [57] propuseram uma subdivisaoespacial adaptativa de modo que ocorra mais subdivisoes do espaco onde asuperfıcie tenha mais variacoes de forma. Normalmente, esta subdivisao adap-tativa do espaco e representada atraves de octrees, estando a subdivisao de cadasubespaco dependente da existencia de algum fragmento da superfıcie que lheseja transversal. Em geral, o teste de existencia e feito utilizando aritmeticaintervalar [51, 19].

1.4 Problemas numericos na representacao de cur-

vas e superfıcies implıcitas

Como o objectivo e determinar o objecto definido por Nf = x ∈ Ω : f(x) = 0 ouseja, determinar os pontos x ∈ Ω tal que f(x) = 0, entao a funcao f pode ser vistacomo uma funcao que permite classificar os pontos do domınio de f relativamenteao objecto Nf . Deste modo e possıvel saber, em princıpio, se um ponto x pertenceou nao ao objecto, atraves do valor de f(x). Se f(x) = 0, entao x pertence aoobjecto. Se f(x) < 0, diz-se que o ponto esta aquem (dentro) do objecto; no caso def(x) > 0 o ponto esta alem (fora) do objecto (Figura 1.6(a)). Na maior parte dostrabalhos, toma-se esta classificacao de pontos como correcta, pois assim os algoritmosnumericos actuais para a determinacao de raızes podem ser aplicados. Isto significa

1.4. PROBLEMAS NUMERICOS 13

que se a funcao f variar de sinal, entre dois pontos do espaco, entao atraves dosalgoritmos numericos para a determinacao de zeros e possıvel determinar a solucao def(x) = 0, desde que f seja uma funcao contınua e verifique as condicoes do Corolariodo Teorema de Bolzano.

(a) x2 + y2 − 1 = 0 (b) (x2 + y2 − 1)2 = 0

(c) −(x2 + y2 − 1)2 = 0 (d) (x2 + y2 − 1)2..((x− 2)2 + (y − 2)2 − 1) = 0

Figura 1.6: (a) Funcao com variacao de sinal (b)(c) Funcao com potencial positivo enegativo (d) Funcao com e sem variacao de sinal.

No entanto esta classificacao nem sempre e valida, pois pode acontecer que umponto esteja aquem ou alem do objecto e f(x) mantenha sempre o mesmo sinal. Porexemplo em f(x, y) = (x2 + y2 − 1)2 (Figura 1.6(b)), o Nf e uma circunferenciacentrada em (0, 0) e raio 1. Repare-se que, por exemplo, o ponto (0, 0) esta dentroda circunferencia no entanto f(0, 0) > 0, o ponto (2, 2) esta fora e f(2, 2) > 0, ouseja o facto de f(x) ser positivo ou negativo podera nao indicar se o ponto x estadentro ou fora do objecto. Repare-se que entre os pontos interiores e exteriores acurva representada na Figura 1.6(c) nao existe variacao de sinal de f . Isto significa


que os algoritmos numericos classicos de bisseccao, falsa posicao, secante e Newtonnao podem ser aplicados, pois nao e satisfeita a pre-condicao de variacao sinal de fem torno de cada zero da funcao. Logo, a utilizacao destes metodos numericos tornaimpossıvel a representacao deste tipo de curvas.

Outro exemplo e dado pela curva definida pela funcao f(x, y) = (x2 + y2 −1)2((x − 2)2 + (y − 2)2 − 1) = 0 (Figura 1.6(d)). Esta curva tem duas componentes.A primeira componente (x2 + y2 − 1)2 = 0 nao tem variacao de sinal, ao passo que asegunda componente (x− 2)2 +(y− 2)2− 1 = 0 tem variacao de sinal. Logo, qualqueralgoritmo que se apoie nos algoritmos numericos descritos nao podera representar aprimeira componente.

Um dos objectivos do trabalho de investigacao que conduziu a elaboracao destatese foi, como se vera nos capıtulos subsequentes, poder representar curvas e superfıciesimplıcitas na sua totalidade, quer tenham componentes com ou sem varicao de sinal.

1.5 Sumario

Com este capıtulo pretendeu-se fazer um levantamento sumario do estado-da-artedos algoritmos conhecidos de representacao de curvas e superfıcies implıcitas. Foramtambem apresentadas breves descricoes destes algoritmos por forma a evidenciar assuas caracterısticas, assim como as suas capacidades e limitacoes.

Uma das caracterısticas destes algoritmos reside no facto de serem desenvolvi-dos para tipos particulares de funcoes, nao sendo metodos gerais aplicaveis a qual-quer funcao real. Condicoes como a diferenciabilidade, continuidade, preservacaotopologica, e variacao de sinal das funcoes reais que descrevem curvas e superfıciesimplıcitas tem levantado bastantes dificuldades no desenvolvimento de algoritmosgenericos de representacao geometrica de curvas e superfıcies.

Por ser objectivo desta tese encontrar algoritmos gerais de representacao de cur-vas e superfıcies implıcitas, tornou-se inevitavel procurar, desenvolver e implementarum metodo numerico capaz de lidar com as dificuldades acima mencionadas, como sedescreve no Capıtulo 2.

Capıtulo 2

Computacao de Zeros e Extremosde Funcoes Reais

Este capıtulo surgiu da constatacao das limitacoes dos metodos numericos actuais naamostragem de curvas e superfıcies implıcitas. A amostragem e feita pela determinacaodos zeros da funcao associada a uma curva (respectivamente, superfıcie) ao longo decada segmento de recta (respectivamente, plano) que particiona o subespaco de R2

(respectivamente, R3) onde se encontra a curva (respectivamente, superfıcie).

Os algoritmos numericos actuais so determinam os zeros de uma funcao se forsatisfeita a pre-condicao de variacao de sinal da funcao (Corolario do Teorema deBolzano) no intervalo onde se procura o zero.

A contribuicao fundamental deste capıtulo e a de propor um novo metodonumerico de determinacao dos zeros de uma funcao que nao tem as pre-condicoeshabituais, entre as quais a da variacao de sinal. Este metodo e uma generalizacaodo metodo da falsa posicao (FP), pois a mesma formula de interpolacao e usada naoso para determinar os zeros, mas tambem os extremos (mınimos e maximos) dumafuncao. Como se vera no Capıtulo 3, a capacidade de determinar zeros sem variacao desinal permite representar componentes que sao indetectaveis pelos metodos numericosconvencionais.

2.1 Introducao

Os metodos numericos de bisseccao, falsa posicao, secante, Newton, Ridders e Brent[40] sao metodos que permitem determinar raızes ou zeros de funcoes reais. Todoseles se baseiam no Corolario do Teorema de Bolzano, o qual estabelece que existepelo menos um zero no intervalo [A, B] se f(A) e f(B) tiverem sinais contrarios, i.e.f(A).f(B) < 0. Nesta tese, qualquer zero que verifique a condicao anterior e denom-inado por zero secante ou zero de sinal variante. Como se ilustra na Figura 2.1(a),

15

16 CAPITULO 2. COMPUTACAO DE ZEROS E EXTREMOS DE FUNCOES

o zero e o resultado da interseccao entre a curva e o segmento [A,B], de modo queos valores de f imediatamente a esquerda e a direita do zero em [A,B] tem sinaiscontrarios.

No entanto, os metodos numericos convencionais falham na determinacao deum zero que e extremo, aqui designado por zero extremo ou zero de sinal invariante.Neste caso, a funcao a direita e a esquerda do zero em [A,B] mantem o mesmo sinal(Figura 2.1(b)). Isto ocorre com funcoes do tipo fN onde N e um numero par. Porexemplo, a funcao f(x, y) = (x + y)2 = 0 descreve uma recta, a esquerda e a direitada qual nao ha variacao de sinal de f .

Dizemos que existe um extremo no intervalo [A,B] se f(A) e f(B) tem sinaisiguais e existe variacao da monotonia de f , ou seja f(A).f(B) > 0 e f ′(A).f ′(B) < 0[40], como se ilustra na Figura 2.1(c). Se f(A) > 0 e f(B) > 0, subdivide-se ointervalo ate que se atinja o extremo, obtendo-se neste caso um mınimo. Depois, enecessario verificar se o extremo e um zero de f , tal como acontece na Figura 2.1(b).Analogamente, se f(A) < 0 e f(B) < 0, pode haver um maximo em [A,B], sendonecessario verificar se o mesmo e zero da funcao.

Apesar de existirem algoritmos de procura de extremos (maximos e mınimos)na literatura de analise numerica, nao existe um metodo numerico unico capaz dedeterminar zeros e extremos de uma funcao. O proposito deste capıtulo e precisamenteintroduzir um metodo numerico que determine zeros e extremos, o qual e designadopor metodo generalizado da falsa posicao (GFP).

2.2 Metodo generalizado da falsa posicao

O metodo generalizado da falsa posicao (GFP) e uma generalizacao do metodo dafalsa posicao. A sua novidade principal e a de poder calcular zeros e extremos atravesde uma unica formula de interpolacao.

A B

f(B)

f(A)

X A B

f(B)

-f(A)

X

f(A)

(b) (a)

A B

f(B)

-f(A)

X

f(A)

(c)

Figura 2.1: Ilustracao do metodo generalizado da falsa posicao.

Atente-se novamente na Figura 2.1 em que |f(A)| e |f(B)| sao os valores abso-

2.2. METODO GENERALIZADO DA FALSA POSICAO 17

lutos da funcao f em A e B, respectivamente. Os dois extremos do intervalo [A,B]sao considerados duas estimativas iniciais da raiz de f(x) = 0. A proxima estimativaX e a interseccao dos segmentos de recta (A, 0)(B, 0) e (A,−|f(A)|)(B, |f(B)|), quee dada por

X = A +|f(A)|

|f(A)|+ |f(B)| .(B − A) (2.1)

De facto, usando a equacao vectorial do segmento de recta definido por doispontos A e B em R, tem-se:

X(t) = A + t . (B − A) com t ∈ [0, 1] (2.2)

Repare-se que A = X(0), B = X(1), de modo que a proxima estimativa X

ocorre quando t = |f(A)||f(A)|+|f(B)| , o que resulta da seguinte igualdade de triangulos:

|f(A)|t

=|f(B)|1− t

(2.3)

Como se vera de seguida, a equacao (2.1) permite determinar quer zeros querextremos de uma funcao.

2.2.1 Determinacao dos zeros

A Figura 2.1(a) ilustra a determinacao dum zero. Se f(A).f(B) < 0, assumindo quef(A) < 0, e substituindo |f(A)| = −f(A) e |f(B)| = f(B) em (2.1), obtem-se

X = A +f(A)

f(B)− f(A).(B − A) (2.4)

que e a conhecida formula de interpolacao do metodo da falsa posicao [14], usada nadeterminacao de zeros secantes.

2.2.2 Determinacao dos extremos

A determinacao de um extremo (mınimo ou maximo) requer que se verifique a condicaof(A).f(B) > 0. Assumindo que f(A) > 0, e substituindo |f(A)| = f(A) e |f(B)| =f(B) em (2.1), obtem-se

X = A +f(A)

f(B) + f(A).(B − A) (2.5)


que permite calcular a proxima estimativa do extremo. Similarmente, se f(A) < 0,e substituindo |f(A)| = −f(A) e |f(B)| = −f(B) em (2.1), obtem-se a proximaestimativa X atraves de (2.5).

Recorde-se que a formula (2.1) funciona igualmente para mınimos e maximos,independentemente de f(A) e f(B) serem ambos positivos ou negativos. De facto,olhando para as equacoes (2.4) e (2.5), ve-se que e suficiente trocar o sinal de f(A)para determinar a proxima estimativa de um extremo atraves da tecnica da falsaposicao. Isto e ilustrado na Figura 2.1(c)) para um mınimo.

Mas a condicao f(A).f(B) > 0 so indica que e provavel que haja um extremoentre A e B. Para verificar se existe de facto um extremo em [A,B], tem de se calcular

os valores das diferencas finitas (ou derivadas discretas) d(x) = f(x+h)−f(x)h

em A e Be, depois, testar se tem sinais opostos, i.e. d(A).d(B) < 0. (Note-se que h e uminfinitesimo.) Deste modo, garante-se que existe um extremo entre A e B. Usa-se adiferenca finita de uma funcao num dado ponto, devido a duas razoes principais:

- Coerencia. O metodo da falsa posicao nao faz uso da derivada. Assim, poruma questao de coerencia metodologica, o metodo generalizado da falsa posicaotambem nao faz uso da derivada.

- Diferenciabilidade. Podem existir zeros e extremos onde a derivada nao existe.Por exemplo, a curva representada por f(x) = |x|1/3 na Figura 2.5(f) tem umapice em x = 0 que tambem e um mınimo. Nao se pode usar a derivada, poisneste ponto ela nao existe.

Se a diferenca finita d(X) da estimativa X e muito proxima de zero, entao X eum valor aproximado do extremo. No caso em que f(X) ≈ 0, com uma precisao de10−6 (erro numerico maximo admissıvel usado na implementacao), entao considera-seque X e um zero da funcao.

Portanto, a formula (2.1) e usada para calcular a proxima estimativa quer deum zero, quer de um extremo. Note-se que na literatura existem metodos numericosdistintos para calcular zeros e extremos. A maior novidade deste metodo consiste emutilizar a mesma formula para determinar zeros e extremos independentemente dascondicoes de diferenciabilidade da funcao.

2.2.3 Analise da convergencia

A convergencia do metodo GFP e garantida pelo seguinte teorema:

Teorema 2.1 Seja f uma funcao real contınua em [A,B] ⊂ R e f(A) . f(B) < 0 oud(A) . d(B) < 0. Se P ∈ [A, B] e tal que f(P ) = 0 ou d(P ) = 0, entao o metodoGFP gera uma sequencia Pn+∞

n=1 de pontos Pn que converge para P , dadas duasaproximacoes iniciais P1 = A e P2 = B. ¤

2.2. METODO GENERALIZADO DA FALSA POSICAO 19

Demonstracao. Sem perda de generalidade, isto e equivalente a provar a existencia deuma sequencia In+∞

n=1 de intervalos In = [Pn, Pn+1] convergindo para P , com P ∈ In

para todo o n.

Seja l(In) a amplitude do intervalo In de forma que l(I1) = |P2 − P1|, l(I2) =|P3 − P2|, e

l(I2)

l(I1)=|P3 − P2||P2 − P1|

Logo, de (2.1) obtem-se

l(I2)

l(I1)=|P2 + |f(P2)|

|f(P1)|+|f(P2)| . (P1 − P2)− P2||P2 − P1| =

|f(P2)||f(P1)|+ |f(P2)|

ou

l(I2) = l(I1) .|f(P2)|

|f(P1)|+ |f(P2)|

Em geral, tem-se

l(In) = l(In−1) .|f(Pn)|

|f(Pn−1)|+ |f(Pn)|

Tendo em conta que

0 <|f(Pn)|

|f(Pn−1)|+ |f(Pn)| < 1

temos l(In) < l(In−1) para todo o n. Assim, a sucessao l(In)+∞n=1 e monotona

decrescente, pelo que 0 ≤ l(In) ≤ |B−A|. Assim, pelo Teorema da Sucessao Monotonae Limitada, pode concluir-se que l(In)+∞

n=1 e convergente. Tendo em conta que P ∈ In

para todo n, pode tambem concluir-se que o metodo GFP converge para P , sendo Pum zero ou um extremo de f . ¥

A rapidez com que a sequencia Pn+∞n=1 converge para P e dada pela ordem de

convergencia α em

limn→+∞

|Pn+1 − P ||Pn − P |α = C

sendo C, denominada constante assimptota. Como se sabe, a ordem de convergenciado metodo da falsa posicao e 1 (veja-se [41]). Tendo em conta que o metodo GFPusa a mesma formula de iteracao, ainda que com uma troca de sinal, conclui-se que aordem de convergencia e tambem 1.

2.2.4 Implementacao do metodo GFP

A funcao em C++ que implementa o metodo GFP e a seguinte:


#include <math.h>

#define ACY 1.0E-6 precisao para a estimativa, funcao, e para a derivada discreta#define h 1.0E-7 infinitesimo da diferenca finita

Real GFP(Real (*func)(Real), Real A, Real B, int n) Metodo generalizado dafalsa posicao para determinar o zero ou o extremo da funcao func entre A e B. A solucaoe procurada ate a precisao ±ACY e nıvel de recursividade n.

if (n==0) return NO;

if (fabs(B-A)< ACY) return (A+B)/2;

Real fA=(*func)(A);

Real fB=(*func)(B);

Real X=FormulaIteracao(func,A,B,fA,fB);

Real fX=(*func)(X);

if (fabs(fX)< ACY) return X;

if (fA*fX<0) return GFP(*func,A,X,n-1);

if (fX*fB<0) return GFP(*func,X,B,n-1);

Real dX=((*func)(X+h)-(*func)(X))/h;

if (fabs(dX)< ACY) return X;

Real dA=((*func)(A+h)-(*func)(A))/h;

Real dB=((*func)(B+h)-(*func)(B))/h;

if (dA*dX<0) return GFP(*func,A,X,n-1);

if (dX*dB<0) return GFP(*func,X,B,n-1);

Real FormulaIteracao(Real (*func)(Real), Real A, Real B, Real fA, Real fB)

Permite determinar a nova estimativa, dada pela formula (2.1).

Real X=A+fabs(fA)/(fabs(fA)+fabs(fB))*(B-A);

return X;

2.3. RESOLUCAO DAS DESCONTINUIDADES 21

2.3 Resolucao das descontinuidades

A juntar a primeira pre-condicao dos algoritmos numericos —variacao de sinal—, haainda uma segunda pre-condicao que e a continuidade da funcao. Por exemplo, afuncao f(x, y) = y− 1

xtem uma descontinuidade em x = 0, pelo que ocorrera um erro

no calculo de f em (0, y), ∀y ∈ R. De modo a evitar este problema, o operador / (dedivisao) da classe Real foi implementado por forma a contornar este problema. Alias,como se pode constatar pelo codigo da classe Real, precaucoes semelhantes foramtomadas na implementacao do operadorˆe na implementacao das seguintes funcoes:sqrt, log, tg.

class Real

double value;

public:

...

Real Real::operator/(const Real &r1)

if (r1.zero())

DOM_VALIDITY = 1;

return MAXREAL;

else return Real(value / r1.value);

Real sqrt(const Real &r)

if (r < Real(0))

DOM_VALIDITY = 2;

return MAXREAL;

else return r.apply(sqrt);

Real log(const Real &r)

if (r <= Real(0))

DOM_VALIDITY = 3;

return MAXREAL;

else return r.apply(log);


Real tg(const Real &r)

// r != PI/2+K*Pi ; com K inteiro

// (r-PI/2)/Pi != K (No inteiro)

Real resto = fmod(r-PI/2, Real(PI));

if (resto.zero())

DOM_VALIDITY = 4;

return MAXREAL;

else return r.apply(tan);

Real Real::operator^ (const Real &Y)

if (zero() && (Y >= 0)) return Real(0);

if (*this > Real(0))

return exp(Y*log(*this));

else

int yy = Y.integer();

if (Y == yy) // expoente inteiro!

return (*this)^(yy);

DOM_VALIDITY = 5;

return MAXREAL;

O controlo da validade do domınio e feito atraves de DOM VALIDITY (variavelinteira), cujos valores possıveis sao os seguintes: 0 (domınio valido por defeito), 1(divisao por zero), 2 (raiz de numero negativo), 3 (logaritmo de numero negativo), 4(domınio da tangente) e 5 (potenciacao xy).

Por defeito, nao existem restricoes sobre o domınio, e neste caso DOM VALIDITY

= 0, o que significa que se assume que nao existe descontinuidades no intervalo emanalise. Em todas as situacoes, em que se efectua o calculo de f(x) (ou d(x)), eefectuado o teste ao valor de DOM VALIDITY, e so no caso em que DOM VALIDITY = 0

e que se deve usar o valor de f(x) (ou d(x)). Caso contrario, abandona-se o intervaloem analise porque existem pontos do intervalo onde a funcao ou a derivada discretanao esta definida.

2.4 Condicoes de paragem

O sub-programa que implementa as condicoes de paragem do algoritmo GFP e oseguinte:

2.4. CONDICOES DE PARAGEM 23

int CondicaoParagem(Real (*func)(Real), Real A, Real B, Real *X,

Real fX, Real dX, int n)

if (n==0) return STOP1

if (fabs(fX)< ACY) return 1; // termino

if (fabs(dX)< ACY) return 1; // termino

if (fabs(B-A)< ACY)

// O intervalo e pequeno, logo pode haver uma assimptota vertical

if (|fX| > Maior(|fA|, |fB|))

return STOP2; // assimptota vertical

else

*X = (A+B)/2;

return 1; // termino

return 0; // o algoritmo deve continuar

Como se pode constatar, as condicoes de paragem sao as seguintes:

(1) O nıvel maximo de recursividade foi alcancado. Este valor e alcancado pordecrementacao e teste da condicao (n==0).

(2) Um zero da funcao foi encontrado. Isto acontece quando a condicao (fabs(fX)<ACY) e verdadeira.

(3) Um zero da derivada foi encontrado. Isto acontece quando a condicao (fabs(dX)<ACY) e verdadeira.

(4) Uma assimptota vertical foi encontrada. A deteccao da assimptota vertical efeita atraves do teste (|fX| > Maior(|fA|, |fB|)).

A condicao |fX| > Maior(|fA|, |fB|) decorre da definicao matematica deassimptota vertical num dado ponto x = a,

x = a e assimptota vertical se limx→a

|f(X)| = +∞

ou seja, a medida que o valor de x tende para a o valor de |f(x)| tende para +∞(Figura 2.2). Assim, e de esperar que o valor de |f(x)| seja maior que o maior dosvalores |f(A)| e |f(B)|. Repare-se que esta condicao so e testada no caso em que ointervalo e relativamente pequeno |B-A| < ACY. Caso contrario, nao ha garantia deexistir assimptota.


B A

X

f(B)

f(A)

f(X)

a B A X

f(B) f(A)

f(X)

a

Figura 2.2: Assimptota vertical

2.5 Metodo TGFP

Como se pode ver pela Figura 2.3(a), o metodo da falsa posicao pode convergirmuito lentamente para um zero, pois sao necessarias muitas iteracoes para chegara solucao. Por exemplo, o zero da funcao f(x) = tg(x)tg(x)− 103 no intervalo [1.3, 1.4](Figura 2.5(d)) foi encontrado apos 221 aproximacoes como se mostra na Tabela 2.4.Pode ate acontecer que um zero nao seja encontrado, mesmo que o algoritmo convirja.Por exemplo, conforme se mostra na Tabela 2.5, o algoritmo converge para o zero dafuncao f(x) = x ex−10 (Figura 2.5(e)) no intervalo [−10.0, 10.0], mas a convergencia etao lenta que o zero so pode ser encontrado para alem do numero maximo de iteracoespre-definido (MAXIT=500).

De modo a superar esta dificuldade, e seguindo o espırito do metodo da biseccaoque se encontra nos algoritmos de Ridders [44, 40] e de Brent [11, 40], considera-seque a estimativa X e dada pela media

X =P + Q

2(2.6)

de duas estimativas previamente calculadas P e Q, em que P e a estimativa calculadapelo metodo GFP e Q ∈]A,B[ e a abcissa do ponto que resulta da interseccao das duasrectas tangentes nos pontos (A, f(A)) e (B, f(B)), como se mostra na Figura 2.3(b) e(c). Assim sendo, a abcissa Q e a solucao x do seguinte sistema de equacoes:

y − f(A) = d(A).(x− A)

y − f(B) = d(B).(x−B)

ou seja,

2.5. METODO TGFP 25

X = Q =f(B)− f(A) + A.d(A)−B.d(B)

d(A)− d(B)(2.7)

onde d(A) e d(B) sao as diferencas finitas em A e B, respectivamente. Mas, comose mostra na Figura 2.3(d), pode acontecer que Q 6∈]A,B[. Neste caso, e de modo agarantir a convergencia do metodo, a proxima estimativa X e dada por P .

(b)(a)

A

B

f(B)

f(A)

P A

B

f(B)

f(A)

P QX

A B

f(B)

f(A)

P

Q

X

A

B

f(B)

f(A)

X=P

Q

(d)(c)

Figura 2.3: Metodo TGFP.

Analisando as Tabelas 2.4 e 2.5, observa-se que o metodo TGFP e muito maisrapido que o metodo GFP. Em geral, o metodo TGFP e mais rapido que o metodoGFP porque a estimativa X e interpolada por duas estimativas preliminares P e Q.Deste modo, a diferenca entre o algoritmo GFP e o algoritmo TGFP esta na formulade iteracao. Assim, e em conformidade com as equacoes (2.6) e (2.7), a formula deiteracao do algoritmo TGFP pode ser implementada da seguinte forma:

Real FormulaIteracao(Real (*f)(Real), Real A, Real B, Real fA, Real fB)

Real P=A+fabs(fA)/(fabs(fA)+fabs(fB))*(B-A);

if (fabs(d(A)-d(B))>DACC)

Real Q=(fB-fA+A.d(A)-B.d(B))/(d(A)-d(B));

if ((Q>A)&&(Q<B))


X=(P+Q)/2.0;

else X = P;

return X;

Logo o algoritmo TGFP pode ser implementado da seguinte maneira:

Real TGFP(Real (*func)(Real), Real A, Real B, Real fA, Real fB,

Real dA, Real dB, int n)

Real X=FormulaIteracao(func,A,B,fA,fB);

Real fX=(*func)(X);

if (DOM_VALIDITY > 0) return ERRO;

Real dX=((*func)(X+h)-fX)/h;

if (DOM_VALIDITY > 0) return ERRO;

int Tipo = CondicaoParagem(func,A,B,&X,fX,dX,n);

if ((Tipo == STOP1) || (Tipo == STOP2))

// Nıvel maximo de recursividade ou assimptota vertical

return NADA;

if (Tipo == 1) // Zero em X

return X;

if (fA*fX<0) return TGFP(*func,A,X,fA,fX,dA,dX,n-1);

if (fX*fB<0) return TGFP(*func,X,B,fX,fB,dX,dB,n-1);

if (dA*dX<0) return TGFP(*func,A,X,fA,fX,dA,dX,n-1);

if (dX*dB<0) return TGFP(*func,X,B,fX,fB,dX,dB,n-1);

2.5.1 Convergencia do metodo TGFP

A convergencia do metodo TGFP e garantida pelo seguinte teorema:

Teorema 2.2 Seja f uma funcao real contınua em [A,B] ⊂ R e f(A) . f(B) < 0 oud(A) . d(B) < 0. Se X ∈ [A,B] e tal que f(X) = 0 ou d(X) = 0, entao o metodoTGFP gera uma sequencia Xn+∞

n=1 de pontos Xn que converge para X, dadas duasestimativas iniciais X1 = A e X2 = B. ¤

Demonstracao. A proxima estimativa Xi calculada pelo metodo TGFP e dada porXi = Pi+Qi

2em que Pi, Qi sao pontos de duas sequencias Pn+∞

n=1 e Qn+∞n=1 respec-

tivamente. O calculo de Pi e dado pelo metodo GFP, que gera a sequencia Pi∞i=1

convergente para X (provado pelo Teorema 2.1).

2.6. METODO NGFP 27

E entao suficiente provar que a sequencia Qn∞n=1 gerada pelo calculo de Qtambem converge para X, porque a soma de duas sucessoes convergentes e umasucessao convergente. Sem perda de generalidade, pode assumir-se que d(Xn−2) > 0e d(Xn−1) < 0 de forma que |d(Xn−2)| = d(Xn−2) e |d(Xn−1)| = −d(Xn−1). Entao aformula (2.7) pode ser escrita do seguinte modo:

Qn =f(Xn−1)− f(Xn−2) + Xn−2.|d(Xn−2)|+ Xn−1.|d(Xn−1)|

|d(Xn−2)|+ |d(Xn−1)|ou, de outro modo,

Qn = Xn−1 + t.(Xn−2 −Xn−1) (2.8)

que e a equacao vectorial da recta definida pelos dois pontos Xn−1 e Xn−2, sendo

t =f(Xn−1)− f(Xn−2)

(|d(Xn−1)|+ |d(Xn−2)|).(Xn−2 −Xn−1)+

|d(Xn−2)||d(Xn−1)|+ |d(Xn−2)|

Seja l(In) a amplitude do intervalo In de forma que l(In) = |Qn−Xn−1|, l(In−1) =|Xn−1 −Xn−2|. Entao, de (2.8) obtem-se

l(In)

l(In−1)=

|Qn −Xn−1||Xn−1 −Xn−2| =

|Xn−1 + t.(Xn−2 −Xn−1)−Xn−1||Xn−1 −Xn−2| = t (2.9)

Mas, como foi descrito anteriormente, a nova estimativa Qn so e tida em conta noprocesso de convergencia se esta pertence ao intervalo ]Xn−2, Xn−1[; caso contrario, Qn

e igual a estimativa Pn ∈]Xn−2, Xn−1[ determinada pelo metodo GFP, ou seja Qn = Pn.De (2.8) e (2.9) tem-se que l(In) < l(In−1), pois 0 < t < 1 para todo n. Assim, asequencia l(In)+∞

n=1 e monotona decrescente e limitada, visto que 0 < l(In) ≤ |B−A|.Pode, pois, concluir-se pelo Teorema da Sucessao Monotona e Limitada que l(In)+∞

n=1

e convergente. Tendo em conta que X ∈ In para todo o n, conclui-se que o metodoTGFP e convergente porque a sequencia Qn∞n=1 converge para X, sendo X um zeroou um extremo de f . Isto completa a demonstracao. ¥

2.6 Metodo NGFP

Uma outra formula de iteracao consiste em determinar os pontos N1 e N2 de inter-seccao de cada uma das rectas tangentes nos extremos do intervalo com o segmento(A, 0)(B, 0), sendo de seguida calculada uma media X entre a media N desses doisvalores e o valor P dado pelo GFP:


X =P + N1+N2

2

2(2.10)

A B

f(B)

f(A)

P

N1

X

N2

(N1+N2)/2

Figura 2.4: Metodo NGFP.

No entanto, as contribuicoes N1 e N2, so sao aceites quando pertencem aointervalo ]A,B[, pois so deste modo se pode garantir a convergencia. No pior doscasos (quando N1 e N2 nao pertencem a ]A, B[), o metodo NGFP tem o mesmocomportamento que o metodo GFP.

Este metodo mostrou-se ainda mais rapido que os metodos GFP e TGFP, comose pode verificar na Tabela 2.7. A diferenca essencial entre o algoritmo NGFP e osalgoritmos GFP e TGFP reside na formula de iteracao. A formula de iteracao 2.10do algoritmo NGFP pode ser implementada da seguinte forma:

Real FormulaIteracao(Real (*func)(Real), Real A, Real B,

Real fA, Real fB)

Real P=A+fabs(fA)/(fabs(fA)+fabs(fB))*(B-A);

int m = 0;

Real N1 = 0;

Real N2 = 0;

N1 = A - fA/d(A);

if ((N1 > A) && (N1 < B)) // Se N1 pertence a ]A,B[

m++;

else N1 = 0;

N2 = B - fB/d(B);

if ((N2 > A) && (N2 < B)) // Se N2 pertence a ]A,B[

m++;

else N2 = 0;

if (m != 0)

Real N = (N1 + N2)/Real(m);

X = (P + N)/Real(2);

2.7. RESULTADOS EXPERIMENTAIS 29

else X = P; // Ponderac~ao dada pelo metodo GFP

return X;

2.7 Resultados experimentais

2.7.1 Resultados experimentais dos metodos GFP, TGFP eNGFP

Note-se que, na literatura existem metodos distintos para determinar zeros e extremosde uma funcao. Os metodos GFP, TGFP e NGFP sao metodos capazes de determinarzeros e extremos de uma funcao real, atraves de uma unica formula de iteracao.

Da discussao tida anteriormente sobre os metodos GFP, TGFP e NGFP, eexpectavel que o metodo NGFP convirja mais rapidamente para a solucao do queo metodo TGFP e que este, por sua vez, convirja mais rapidamente que o metodoGFP. Contudo, uma analise breve das Tabelas 2.1, 2.2, 2.3, 2.4, 2.5 e 2.6 referentes ascurvas da Figura 2.5 mostra que nem sempre assim e.

A Tabela 2.1 mostra a sequencia de estimativas que convergem para o mınimoda funcao (Figura 2.5(a)) em x = 0. Neste caso o metodo TGFP tem pior desempenhodo que o metodo GFP. Isto deve-se, em parte, ao facto de a primeira estimativa dadapelo TGFP estar mais afastada da solucao do que a dada pelo metodo GFP. Alemdisso, a estimativa dada pela interseccao das rectas tangentes esta mais afastada dasolucao que a dada pelo GFP.

Na Tabela 2.2 mostram-se as estimativas a convergir para um maximo em ]1, 1.8[da funcao representada na Figura( 2.5(b)). Repare-se que, neste caso, o metodo NGFPe precisamente o metodo GFP, uma vez que os valores N1 e N2 que sao calculados emcada aproximacao nunca se encontram dentro do intervalo, e logo nao sao usados naproxima estimativa.

As Tabelas 2.3 e 2.4 apresentam resultados que convergem para um zero secanteem x ≈ −1.44732586 e x ≈ 1.35471068, respectivamente. Neste caso, o desempenhodos algoritmos esta de acordo com as expectativas, sendo notorio o melhor desempenhodo algoritmo NGFP que encontrou a solucao apos 7 estimativas. Por contraposicao, ometodo GFP so encontrou a solucao apos 221 estimativas, quando o maximo admissıvelde estimativas e de 500 (MAXIT). Estes resultados poem em evidencia a conhecidadeficiencia do metodo da falsa posicao (FP). Esta deficiencia tem que ver com o factodo algoritmo FP ter convergencia muito lenta quando o zero da funcao esta muitoafastado do extremo do intervalo com menor valor absoluto de f .

Este problema de convergencia do metodo FP torna-se ainda mais problematico


Tabela 2.1: f(x) = x2 + 5 em]− 1, 8[, mınimo local

GFP TGFP NGFP-0.28000000 1.61000000 1.281875040.28763040 0.22856995 0.082894070.00369422 -0.35937681 -0.40969925-0.13704965 -0.06429310 -0.15950574-0.06654588 0.08178817 -0.03808148-0.03141033 0.00873823 0.02237352-0.01385634 -0.02777000 -0.00785110-0.00508090 -0.00951522 0.00726054-0.00069333 -0.00038845 -0.000295270.00150044 0.00417489 0.003482610.00040355 0.00189324 0.00159366-0.00014489 0.00075241 0.000649190.00012933 0.00018200 0.00017696-0.00000778 -0.00010320 -0.000059150.00006077 0.00003942 0.000058900.00002649 -0.00003186 0.000000120.00000935 0.000003800.00000078 -0.00001400-0.00000349 -0.00000507-0.00000135 -0.00000061-0.00000028 0.00000162

0.00000053-0.00000001

Tabela 2.2: f(x) = sin2(x)+1 em]1, 1.8[, maximo local

GFP TGFP NGFP1.37371159 1.40667912 1.373711591.58757965 1.60311553 1.587579651.47961821 1.50500763 1.479618211.53349062 1.55406835 1.533490621.56052763 1.57859119 1.560527631.57405423 1.56632991 1.574054231.56729077 1.57246056 1.567290771.57067250 1.56939526 1.570672501.57236337 1.57092794 1.572363371.57151794 1.57016162 1.571517941.57109522 1.57054481 1.571095221.57088386 1.57073640 1.570883861.57077818 1.57083219 1.570778181.57083102 1.57078432 1.570831021.57080460 1.57080828 1.570804601.57079139 1.57079632 1.570791391.57079799 1.570797991.57079469 1.570794691.57079634 1.57079634


Tabela 2.3: f(x) = x3

6− x2

2+x+3

em ]− 3.0, 2.0[, zero secante

GFP TGFP NGFP0.37500000 -0.75694444 -0.88112738-0.53319057 -1.58469673 -1.52097775-1.03518015 -1.30362522 -1.43915728-1.27380512 -1.44363272 -1.44718798-1.37669220 -1.48098871 -1.44732408-1.41899245 -1.45482024 -1.44732581-1.43602899 -1.44827342-1.44283263 -1.44663918-1.44554045 -1.44739097-1.44661668 -1.44717043-1.44704420 -1.44730324-1.44721400 -1.44733645-1.44728142 -1.44732281-1.44730820 -1.44732771-1.44731883 -1.44732552-1.44732305 -1.44732620-1.44732473 -1.44732586-1.44732540-1.44732566

Tabela 2.4: f(x) = tg(x)tg(x)−103

em ]1.3, 1.4[, zero secante

GFP TGFP NGFP1.30338855 1.34653282 1.346652741.30662101 1.36879627 1.360433821.30970158 1.35599359 1.354791161.31263431 1.35307793 1.354699231.31542317 1.35461536 1.354710481.31807200 1.35501151 1.354710441.32058464 1.35476199 1.354710441.32296484 1.354699571.32521639 1.354720641.32734307 1.354712981.32934869 1.354711071.33123710 1.354710681.333012231.334678071.336238671.33769814. . . apos 221iteracoes1.354710440

Tabela 2.5: f(x) = x ex − 10 em]− 10.0, 10.0[ zero secante

GFP TGFP NGFP-9.99909196 -0.45409143 -0.45407077-9.99818396 4.31865924 4.31868396-9.99727600 1.60329676 1.43424159-9.99636809 2.60588042 1.98587097-9.99546021 1.94194915 1.74995900-9.99455238 1.75574808 1.74502469-9.99364459 1.71324321 1.74552986-9.99273683 1.73999411 1.74552803-9.99182912 1.74669341-9.99092146 1.74443607-9.99001383 1.74554623-9.98910624 1.74525966-9.98911000 1.74546550-9.98819869 1.74551696. . . 1.74552982. . . 1.74552572. . . 1.74552791. . . 1.74552846>MAXIT 1.74552819

Tabela 2.6: f(x) = |x|1/3 em]− 0.5, 1.5[, zero extremo

GFP TGFP NGFP0.31891712 0.01012557 0.65945851-0.05991142 -0.09913853 0.053002310.07804268 -0.02467953 -0.115323730.00602813 -0.00010021 -0.06316061-0.01490493 0.00170721 -0.00338224-0.00287112 0.00040563 0.009740400.00103121 0.00000471 0.00439930-0.00059016 -0.00002311 0.000338150.00014533 -0.00000559 -0.00075880-0.00013805 -0.00000001 -0.000406640.00000242 -0.00000071 -0.00002280-0.00002656 -0.00000036 0.00006355-0.00000657 0.00002928-0.00000133 0.000002150.00000036 -0.00000503-0.00000060 -0.00000267-0.00000030 -0.00000017

0.000000370.00000010


-2 2 4 6

10

20

30

40

Untitled-8 1

-1 1 2 3 4 5 6

1.2

1.4

1.6

1.8

2

Untitled-8 1

(a) f(x) = x2 + 5 (b) f(x) = sin2(x) + 1

-3 -2 -1 1 2

-8

-6

-4

-2

2

4

Untitled-8 1

1.25 1.3 1.35 1.4

5000

10000

15000

20000

25000

Untitled-8 1

(c) f(x) = x3

6− x2

2+ x + 3 (d) f(x) = tg(x)tg(x) − 103

-10 -7.5 -5 -2.5 2.5 5

20

40

60

numerica.nb 1

-10 -5 5 10

0.5

1

1.5

2

(e) f(x) = x ex − 10 (f) f(x) = |x|1/3

Figura 2.5: Graficos de funcoes reais com zeros e extremos.

quando a solucao nao e atingida dentro do numero maximo de estimativas (MAXIT),como acontece com a funcao da Tabela 2.5. Recorde-se que o algoritmo TGFPultrapassa esta dificuldade porque calcula a proxima estimativa X atraves da mediaaritmetica da estimativa P dada pelo algoritmo GFP e a estimativa Q dada pelaabcissa do ponto de interseccao das rectas tangentes em (A, f(A)) e (B, f(B)). Istoexplica porque e que o algoritmo TGFP e em geral mais rapido do que o algoritmoGFP. Mas, em alguns casos, o algoritmo TGFP e um pouco mais lento que o algoritmoGFP, como ja se viu. Isto acontece porque a proxima estimativa P esta mais proximado zero do que a estimativa X.

A Tabela 2.6 refere-se ao calculo de um extremo que e um zero. Os resultadosapresentados nesta tabela sao particularmente importantes porque mostram que os tresmetodos numericos conseguem determinar extremos onde a funcao nao e diferenciavel.


Isto significa que os metodos funcionam independentemente de a funcao ser ou naodiferenciavel. Neste caso concreto o metodo NGFP torna-se mais lento, devido ao factode que uma das estimativas N1 ou N2 nao se encontrar dentro do intervalo ]− 0.5, 1.5[em analise. Digno de nota e tambem o facto de o algoritmo NGFP ser particularmenterapido na determinacao dos zeros secantes. Isto deve-se ao facto de a probabilidadede N1 e N2 pertencerem ao intervalo onde se procura o zero ou o extremo da funcaoser muito elevada.

2.7.2 Comparacao com outros metodos

No estudo comparativo entre os metodos aqui propostos (GFP, TGFP e NGFP) e osmetodos numericos convencionais (bisseccao, falsa posicao, Newton, Ridders e Brent),so se teve em conta funcoes contınuas e os zeros secantes (com variacao de sinal), poisos zeros extremos so podem ser calculados pelos metodos GFP, TGFP e NGFP. Note-seque, nestas condicoes, o metodo GFP e equivalente ao metodo classico da falsa posicao.Por esta razao, o metodo FP nao foi incluıdo na comparacao. Esta comparacaorequereu a implementacao dos metodos classicos descritos no livro ”Numerical Recipesin C”de Press et al. [40].

Funcao Intervalo GFP TGFP NGFP Bisseccao Newton Ridders BrentNf Nf Nf Nf Nf Nf Nf

ex − 1 [−5, 10] > 500 13 9 23 147 − 25ln(x) [0.1, 5] 52 14 6 21 5 − 24x3 [−1, 5] > 500 4 4 8 12 − 24xex − 10 [−5, 6] > 500 18 9 25 Divergiu − 25tg(x)tg(x)−103

[1.3, 1.5] > 500 17 9 19 Divergiu − 19

x2 − 1 [0, 5] 59 15 7 21 7 6 24x3−2x−5 [2, 5] 64 14 7 23 3 5 23cos(x) [0.1, 3] 25 7 5 21 5 3 23sin(x) [−0.2, 0.5] 21 11 5 17 3 4 21Nf −−− −−− 12.5 6.7 19.7 −−− −−− 23.1α 1 > 1 > 1 1 2 > 1 > 1

Tabela 2.7: Comparacao entre metodos numericos em termos do numero Nf deestimativas.

Pela analise da Tabela 2.7, verifica-se que os metodos aqui propostos sao os queapresentam melhores resultados, uma vez que necessitam de um menor numero Nf deaproximacoes ou estimativas para encontrar um zero duma funcao. Alem disso, paraas funcoes indicadas na Tabela 2.7, o numero medio Nf de estimativas usadas pelometodo NGFP e cerca de metade do numero medio de estimativas dadas pelo metodoTGFP. No que diz respeito ao metodo GFP, constatou-se que em quatro das funcoeso processo de convergencia necessitava de um maior numero de estimativas do que oadmitido por MAXIT, pelo que deixou de fazer sentido o calculo de Nf .


Relativamente aos metodos convencionais pode dizer-se que o metodo da bis-seccao foi aquele que apresentou melhores resultados. O metodo de Newton divergiuem duas das funcoes porque alguma das seguintes pre-condicoes usuais nao estavamsatisfeitas, nomeadamente: (i) continuidade de f ′ e f ′′ no intervalo em analise; (ii)anulamento de f ′ e f ′′; (iii) proximidade da estimativa inicial relativamente a solucao.Note-se que o metodo de Ridders chegou mesmo a quebrar em cinco das funcoeslistadas.

Atraves dos resultados apresentados, pode-se concluir que, nos casos onde existevariacao de sinal da funcao, os metodos TGFP e NGFP apresentam uma media deresultados melhor do que os restantes metodos. E interessante verificar que o metodoNGFP e ate mais rapido que o metodo de Newton que tem convergencia quadratica,o que sugere que a convergencia do metodo NGFP podera ser superior a 2.

2.8 Sumario

Os algoritmos GFP, TGFP e NGFP foram desenhados e implementados para resolverum problema de computacao grafica, o qual consiste em representar componentes deuma curva sem que estas tenham variacao de sinal. Por exemplo, a curva descrita pelafuncao f(x, y) = (x2 + y2− 1)2 em R2 nao pode ser desenhada pelos algoritmos que seapoiam no Corolario de Bolzano porque a funcao e sempre positiva ou nula. Mas, osmetodos foram concebidos duma forma mais geral, pois permitem determinar nao sozeros (com ou sem variacao de sinal), mas tambem extremos (mınimos e maximos).

As vantagens principais dos algoritmos GFP, TGFP e NGFP relativamente aosalgoritmos numericos classicos sao as seguintes:

- Cada um deles permite determinar zeros e extremos de uma funcao real atravesde uma unica formula de interpolacao.

- Nenhum deles requer as pre-condicoes usuais dos algoritmos classicos, nomeada-mente: continuidade, variacao de sinal (Corolario de Bolzano) e diferenciabili-dade.

- Cada um deles pode ser adaptado para determinar os pontos de inflexao, emboraeste assunto nao tenha sido objecto de discussao nesta tese. Para isso, serianecessario calcular a segunda diferenca dividida.

Como se vera nos proximos capıtulos, as caracterısticas destes novos metodosserao fundamentais para poder representar curvas e superfıcies implıcitas na suatotalidade. De facto, a grande maioria dos algoritmos actuais apresentam um deficeassinalavel na representacao destes objectos geometricos, pois nao preservam a suaforma topologica: componentes, furos, pontos de auto-interseccao, dimensao, etc.

Capıtulo 3

Curvas Implıcitas Nao-Homogeneas

Este capıtulo descreve um novo algoritmo de amostragem e poligonizacao de curvasdefinidas implicitamente por funcoes reais. A amostragem de pontos duma curvaimplıcita e levada a cabo pela articulacao de:

• Subdivisao ou particao do domınio da funcao que descreve a curva.

• Calculo dos pontos de interseccao entre a curva e cada segmento de recta queparticiona o domınio da funcao, o que pode ser feito atraves de qualquer dos tresmetodos numericos introduzidos no capıtulo anterior.

Estes metodos numericos nao fazem uso do calculo de derivadas, o que permitemaior robustez na amostragem de curvas com singularidades (i.e. apices e auto-interseccoes). Alem disso, o metodo pode e e ainda usado na deteccao e calculo depontos isolados e de varias componentes que a curva pode eventualmente ter. Comose vera, a capacidade do metodo numerico em determinar zeros e extremos dumafuncao, atraves duma unica formula de interpolacao, permite determinar componentesda curva com e sem variacao de sinal.

O calculo dos pontos isolados e das componentes sem variacao de sinal sao,porventura, as principais inovacoes do algoritmo descrito neste capıtulo.

3.1 Subdivisao do domınio

Ha varias formas de particionar o domınio duma funcao associada a uma curvaimplıcita em R2. A subdivisao quadtree tem sido bastante usada na literatura. Nestecapıtulo, usa-se um metodo diferente de subdivisao que e o metodo BSP (Binary SpacePartitioning). O metodo BSP baseia-se na particao binaria de um subespaco de R2,sendo esta particao um caso particular de subdivisao binaria multi-dimensional, emque linhas subdividem R2, planos subdividem R3 e hiperplanos subdividem Rn (n > 3).

35

36 CAPITULO 3. CURVAS IMPLICITAS NAO-HOMOGENEAS

Uma condicao necessaria, mas nao suficiente, para subdividir um subespaco (parte dodomınio) e que a curva ou parte dela se encontre nesse subespaco. Como se vera maisadiante, esta subdivisao depende essencialmente da variacao da curvatura da curva.Quanto maior for a variacao da curvatura, mais subdivisoes havera no subespaco quecontem a curva, i.e. a particao binaria e adaptativa relativamente a forma da curva.

3.1.1 Estrutura de dados: arvore BSP

As arvores BSP foram introduzidas por Fuchs, Kedem e Naylor [21] em 1980, sendousadas na determinacao das faces visıveis de um objecto geometrico sem usar o metodode z-buffer. Desde entao, as arvores BSP tem sido usadas em varias sub-areas dacomputacao grafica [52, 53, 54].

Em termos geometricos, o processo de criar uma arvore BSP consiste em dividirum d-espaco em dois d-subespacos atraves de um hiperplano de dimensao d − 1. Oprocesso e repetido ate que o espaco inicial esteja dividido numa hierarquia de regioesou subespacos convexos, i.e. uma arvore BSP e gerada por subdivisao binaria de umespaco d-dimensional em subespacos convexos, em que cada subespaco fica associadoa um no da arvore BSP. No espaco de duas dimensoes, o hiperplano e uma linha recta.Assim, se se dividir uma regiao convexa 2-dimensional por uma linha recta, obtem-seduas sub-regioes convexas, sendo esta a unica operacao que e necessaria para construiruma arvore BSP.

Atente-se a Figura 3.1. Seja Ω uma regiao aberta em R2 e l uma recta quesubdivide Ω atraves de algum criterio. Entao Ω e dividido por l em duas novas regioesabertas, Ω+ = Ω ∩ h+ e Ω− = Ω ∩ h−, onde h+ e h− sao os semi-espacos a direitae a esquerda de l em R2, respectivamente. Note-se que uma regiao 1-dimensionalli = Ω∩ l, (i = 1, 2, ...), designada por segmento de bisseccao de Ω, tambem e gerada.De um modo semelhante, qualquer uma destas novas regioes pode ser dividida ate quealgum criterio de paragem seja satisfeito. Terminado este processo, a reuniao de todosos subespacos associados aos nos terminais da arvore BSP sera o espaco inicial Ω. Acorrespondente estrutura de suporte da arvore BSP em C++ e a seguinte:

class BSPtree

BSPnode *root;

class BSPnode

List<BSPsubline> *Fr;

List<Point> *F;

List<Point> *I;

BSPsubline *l;

BSPnode *left, *right, *next;

3.1. SUBDIVISAO DO DOMINIO 37

class BSPsubline

Point *start;

Point *end;

class Point

double *x;

double *y;

double *z;

BSPnode *left, *right;

A classe BSPnode caracteriza cada no da arvore BSPtree. A variavel Fr denotaa fronteira Fr(Ω) do subespaco convexo Ω, a qual consiste numa lista de segmentosde recta que delimitam o espaco Ω; F e a lista de pontos da curva resultantes dainterseccao Fr(Ω) ∩ C; I e a lista de pontos isolados da curva em cada subespaco;l e a linha responsavel pela subdivisao do espaco Ω em dois novos subespacos, Ω−e Ω+; left denota Ω−, enquanto que right denota Ω+; next e usado aquando davisualizacao da curva C, sendo um apontador para o proximo no terminal mais aesquerda da arvore BSP. A visualizacao da arvore BSP atraves do apontador next

permite reduzir a sua complexidade a travessia duma lista simplesmente ligada, cujosnos sao os nos terminais da arvore BSP.

A classe BSPsubline serve para definir cada segmento de recta l de particao doespaco Ω, assim como cada segmento da fronteira Fr de Ω.

A classe Point caracteriza cada ponto amostrado de uma curva implıcita, mas,como se vera no proximo capıtulo, tambem serve para caracterizar cada ponto amostradode uma superfıcie implıcita.

3.1.2 Linha de bisseccao

Dado um subespaco Ω, a determinacao da sua linha de bisseccao BSP e feita doseguinte modo:

• Determinam-se os dois pontos mais afastados P , Q de Fr(Ω) ∩ C.

• Determina-se a mediatriz l do segmento PQ.

A linha de bisseccao e precisamente a mediatriz de PQ (Figura 3.2). O facto dese escolher os pontos mais afastados P e Q permite que seja efectuada uma subdivisaomais equitativa do subespaco que se esta a processar. Este processo permite que aarvore BSP seja mais equilibrada, ficando deste modo minimizados os problemas destack overflow em espaco de memoria.


l1

+ -




l1

l3 l2

l1

l3 l2

l6 l5 l7 l4

(d) Terceira partição

l1 +

-

l2 l1

l3

l2 l1

l3

l6

l7

l5

l4

Figura 3.1: Subdivisao e arvore BSP.

3.1. SUBDIVISAO DO DOMINIO 39

P

Q

l

P

Q

P

Q

l

(a) (b) (c)

l

l1

l1

l1

Figura 3.2: Determinacao da linha BSP para tres subespacos.

No caso particular da linha de bisseccao do espaco inicial Ω, o processo dedeterminacao e um tudo nada diferente, como se descreve de seguida:

• Determina-se os pontos de F = Fr(Ω) ∩ C.

• Se #F ≥ 2, determina-se os pontos mais afastados P , Q de F ; a linha debisseccao e a mediatriz l do segmento PQ.

• Caso contrario, determina-se um ponto X de um dos segmentos AB de Fr(Ω)que seja um mınimo absoluto nao nulo de f , i.e. |f(X)| < |f(Y )| para todo oY ∈ AB, com f(X) 6= 0. Neste caso, a linha de bisseccao e precisamente a linhaperpendicular a AB em X.

A existencia de um mınimo absoluto X no segmento AB de Fr(Ω) indica quea curva esta mais proxima de X do que qualquer outro ponto de AB. Deste modo,a linha de bisseccao l perpendicular a AB tem grande probabilidade de intersectar acurva, a nao ser que o atractor de X seja um ponto isolado. Este processo permite,pois, amostrar curvas com componentes muito pequenas, mesmo que o espaco inicialseja grande.

3.1.3 Determinacao dos subespacos Ω− e Ω+

Os subespacos Ω− e Ω+ que resultam da bisseccao do espaco Ω atraves da linha l saodeterminados pelo seguinte algoritmo:

Algoritmo 3.1 Subspaces

INPUT:

(a) Ω: subespaco de R2

(b) l: linha de subdivisao de Ω

OUTPUT:


(a) Ω−: Ω− ⊂ Ω

(b) Ω+: Ω+ ⊂ Ω

Begin

1. Ω− ← ∅

2. Ω+ ← ∅

3. Determinar os pontos A e B que resultam da interseccao Fr(Ω) ∩ l, i.e. a linhali que subdivide Fr(Ω).

4. Criar o segmento de bisseccao li, atraves dos dois pontos anteriores.

5. Redistribuir os segmentos da Fr(Ω) em Fr(Ω−) e Fr(Ω+). Isto requer a classi-ficacao de cada um dos segmentos em relacao a li.

6. Adicionar li a fronteira de Ω− e Ω+, completando deste modo as fronteiras dossubespacos.

7. Determinar os pontos da curva C ∩ li em Ω.

8. Redistribuir os pontos de Fr(Ω) ∩ C por Fr(Ω−) ∩ C e Fr(Ω+) ∩ C. Isto requer aclassificacao de cada um dos pontos em relacao a li.

9. Adicionar os pontos da curva C ∩ li a Fr(Ω−) ∩ C e Fr(Ω+) ∩ C.

End

Basicamente, este algoritmo cria dois novos nos BSP, left e right, associadosa Ω− e a Ω+, respectivamente (passos 1 e 2 do algoritmo Subspaces). Estes nosficam hierarquicamente ligados ao no de Ω. Os restantes passos do algoritmo saoexclusivamente dedicados a actualizar as fronteiras Fr(Ω+) e Fr(Ω−) dos dois novossubespacos, como se ilustra na Figura 3.3(a) e (b). Apos determinar os pontos A e Bque resultam da interseccao Fr(Ω)∩l (passo 3), cria-se o segmento de bisseccao li = AB(passo 4). Na Figura 3.3(c), o ponto A subdivide s5 em dois novos segmentos, s51 e s52,ao passo que o ponto B subdivide s3 em dois novos segmentos, s31 e s32. Os segmentosde Fr(Ω) = s1, s2, s3, s4, s5 sao, entao, redistribuıdos pelas fronteiras dos subespacosΩ+ e Ω− (passo 5), i.e. Fr(Ω+) = li, s31, s4, s51 e Fr(Ω−) = li, s52, s1, s2, s32. Note-se que li tem de ser adicionado a fronteira de ambos os subespacos (passo 6).

Para definir as fronteiras dos novos subespacos Fr(Ω+) e Fr(Ω−) e necessarioefectuar a classificacao dos segmentos de recta que estao na fronteira do subespacoFr(Ω) em relacao a linha l. A classificacao de cada segmento reduz-se a classificacaodos pontos terminais de cada segmento em relacao a l. Seja M um ponto qualquer dosegmento li; por exemplo, seja M o ponto medio de PQ (Figura 3.3(c)). A classificacaode qualquer ponto X na fronteira de Ω podera ser efectuada do seguinte modo:

3.2. ALGORITMO BSP E CONDICOES DE PARAGEM 41

l Ω-

Ω+

s1

s2

s4

s5

(a) (b)

s3

P

Q

li

Ω-

Ω+

s1

s2

s4

s51 s31

P

Q

s32

s52

li

s1

s2

s4

s51 s31

P

Q

s32

s52 M

G

C

u

(c)

A

B

D

E

F

Figura 3.3: Classificacao de pontos em relacao a linha BSP.

• Determina-se o vector unitario −→u aplicado em M com a direccao de l;

• Determina-se o vector unitario −→w perpendicular a Ω em M ;

• Determina-se o vector unitario −→v =−−→MX;

• Se o produto misto (−→v ×−→u ).−→w < 0, X pertence a Ω−;

• Se o produto misto (−→v ×−→u ).−→w > 0, X pertence a Ω+.

Na Figura 3.3(c), os pontos F,G pertencem a Ω+, enquanto que C, D,E per-tencem a Ω−. Se o produto misto (−→v × −→u ).−→w = 0, o ponto X pertence quer a Ω+

quer a Ω−. Isto acontece quando X ∈ li (e.g. A e B) e o ponto de interseccao entre afronteira de Ω e l.

Este processo aparentemente complexo e necessario, dado que Ω podera ser umsubespaco de qualquer plano Ax+By+Cz+D = 0 em R3 (nao necessariamente Z = 0).Se se tivesse somente considerado o subespaco Ω em R2, bastaria ter a equacao darecta ax + by + c = 0 para classificar qualquer ponto X = (x0, y0) relativamente aosegmento li. Se ax0 + by0 + c > 0, entao X ∈ Ω+. Se ax0 + by0 + c < 0, entao X ∈ Ω−.No caso de ax0 + by0 + c = 0, entao X pertence a fronteira dos dois subespacos (Ω− eΩ+) simultaneamente, uma vez que e um ponto do segmento li. Mas como se vera noproximo capıtulo, a classificacao de pontos dum plano em R3 em relacao a uma linharecta tem a vantagem de poder ser usada na amostragem de pontos de superfıciesimplıcitas.

A amostragem dos pontos da curva (passo 7) no segmento de bisseccao li e expli-cada mais adiante. Os passos 8 e 9 sao semelhantes aos passos 5 e 6, respectivamente.

3.2 Algoritmo BSP e condicoes de paragem

A particao binaria do subespaco Ω em dois subespacos mais pequenos Ω+ e Ω− edescrito atraves do seguinte algoritmo (BSP):


Algoritmo 3.2 BSP

INPUT:

(a) Ω: subespaco de R2

(b) ∆: distancia maxima admissıvel entre pontos

Begin

1. if (nenhum criterio de paragem e verdadeiro)

- l ← Determina a linha de bisseccao de Ω

- Subspaces(Ω, l, Ω−, Ω+);

- BSP(Ω−,∆);

- BSP(Ω+,∆);

End

A particao de Ω e recursiva visto que o algoritmo se chama a si proprio para Ω−e Ω+. O algoritmo BSP para se algum dos seguintes criterios de paragem se verificar:

- 1o criterio. Se nao existe nenhum segmento 1-dimensional da curva em Ω. Note-se que nao se esta aqui a considerar os pontos isolados porque o seu tratamentosera feito noutro algoritmo.

- 2o criterio. Este criterio tem que ver com o controlo da curvatura. Em primeirolugar, testa-se se as condicoes: d(P, Q) < τ (com τ = 3∆) e #(Fr(Ω) ∩ C) == 2sao verdadeiras ou nao. Se forem verdadeiras, determina-se um ponto M =l ∩ C entre P e Q, e calcula-se o angulo ∠(PM, MQ). Se este angulo foraproximadamente 180o, a particao do subespaco termina pois a variacao dacurvatura nao e significativa.

- 3o criterio. Se a distancia entre os dois pontos da curva mais afastados quepertencem a fronteira de Ω e menor que ∆, entao nao e necessario subdividirmais o subespaco.

3.3 Determinacao dos pontos da curva

Uma curva implıcita pode possuir uma ou mais componentes. Por exemplo, as curvasrepresentadas nas Figuras 3.12(b), (d), e (a) tem uma, duas e quatro componentes,respectivamente. Esta seccao e dedicada ao calculo dos pontos de interseccao entre

3.3. DETERMINACAO DOS PONTOS DA CURVA 43

cada segmento li e a curva C definida pela funcao real f (passo 7 do algoritmoSubspaces). Este calculo reduz-se a determinacao de zeros da funcao ao longo de li.Contudo, os algoritmos numericos de bisseccao, falsa posicao, Newton, Ridders e deBrent tem duas pre-condicoes que sao bastante limitativas no processo de amostragemduma curva implıcita, nomeadamente:

- Raiz unica. E assumido que so existe uma raiz num determinado intervalo [A,B],i.e. ∃1X ∈ [A,B] : f(X) = 0.

- Corolario de Bolzano. E assumido que f muda de sinal entre A e B, i.e.f(A).f(B) < 0.

Ora, acontece que a curva pode intersectar li em um ou mais pontos, ou seja,pode haver mais que um ponto da curva em li = AB como se mostra na Figura 3.3(b).A solucao obvia para contornar a primeira pre-condicao e subdividir li em pequenossegmentos, como se mostra no seguinte algoritmo:

Algoritmo 3.3 CurvePoints

INPUT:

(a) f : Rn → R

(b) li: segmento de bisseccao AB

(c) ε : erro numerico maximo admissıvel

OUTPUT:

(a) Ip: lista de pontos de interseccao do segmento li com C

Begin

1. n ← Determinar o numero de segmentos de li

2. Ip ← ∅

3. foreach i = 0, 1, . . . , n− 1

(a) Xi ← A + in.(B − A)

(b) Xi+1 ← A + i+1n

.(B − A)

(c) Zeros(f, Xi, Xi+1, Ip, ε)

4. retorna Ip

End


A subdivisao de li (passos 3(a) e 3(b) do Algoritmo 3.3) em pequenos segmentose necessaria para garantir a priori que o algoritmo possa determinar a maior parte doszeros da funcao, ou seja, os pontos da curva em li. A subdivisao de li em pequenossegmentos e adaptativa, na medida em que depende do comprimento d do segmento li.Seja D o comprimento maximo admissıvel de li que e o comprimento da diagonal dosubespaco inicial Ω, e N = 10, o numero maximo de segmentos em que se subdividea referida diagonal. Se d ≤ D/N , nao ha necessidade de subdividir mais li; casocontrario, subdivide-se li num numero de segmentos dado por

n = N .d

D(3.1)

Este numero de segmentos e calculado no passo 1 do Algoritmo 3.3). Depois, aplica-se a formula de interpolacao do metodo numerico TGFP a cada segmento XiXi+1

de modo a determinar os pontos li ∩ C, ou seja, todo o ponto X da curva C quesatisfaz f(X) = 0 ao longo de li. O metodo TGFP esta embutido no algoritmo Zeros

(passo 3) e, subsidiariamente, nos algoritmos ZeroSecante e ZeroExtremo que saodescritos de seguida. O metodo TGFP e —tanto quanto se julga saber— o primeirometodo numerico capaz de determinar ambos os tipos de zeros (zeros secantes e zerosextremos) no intervalo ou no segmento XiXi+1.

No entanto, como se vera de seguida, o metodo TGFP nao e usado directamentenos algoritmos ZeroSecante e ZeroExtremo. O que e realmente usado e a sua formulade interpolacao para calcular a proxima estimativa de cada zero. Mas isso requer queo segmento XiXi+1 em analise seja particionado em sub-segmentos de tal forma quecada segmento contenha um unico zero (ponto da curva).

3.3.1 Pontos com variacao de sinal

Qualquer ponto com variacao de sinal e um zero secante de f . Um zero secante e umzero de f no segmento XiXi+1 tal que f(Xi).f(Xi+1) < 0, ou seja, f(Xi) e f(Xi+1)tem sinais contrarios. Por exemplo, a curva definida por f(x, y) = y − x = 0 naFigura 3.4(a) tem uma so componente que e uma linha recta, aquem e alem da qualos respectivos potenciais da funcao tem sinais contrarios. Uma componente deste tipoe designada por componente com variacao de sinal, em que cada um dos seus pontosX satisfaz f(X) = 0. Toda a componente com variacao de sinal e amostrada atravesdo seguinte algoritmo:


Algoritmo 3.4 ZeroSecante

INPUT:

(a) f : Rn → R

(b) Xi, Xi+1 : extremos do segmento XiXi+1

(c) ε : erro numerico admissıvel

INPUT/OUTPUT:

(a) Ip: lista de pontos de interseccao do segmento XiXi+1 com f

Begin

1. X ← FormulaIteracao(Xi, Xi+1, f(Xi), f(Xi+1))

2. if (|f(X)| < ε) // X e zero ?

Ip ← Ip + X∆x ← ε

X− ← X − |∆x|. (Xi+1−Xi)||Xi+1−Xi||

X+ ← X + |∆x|. (Xi+1−Xi)||Xi+1−Xi||

if (X− ∈]Xi, X[) Zeros(f,Xi, X−, Ip, ε)

if (X+ ∈]X, Xi+1[) Zeros(f,X+, Xi+1, Ip, ε)

3. else // X nao e zero.

if (AssimpVertical(Xi, X, Xi+1, f(Xi), f(X), f(Xi+1))) return

if (||XiXi+1|| < ε)

Ip ← Ip + (Xi + Xi+1)/2return

if (f(Xi).f(X) < 0)

ZeroSecante(f, Xi, X, Ip, ε)

Zeros(f,X,Xi+1, Ip, ε)

else

ZeroSecante(f, X,Xi+1, Ip, ε)

Zeros(f,Xi, X, Ip, ε)

End


(a) y − x = 0

(b) (x2 + y2 − 3)2 = 0

(c) (y − x)(x2 + y2 − 3)2 = 0

Figura 3.4: Curvas com componentes de tipos distintos.

Como se vera mais adiante, este algoritmo e invocado pelo algoritmo Zeros,quando f(Xi).f(Xi+1) < 0. A proxima estimativa X e calculada atraves da formulade interpolacao do metodo TGFP (passo 1). Quando se detecta um zero da funcaoatraves do teste (|f(X)| < ε) (passo 2) e importante nao abandonar imediatamenteo intervalo onde o mesmo se encontra, pois podera haver mais zeros a direita ou aesquerda de X, ou seja, e necessario verificar se existem mais zeros nos intervalos[Xi, X[ e ]X,Xi+1], e nao nos intervalos [Xi, X] e [X,Xi+1], pois X ja e zero de f . Napratica, o calculo destes hipoteticos zeros faz-se nos intervalos [Xi, X−] e [X+, Xi+1],em que X− e X+ sao dois pontos muito proximos de X a distancia de ∆x = ε, aesquerda e a direita, respectivamente.


O calculo de X− e X+ e efectuado atraves da equacao vectorial da recta doseguinte modo:

X− = X − |∆x|. (Xi+1 −Xi)

||Xi+1 −Xi||

X+ = X + |∆x|. (Xi+1 −Xi)

||Xi+1 −Xi||

E importante garantir que o valor de ∆x e igual ou maior ao erro ε admissıvel;caso contrario, tendo em conta o teste |f(X)| < ε (passo 2), calcular-se-ia o mesmozero (a menos de ε).

O passo 3 do algoritmo ZeroSecante faz o tratamento dos segmentos XiX eXXi+1 quando a estimativa X nao e um zero. Se f(Xi).f(X) < 0 entao o zero secanteesta contido no segmento XiX; caso contrario, o zero pertence ao segmento XXi+1.Em qualquer dos casos, o algoritmo ZeroSecante chama-se a si proprio por forma acalcular uma nova estimativa X do zero. Note-se que o segmento que nao contemo zero nao e logo abandonado. E por isso que o algoritmo Zeros e invocado peloalgoritmo ZeroSecante, por forma a verificar se o segmento contem um numero parde zeros secantes ou um ou mais zeros extremos.

Repare-se que pode haver variacao de sinal sem que haja qualquer zero de fno segmento XiXi+1. Este caso pode ocorrer quando se esta perante uma assimptotavertical, como e o caso da curva y − 1/x = 0 da Figura 3.12(d). O teste a assimptotavertical e o primeiro teste do passo 3. Se este teste nao fosse efectuado, a assimptotaera identificada como fazendo parte da curva, o que nao e verdade.

O segundo teste do passo 3 diz respeito ao controlo da precisao no domınio de f .Quando o tamanho do segmento XiXi+1 e inferior a ε, entao assume-se que o pontomedio do segmento e um zero de f (ou ponto da curva).

3.3.2 Pontos sem variacao de sinal

Um ponto sem variacao de sinal e um zero extremo de f . Um zero extremo X e umaraiz de f no segmento XiXi+1 tal que f(Xi).f(Xi+1) > 0 e d(Xi).d(Xi+1) < 0, o quesignifica que a funcao f nao tem variacao de sinal no segmento XiXi+1, mas tem umextremo (maximo ou mınimo), pelo que e necessario depois verificar se esse extremoX e zero da funcao f .

O calculo destes zeros e importante para amostrar componentes da curva semvariacao de sinal como, por exemplo, a componente da curva f(x, y) = (x2+y2−3)2 = 0representada na (Figura 3.4(b)). Neste caso, a funcao e positiva aquem e alem dacomponente da curva. Qualquer componente sem variacao de sinal e pois amostradapelo seguinte algoritmo:


Algoritmo 3.5 ZeroExtremo

INPUT:

(a) f : Rn → R

(b) Xi, Xi+1 : extremos do segmento Xi, Xi+1


INPUT/OUTPUT:


Begin

1. X ← FormulaIteracao(Xi, Xi+1, f(Xi), f(Xi+1))

2. if (|f(X)| < ε) // X e zero ?

Ip ← Ip + X∆x ← ε

X− ← X − |∆x|. (Xi+1−Xi)||Xi+1−Xi||

X+ ← X + |∆x|. (Xi+1−Xi)||Xi+1−Xi||

if (X− ∈]Xi, X[) Zeros(f, Xi, X−, Ip, ε)

if (X+ ∈]X,Xi+1[) Zeros(f,X+, Xi+1, Ip, ε)

3. else // X nao e zero.

if (AssimpVertical(Xi, X,Xi+1, f(Xi), f(X), f(Xi+1))) return

if (||XiXi+1|| < ε)

X ← (Xi + Xi+1)/2

if (|f(X)| < ε)

Ip ← Ip + Xreturn

if (f(Xi) . f(X) < 0)

ZeroSecante(f, Xi, X, Ip, ε)

ZeroSecante(f, X, Xi+1, Ip, ε)

else

if (d(Xi) . d(X) < 0)

ZeroExtremo(f,Xi, X, Ip, ε)


Zeros(f, X, Xi+1, Ip, ε)

else

ZeroExtremo(f, X, Xi+1, Ip, ε)

Zeros(f, Xi, X, Ip, ε)

End

Portanto, o algoritmo ZeroExtremo so e invocado pelo algoritmo Zeros quandoexiste variacao da monotonia e nao ha variacao de sinal de f . Apos determinar a novaestimativa X (passo 1) atraves da formula de interpolacao, ha que verificar se X e umzero (passo 2). No caso de X ser zero, e necessario verificar se existem mais zeros nosintervalos [Xi, X−] e [X+, Xi+1] atraves do algoritmo Zeros, a semelhanca do que fazo algoritmo ZeroSecante.

Se X nao e zero (passo 3), ha que verificar se a condicao f(Xi).f(X) < 0 everdadeira ou nao. Se for verdadeira, entao existem dois zeros no segmento XiXi+1, oprimeiro no sub-segmento XiX e o segundo no sub-segmento XXi+1. Caso contrario,ha que verificar a existencia de um zero extremo num dos segmentos e chama-se oalgoritmo Zeros para o outro segmento, pois podera haver aı mais zeros.

Repare-se que o teste a assimptota tambem e efectuado no caso em que nao havariacao de sinal de f no segmento XiXi+1. Por exemplo, a curva (y − 1/x)2 = 0 daFigura 3.5(a) tem uma assimptota vertical em x = 0, a esquerda e a direita da quala funcao tem sempre o sinal positivo. A curva (y − tg(x))2 = 0 da Figura 3.5(b) temsempre sinal positivo, tendo assimptotas verticais em x = π/2 + kπ onde k ∈ Z. Aszonas mais claras da Figura 3.5(a) e (b) indicam a presenca de assimptotas verticais.

(a) (y − 1/x)2 = 0 (b) (y − tg(x))2 = 0

Figura 3.5: Curvas sem variacao de sinal.

O segundo teste do passo 3 tem a ver com o controlo da precisao no domınio def . Quando o tamanho do segmento XiXi+1 e inferior a ε, entao verifica-se se o ponto


medio do segmento e um zero de f (ou ponto da curva), nao havendo lugar a maissubdivisoes do segmento.

3.3.3 Pontos: algoritmo geral

Alem das componentes com (Figura 3.4(a)) e sem (Figura 3.4(b)) variacao de sinal,ha ainda um terceiro tipo de componentes: componentes mistas (Figura 3.4(c)). Umacomponente mista de uma curva possui segmentos com variacao de sinal e segmentossem variacao de sinal. Por exemplo, combinando as curvas descritas na Figura 3.4(a)e (b) obtem-se a curva definida por f(x, y) = (y−x)(x2 +y2−3)2 = 0 na Figura 3.4(c)que tem uma componente mista com tres segmentos com variacao de sinal (referentesa curva na Figura 3.4(a)) e dois segmentos sem variacao de sinal (referentes a curvana Figura 3.4(b)). Isto mostra que os algoritmos ZeroSecante e ZeroExtremo podemser ambos chamados para determinar os pontos de uma unica componente (mista) dacurva. Daı a necessidade de ter um algoritmo mais geral, chamado Zeros, capaz deamostrar os tres tipos de componentes:

Algoritmo 3.6 Zeros

INPUT:

(a) f : Rn → R

(b) Xi, Xi+1 : extremos do segmento XiXi+1


INPUT/OUTPUT:


Begin

1. if (f(Xi).f(Xi+1) < 0)

ZeroSecante(f,Xi, Xi+1, Ip, ε)

return

2. if ((f(Xi) > 0∧d(Xi) < 0∧d(Xi+1) > 0)∨(f(Xi) < 0∧d(Xi) > 0∧d(Xi+1) < 0))

ZeroExtremo(f,Xi, Xi+1, Ip, ε)

return

3. if (Cond1||Cond2||Cond3||Cond4)


// Subdivide XiXi+1 em XiX e XXi+1

X ← Xi+Xi+1

2

Zeros(f, Xi, X, Ip, ε);

Zeros(f, X, Xi+1, Ip, ε);

return

4. if (Cond5||Cond6)

Ip ← Ip + Xi+1∆x ← ε

X ← Xi+1 − |∆x|. (Xi+1−Xi)||Xi+1−Xi||

if (X ∈]Xi, Xi+1[) Zeros(f, Xi, X, Ip, ε)

return

5. if (Cond7||Cond8)

Ip ← Ip + Xi∆x ← ε

X ← Xi + |∆x|. (Xi+1−Xi)||Xi+1−Xi||

if (X ∈]Xi, Xi+1[) Zeros(f, X, Xi+1, Ip, ε)

return

End

No Algoritmo 3.3, a divisao de li em pequenos segmentos facilita a procurados zeros da funcao, mas nao garante que se determinem todos os pontos da curvapertencentes a li. Por sua vez, o Algoritmo 3.6 permite, na pratica, determinartodos os pontos da curva a menos de ε usando a seguinte estrategia: subdivisao eescansao(sondagem) de ambos os sub-segmentos XiX e XXi+1 de XiXi+1, sendo X aproxima estimativa do zero (ponto da curva).

Como ja foi referido acima, os algoritmos ZeroSecante e ZeroExtremo fazemuso da formula de interpolacao do metodo TGFP. Os algoritmos ZeroSecante eZeroExtremo aparecem nos passos (1) e (2) do algoritmo Zeros, respectivamente.Todos os zeros sao calculados neste dois primeiros passos. No entanto, ha situacoes emque e necessario subdividir mais o intervalo XiXi+1 por forma a isolar cada zero numsub-intervalo, o que e feito pela analise das condicoes 1-8 ilustradas nas Figuras 3.6 e3.7. Esta subdivisao do intervalo XiXi+1 pode acontecer em duas situacoes:

- Quando a funcao f nao e aproximadamente linear ao longo do segmento XiXi+1,como se mostra na Figura 3.6. Este problema e resolvido no passo (3).


Xi Xi+1

Cond2

Xi Xi+1

Cond3

Xi Xi+1

Cond4

Xi Xi+1

Cond1

(a) (b) (c) (d)

Figura 3.6: Possibilidade de haver mais zeros em XiXi+1.

Xi Xi+1

Cond5

Xi Xi+1

Cond6

Xi

Xi+1

Cond7

Xi

Xi+1

Cond8

(a) (b) (c) (d)

Figura 3.7: Existencia de outro zero no interior de XiXi+1 quando f(Xi).f(Xi+1) = 0.

- Quando a funcao f tem um zero num dos extremos do segmento XiXi+1, comose mostra na Figura 3.7. Este problema e resolvido nos passo (4) e (5).

Na primeira situacao, ha que ser particularmente cuidadoso quando f(Xi+1) ef(Xi) tem sinais identicos. Seja ∆f = f(Xi+1) − f(Xi) a diferenca dos valores de fnos pontos extremos de XiXi+1. Com se ilustra na Figura 3.6, quatro casos podemocorrer:

- Cond1 = (f(Xi) > 0) ∧ (d(Xi) < 0) ∧ (d(Xi+1) < 0) ∧ (∆f > 0). Neste caso opotencial da funcao oscila em XiXi+1, pelo que ha a possibilidade de existiremum ou mais zeros no segmento. Subdivide-se entao XiXi+1 pelo ponto medio Xe invoca-se o algoritmo Zeros para os sub-segmentos XiX e XXi+1.

- Cond2 = (f(Xi) > 0) ∧ (d(Xi) > 0) ∧ (d(Xi+1) > 0) ∧ (∆f < 0). Neste caso opotencial da funcao tambem oscila em XiXi+1, mas agora as diferencas finitasnos pontos extremos sao positivas. Assim, como no caso anterior, subdivide-seXiXi+1 pelo ponto medio X e invoca-se o algoritmo Zeros para os dois sub-segmentos.

- Cond3 = (f(Xi) < 0) ∧ (d(Xi) < 0) ∧ (d(Xi+1) < 0) ∧ (∆f > 0). Caso simetricoda condicao Cond1, com f(Xi) < 0 e f(Xi+1) < 0.

- Cond4 = (f(Xi) < 0) ∧ (d(Xi) > 0) ∧ (d(Xi+1) > 0) ∧ (∆f < 0). Caso simetricoda condicao Cond2, com f(Xi) < 0 e f(Xi+1) < 0.

3.4. DETERMINACAO DE PEQUENAS COMPONENTES 53

A segunda situacao ocorre quando um dos extremos de XiXi+1 e zero da funcao,i.e., f(Xi).f(Xi+1) = 0. A Figura 3.7 ilustra esta situacao, mostrando que alguns zerospodem ser indetectaveis, a nao ser que o algoritmo seja capaz de tratar as seguintescondicoes:

- Cond5 = (f(Xi) < 0) ∧ (f(Xi+1) = 0) ∧ (d(Xi+1) < 0). Neste caso, a funcaomuda de sinal em [Xi, Xi+1]. Assim, actualiza-se X por uma nova estimativa

dada por X ← Xi+1 − |∆x|. (Xi+1−Xi)||Xi+1−Xi|| a esquerda de Xi+1. Depois, invoca-se o

algoritmo Zeros para o sub-segmento XiX.

- Cond6 = (f(Xi) > 0) ∧ (f(Xi+1) = 0) ∧ (d(Xi+1) > 0). Este caso e semelhanteao anterior, com a diferenca de que d(Xi+1) > 0.

- Cond7 = (f(Xi) = 0) ∧ (f(Xi+1) < 0) ∧ (d(Xi) > 0). A funcao e agora nulaem Xi. Assim, substitui-se X por uma nova estimativa dada por X ← Xi +|∆x|. (Xi+1−Xi)

||Xi+1−Xi|| a direita de Xi. De seguida, invoca-se o algoritmo Zeros para o

sub-segmento XXi+1.

- Cond8 = (f(Xi) = 0) ∧ (f(Xi+1) > 0) ∧ (d(Xi) < 0). Este caso e semelhante aoanterior, com a diferenca de que d(Xi) < 0.

3.4 Determinacao de pontos isolados e outras com-

ponentes de tamanho reduzido

3.4.1 Pontos isolados

Um ponto isolado X e uma componente singular sem variacao de sinal de uma curva.Como tal, o sinal da funcao nao se altera numa pequena vizinhanca de qualquer pontoisolado. Apesar disso, um ponto isolado e tambem uma raiz de f , i.e. f(X) = 0. Istosignifica que um ponto isolado de uma curva pode ser visto como um mınimo localabsoluto nulo de f pertencente a curva.

A ideia subjacente ao calculo de qualquer ponto isolado e entao verificar se f emonotona ao longo de algum segmento AB da fronteira de cada subespaco associadoa cada no terminal da arvore BSP. O estudo da monotonia de f e normalmente feitaatraves das derivadas, mas, tendo em conta que o metodo TGFP nao usa derivadas,recorreu-se ao uso das diferencas finitas direccionais. De facto, o metodo TGFP podeser usado para determinar uma sequencia de mınimos absolutos que tende para umzero de f que e geometricamente um ponto isolado. Isto e possıvel porque o metodoTGFP pode ser usado para determinar zeros e extremos (mınimos e maximos) de umafuncao real.


A Figura 3.8 ilustra a determinacao de uma sequencia de mınimos convergentepara um zero ou um ponto isolado P dentro de um subespaco terminal Ω. Em primeirolugar, determina-se um mınimo local absoluto P0 de f num dos segmentos AB deFr(Ω) com uma precisao numerica ε (Figura 3.8(a)). Em segundo lugar, calcula-se umsegmento de recta em Ω que e perpendicular a AB em P0. Este segmento de rectapassa a ser o segmento AB, sobre o qual se calcula o novo mınimo local absoluto P1

(Figura 3.8(b)). Em terceiro lugar, determina-se o novo segmento AB perpendicularao anterior em P1, assim como o novo mınimo absoluto P2 sobre AB (Figura 3.8(c)).Este processo e repetido ate que Pi (i = 0, 1, . . . , n) seja aproximadamente um zerode f a menos de ε, como se descreve no seguinte algoritmo:

Algoritmo 3.7 IsolatedPoint

INPUT:

(a) Ω: subespaco terminal

(b) Si: um segmento da fronteira de Ω

(c) ρ : nıvel maximo de recursividade

(d) ε : erro numerico maximo admissıvel

OUTPUT:

(a) Pi: um ponto isolado

Begin

1. if (ρ == 0) return NULL

2. if (Si == NULL)

Pi ← NULL

nsi ← 0

while (nsi < #Fr(Ω)) ∧ (Pi == NULL)

Si ← segmento numero nsi de Fr(Ω)

Pi ← mınimo local no segmento Si de C, usando o metodo TGFP

nsi ← nsi + 1

else

Pi ← mınimo local do segmento Si, usando o metodo TGFP

3. if (Pi == NULL) return NULL


4. if (|f(Pi)| ≤ ε) return Pi

5. Si+1 ← segmento perpendicular a Si em Pi dentro de Ω

6. return IsolatedPoint(Ω, Si+1, ρ− 1, ε)

End

Repare-se que o algoritmo IsolatedPoint e recursivo (passo 6), de tal modo quede cada vez que o metodo e invocado e calculado um novo segmento AB (passo 5),assim como um novo mınimo absoluto de f em AB (passo 2). Este algoritmo devolveum ponto isolado P = Pi se ele existir (passo 4) ou NULL no caso contrario (passo 3).

(a) B

A

(b)

P2 B A

(c) (d)

P1

P0 B A

Pi

B

A

P

Figura 3.8: Determinacao de um ponto isolado num subespaco Ω.

A prova de que este processo e convergente para o ponto isolado P e dada peloseguinte teorema:

Teorema 3.1 Dado um subespaco Ω com um ponto isolado P da curva C = (x, y) ∈Ω ⊂ R2 : f(x, y) = 0, entao existe uma sequencia de pontos Pn∞n=0 em Ω queconverge para P . ¤

Demonstracao. A estrategia consiste em demonstrar que existe uma sucessao demınimos absolutos mi∞i=0, com mi = |f(Pi)|, que e decrescente e limitada. Seja


S0 um segmento de Fr(Ω). Usando o metodo TGFP, determina-se um ponto P0 emS0 onde f tem um mınimo absoluto. O proximo ponto P1 e calculado no segmento S1

que e perpendicular a S0 em P0. De um modo geral, o proximo ponto Pi e calculadodo seguinte modo:

Pi = X ∈ Si : |f(X)| < |f(Y )|, ∀Y ∈ Si com Si⊥Si−1 em Pi−1e Pi 6= Pi−1

Se f(Pi) = 0 entao o ponto isolado P foi encontrado. Caso contrario, pordefinicao de Pi, pode concluir-se que |f(Pi)| < |f(Pi−1)|, uma vez que a distancia dosvarios Pi a P vai decrescendo, i.e. d(Pi, P )| < d(Pi−1, P ). Isto prova que a sequenciami∞i=0 e monotona decrescente. Alem disso, f(Pi) ∈ [0, |f(P0)|], o que significaque mi∞i=0 e limitada. Em consequencia, pelo Teorema da Sucessao Monotona eLimitada, pode concluir-se que mi∞i=0 e convergente para 0, o que e equivalente adizer que Pn∞n=0 converge para P . ¥

3.4.2 Teste do ponto isolado

Para verificar se um dado ponto P e um ponto isolado determina-se a interseccao dacurva com uma vizinhanca quadrada inscrita numa pequena circunferencia centradaem P , como se ilustra na Figura 3.9. Se nao existe interseccao, assume-se que P e umponto isolado; caso contrario, P nao e ponto isolado. O algoritmo de teste do pontoisolado e o seguinte:

Algoritmo 3.8 IsolatedPointChecking

INPUT:

Ω : subespaco de R2

P : possıvel ponto isolado

f : funcao real que descreve a curva

r : raio da vizinhanca circular centrada em P

OUTPUT:

devolve true ou false, conforme o ponto seja isolado ou nao

Begin

(1) −→w ← vector ortogonal a Ω em P ;

(2) −→v ← vector ortogonal a −→w em Ω;


(3) −→u ← −→w ×−→v ;

(4) −→u ← r.−→u||−→u || ;

(5) −→v ← r.−→v||−→v || ;

(6) X0 ← P +−→u ;

(7) X1 ← P +−→v ;

(8) X2 ← P −−→u ;

(9) X3 ← P −−→v ;

(10) return (X0X1 ∩ C) ∨ (X1X2 ∩ C) ∨ (X2X3 ∩ C) ∨ (X3X0 ∩ C)

End

E importante calcular o vector −→w (vector ortogonal a Ω em R3), pois ha aintencao de aplicar o algoritmo descrito neste capıtulo na amostragem de superfıcies,como se vera no proximo capıtulo. Daı que o calculo do vector −→w seja feito paraqualquer plano em R3. Neste capıtulo, o subespaco Ω esta contido em R2 (ou Z = 0),pelo que o vector −→w e dado por (0, 0, 1).

P

(a)

v

u

w

(b)

P v

u X0

X1

X2

X3

(c)

X0 X1

X2

X3

Figura 3.9: Interseccao da curva com o quadrado inscrito na vizinhanca circular de P .

3.4.3 Sobre-particao dos subespacos terminais

O Teorema 3.1 assume que existe somente um ponto (isolado) da curva em Ω. Masraramente isso acontece pois, apos terminar o processo de subdivisao BSP, a maiorparte dos subespacos terminais contem segmentos 1-dimensionais da curva. Estessubespacos terminais contem eventualmente um ou mais pontos isolados ou ate mesmocomponentes 1-dimensionais que nao foram detectadas no processo de subdivisaopor serem demasiado pequenas. Assim, torna-se necessario prosseguir o processode subdivisao BSP por mais alguns passos de modo a que cada ponto isolado ousegmento da curva fique num unico subespaco. A Figura 3.10 apresenta os tres casosfundamentais em que e necessario subdividir ainda mais o subespaco terminal Ω:


(a) (b) (c)

Ω P

Q

I

Ω Ω

J

I

Figura 3.10: Tres casos onde e necessario particionar o subespaco.

1. A existencia de pelo menos um ponto isolado em Ω, sendo Ω atravessado porum ou mais segmentos 1-dimensionais da curva (Figura 3.10 (a));

2. A existencia de pelo menos uma componente 1-dimensional da curva em Ω(Figura 3.10 (b));

3. A existencia de pelo menos dois pontos isolados em Ω (Figura 3.10 (c)).

O Algoritmo 3.9 descreve como e efectuada esta subdivisao adicional de cadasubespaco terminal, de modo a processar as componentes 0-dimensionais (pontosisolados) e 1-dimensionais de tamanho reduzido.

Algoritmo 3.9 OverBSP

INPUT:

(a) Ω: subespaco terminal

(b) ρ : nıvel maximo de recursividade


(d) ∆: distancia maxima entre pontos

INPUT/OUTPUT:

(a) Lisol: lista de pontos isolados

Begin

(1) if (ρ == 0) return

(2) if @Si ⊂ Fr(Ω) com pelo menos 1 mınimo absoluto nao-nulo Pi

return


(3) S⊥ ← segmento perpendicular a Si em Pi

(4) Calcular a lista Lmin de mınimos absolutos de f em S⊥

(5) if (#Lmin ≥ 2)

(a) if existe pelo menos um mınimo absoluto nao-nulo

Subspaces(S⊥, Ω, Ω1, Ω2)

OverBSP(Ω−, ρ-1, Lisol, ε)

OverBSP(Ω+, ρ-1, Lisol, ε)

(b) else

BSP(Ω, ∆/10)

(6) else //#Lmin = 1

- if o mınimo absoluto e nao-nulo

Lisol ← Lisol+IsolatedPoint(Ω, S⊥, ρ, ε)

End

De modo a determinar qualquer componente da curva em Ω, que nao tenha sidodetectada anteriormente pelo processo BSP, e necessario determinar um segmentoSi da fronteira de Ω com pelo menos um mınimo absoluto nao-nulo Pi (passo 2 doAlgoritmo 3.9). (Recorde-se que os mınimos absolutos sao determinados pelo metodoTGFP.) A existencia de um ou mais mınimos absolutos nao-nulos num segmentofronteira de Ω dao uma boa indicacao da presenca de uma ou mais componentesda curva em Ω.

Para saber se existe de facto uma ou mais componentes de tamanho reduzidoem Ω, ha que verificar se o potencial da funcao em valor absoluto decresce ate zero nosentido do interior de Ω. Assim, so se consideram os pontos Pi com mınimos absolutosnao-nulos que obedecem a condicao f(Pi).d−→v (Pi) < 0, em que d−→v (Pi) e a derivadadiscreta direccional segundo um dado vector −→v perpendicular ao segmento da fronteirade Ω em Pi no sentido do interior de Ω. Portanto, so existe alguma componente nointerior de Ω em duas situacoes:

- Se f(Pi) > 0 e d−→v (Pi) < 0, entao o potencial da funcao decresce no sentido dointerior de Ω, o que indica que ha grande probabilidade de existir um zero nointerior de Ω.

- Se f(Pi) < 0 e d−→v (Pi) > 0, entao o potencial da funcao cresce no sentido dointerior de Ω, o que indica que ha grande probabilidade de existir um zero nointerior de Ω.


Em qualquer destas situacoes, a procura de componentes da curva em Ω continuanos passos subsequentes do Algoritmo 3.9. Assim, apos determinar o mınimo Pi,calcula-se o segmento S⊥ perpendicular a Si em Pi (passo 3), assim como os mınimosabsolutos de f em S⊥ (passo 4). Se o numero de mınimos e pelo menos dois, entaoe necessario subdividir Ω (passo 5). Se existe pelo menos um mınimo absoluto nao-nulo, entao subdivide-se Ω em Ω− e Ω+ atraves de S⊥. Depois, o algoritmo OverBSP

e invocado para Ω− e Ω+ (casos da Figura 3.10(a) e (c)). No caso contrario, invoca-seo Algoritmo BSP para Ω com uma precisao de ∆/10 (caso da Figura 3.10(b)). Nocaso em que exista um so mınimo e que nao e nao-nulo, entao invoca-se o algoritmoIsolatedPoint (passo 6).

Note-se que a ideia geral do algoritmo e que, no final da subdivisao bintree, cadasubespaco terminal tenha um unico segmento da curva ou um unico ponto isolado.So assim as condicoes do Teorema 3.1 sao satisfeitas, sendo entao possıvel usar oalgoritmo IsolatedPoint ilustrado na Figura 3.8. Claro que e necessario distinguirentre um ponto isolado e um segmento de curva, o que e feito atraves do algoritmoIsolatedPointChecking. Note-se que o facto de se tomarem os pontos Pi onde ftem um mınimo absoluto no sentido do interior de Ω (passo 2) permite abandonarmuito cedo a procura dos pontos isolados no subespaco que nao verifique a condicaof(Pi).d−→v (Pi) < 0.

3.5 Determinacao de pontos de auto-interseccao

O algoritmo BSP nao bloqueia nem quebra na presenca de singularidades da curva, taiscomo apices e auto-interseccoes. Isto deve-se ao facto de a amostragem dos pontosser feita a custa das diferencas finitas, e nao de derivadas. Alem disso, como sedemonstra pelo Teorema 3.2, a subdivisao bintree de um subespaco Ω converge paraas auto-interseccoes da curva, se elas existirem.

Teorema 3.2 Seja Ω um dado subespaco que contem um ponto de auto-interseccao.O processo de particao binaria do espaco converge para um pequeno subespaco quecontem a auto-interseccao. ¤

Demonstracao. Seja Ω um subespaco que contem um ponto de auto-interseccao, e,Ω− e Ω+ os dois subespacos nao-nulos que resultam da particao de Ω = Ω− ∪ Ω+, oque implica que Ω− ⊂ Ω e Ω+ ⊂ Ω. Sem perda de generalidade, assume-se que Ω− eo subespaco que contem o ponto de auto-interseccao apos a particao de Ω, pelo que esuficiente provar que a area de Ω− tende para zero. Isto e o mesmo que dizer que Ω−converge para o ponto de auto-interseccao.

Seja A1(Ω−) = A(Ω)α1

a area do subespaco Ω− que contem o ponto de auto-interseccao, sendo α1 > 1. Note-se que α1 tem de ser maior que 1, pois Ω− ⊂ Ω eΩ+ 6= ∅. Caso contrario, se α1 ≤ 1, entao A1(Ω−) ≥ A(Ω), o que e absurdo pois

3.6. VISUALIZACAO DA CURVA 61

Ω− ⊂ Ω. O valor de α1 e um valor muito proximo de 2 uma vez que o subespaco Ωe particionado pela mediatriz dos pontos mais afastados, i.e. a particao de Ω tende aser pela metade.

Assim, a particao recursiva do subespaco que contem o ponto de auto-interseccaoconduz a obtencao da seguinte sequencia de areas: A1(Ω−) = A(Ω)

α1, A2(Ω−) = A1(Ω)

α2=

A(Ω)α1×α2

,..., An(Ω−) = A(Ω)α1×α2×...×αn

, com α1, α2, . . . , αn > 1. Seja α = minα1, α2, ..., αno mınimo de todos αi (i = 1, 2, ..., n). Portanto, A(Ω)

α1×α2×...×αn< A(Ω)

α×α×...×α= A(Ω)

αn , o que

implica que a area A(Ω)αn do subespaco que contem a singularidade converge para zero

a medida que o numero de subdivisoes n tende para +∞, ou seja:

limn→∞

A(Ω)

αn= 0

Fica assim demonstrado que o processo de particao converge para a singularidade. ¥

Este teorema mostra que o algoritmo BSP pode ser explorado para determinaros pontos de auto-interseccao de uma curva. Para isso, comeca-se por identificar ossubespacos terminais da arvore BSP que tem tres ou mais pontos da curva na fronteirade Ω, i.e. #(Fr(Ω) ∩ C) ≥ 3. Note-se que esta condicao nao e valida para subespacosintermedios da arvore BSP. Por exemplo, na Figura 3.2(b), existem quatro pontosna fronteira de Ω, mas nao existe nenhuma auto-interseccao em Ω. A procura dassingularidades so deve ocorrer nos subespacos terminais. De facto, todo o subespacoterminal Ω e relativamente pequeno, uma vez que satisfaz a condicao d(P,Q) < ∆,i.e. a distancia entre os dois pontos P, Q ∈ Fr(Ω) ∩ C mais afastados e menor que∆. Assim se a condicao #(Fr(Ω) ∩ C) ≥ 3 for verdadeira, assume-se que o subespacocontem uma auto-interseccao, sendo esta aproximada pelo ponto medio de PQ.

Como nota marginal, este metodo podera tambem ser usado para determinar ospontos de interseccao entre duas curvas planares C1 e C2. Sejam f e g as funcoes reaisque descrevem as curvas C1 e C2, respectivamente. Entao os pontos de interseccaodas duas curvas sao pontos de auto-interseccao de uma outra curva C descrita peloproduto das funcoes f × g. Os pontos de interseccao de C1 e C2 sao todos os pontosde auto-interseccao de C que sao pontos de C1 e C2.

3.6 Visualizacao da Curva

Depois de executados os algoritmos BSP e OverBSP que permitem construir a estruturaBSP, assim como determinar os pontos de uma curva, e necessario efectuar a visual-izacao da curva. Representar uma curva resume-se a determinar uma aproximacaopoligonal dos pontos encontrados no processo de particao. Para isso, nao e necessarioimpor qualquer ordem especial na representacao dos pontos. Em geral, e so necessariodesenhar um segmento de recta entre dois pontos de um no terminal associado a umsubespaco. No caso de o subespaco conter um ponto de auto-interseccao, e necessario


desenhar um segmento de recta entre o ponto de auto-interseccao e cada ponto dafronteira do subespaco.

Interessa aqui referir que somente os nos terminais da estrutura BSP interessamno processo de visualizacao, pois os nos intermedios servem somente para controlar oprocesso de subdivisao do subespaco. Assim, torna-se importante ter os nos terminaisligados entre si, de modo a que o processo de visualizacao da curva seja equivalente apercorrer uma lista ligada. E por isto que cada no BSP tem um apontador next. Ofacto de se visitarem somente os nos terminais atraves da lista ligada permite ter umganho temporal aproximado de 50% na representacao da curva em relacao ao tempodispendido para atravessar toda a estrutura BSP, como se demonstra no seguinteteorema:

Teorema 3.3 Representar uma curva implıcita utilizando somente os nos terminaisda arvore BSP permite um ganho aproximado de 50% em relacao a travessia de todaa estrutura BSP. ¤

Demonstracao. Seja Θ a arvore BSP com Ni nos apos particionar i = 0, 1, . . . , nsubespacos, sendo Mi o numero de nos terminais de Θ. No inıcio, Θ tem somenteum no (no raiz), i.e. N0 = M0 = 1. Depois de particionar Ω em dois subespacos,tem-se N1 = 3 nos e M1 = 2 nos terminais. Por recursao, apos particionar o i-esimosubespaco, a arvore BSP tem Ni = Ni−1 +2 nos e Mi = Mi−1 +1 nos terminais. Destemodo, pode concluir-se que

limi→∞

Mi

Ni

= limi→∞

Mi−1 + 1

Ni−1 + 2= lim

i→∞Mi−2 + 1 + 1

Ni−2 + 2 + 2= . . . = lim

i→∞M0 + i

N0 + 2i=

1

2

ou seja o numero de nos terminais e cerca de metade do numero de nos totais daestrutura BSP. ¥

O algoritmo geral para representar uma curva e o seguinte:

Algoritmo 3.10 ImplicitCurve

INPUT:

(a) C: curva definida por f(x, y) = 0

(b) Ω: subespaco rectangular de R2

(c) ∆: distancia maxima admissıvel entre pontos

Begin

1. BSP(Ω,∆);

3.7. COMPARACAO E RESULTADOS EXPERIMENTAIS 63

2. overBSP(Ω,∆/10);

3. foreach subespaco terminal Ω

- if (existe auto-interseccao)

self ← Determina ponto de auto-interseccao

foreach ponto da curva P ∈ Fr(Ω)drawlinesegment(self,P )

- else (#F=2)

P ← Devolve o primeiro ponto da curva F

Q ← Devolve o segundo ponto da curva F

drawlinesegment(P ,Q)

4. Desenha os pontos isolados

End

Note-se que a computacao dos pontos isolados e efectuado atraves do algoritmooverBSP.

3.7 Comparacao e resultados experimentais

3.7.1 Comparacao geral

De seguida, apresenta-se uma comparacao do algoritmo proposto neste capıtulo comoutros algoritmos existentes na literatura, assim como com os algoritmos embutidosem alguns softwares comerciais.

Esta comparacao e feita em termos da resolucao dos seguintes problemas:

- Singularidades. Muitos algoritmos de perseguicao quebram no calculo das sin-gularidades (veja-se, por exemplo, Chandler [13]). Moeller e Yagel propuseramum algoritmo mais geral que e capaz de lidar com curvas que possuem pontos debifurcacao atraves da analise da mudanca de sinal das derivadas parciais numavizinhanca rectangular [34]. Mas, o algoritmo de Moeller e Yagel quebra noutrassingularidades (por exemplo, apices) que pertencem ao domınio da funcao, masnao ao domınio das derivadas parciais. Por exemplo, nao e capaz de representaras singularidades (0, 2), (2, 0), (0,−2), e (−2, 0) da curva |x| + |y| − 2 = 0 (naforma de losango), o que se deve ao facto de as derivadas parciais nao existiremnos pontos anteriores.


(a) ((x2 + y2)2 − 2x2 + 2y2)((x + 2.5)2 + (y + 2.3)2 − 0.3)2((x− 1)2 + (y − 1.5)2) = 0

(b) (x− 2.3 + (y − 1.5)2)((x + 1.3)2 + (y + 1.5)2)((x + 1.6)2 + (y − 1.8)2) = 0

(c) (x3 + y)((x− 2)2 + (y + 1)2)((x + 3)2 + (y − 1.5)2 − 0.5)2 = 0

(d) (y3 − |x|)((x + 1)2 + (y + 1)2) = 0

Figura 3.11: Curvas algebricas: (esquerda) processo de subdivisao BSP e (a direita)poligonizacao. Regioes verdes e roxas e onde f > 0 e f < 0, respectivamente.


(a) (x2 + 4y2 + y3 − 2 + yx2+0.45 − 1)((x− 2.2)2 + (y− 2.2)2)((x + 2.8)2 + (y− 2.2)2 − 0.3)2 = 0

(b) (x− 1y2+1 )(y − 1

x2+1 ) = 0

(c) (x2 − yy2x2+1 )(y − 1

x2+1 ) = 0

(d) y − 1x = 0

Figura 3.12: Curvas racionais: (esquerda) processo de subdivisao BSP e (a direita)poligonizacao. Regioes verdes e roxas e onde f > 0 e f < 0, respectivamente.


(a) sin3(2x) + 4 sin3(y)− 3 sin(2x) sin(y) = 0

(b) (y2 − log2(x2 + 1))(x2 − log2(y2 + 1)) = 0

(c) (sin(x + y)− y2)(log(x2 + 1)− y2) = 0

(d) (tan(x)− sin(y))(|y|+ |x| − 1)(|y2 |+ |x| − 1) = 0

Figura 3.13: Curvas transcendentes.


Os algoritmos de subdivisao do espaco tem a vantagem de nao dependeremda diferenciabilidade (e da regularidade requerida pelos algoritmos baseados naconversao) da funcao que descreve uma determinada curva implıcita. Comoconsequencia, estes algoritmos nao quebram nos apices e nas auto-interseccoes,pois a amostragem dos pontos da curva pode ser feita atraves de um metodonumerico desprovido de derivadas (por exemplo, o metodo de bisseccao). Assim,em princıpio, os algoritmos de subdivisao sao mais robustos que os algoritmosde perseguicao.

Por outro lado, apos varias sessoes de teste de softwares matematicos (porexemplo, Maple, MathLab e Mathematica), constatou-se que estes tem difi-culdades em lidar correctamente com singularidades, em particular com auto-interseccoes. Por exemplo, a curva (y2 − log2(x2 + 1))(x2 − log2(y2 + 1)) = 0—que aparece correctamente representada na Figura 3.13(b)— tem um ponto deauto-interseccao em (0, 0). Porem, nao foi possıvel obter a mesma representacaousando o Maple, pois o ponto de auto-interseccao (0, 0) nao aparece representado.Na verdade, o que existe e um atalho entre cada dois ramos da curva, dandoa ilusao de que a curva tem quatro componentes em vez de uma so. No casodo Mathematica, o resultado obtido na representacao da curva foi ainda maisbizarro.

- Pontos isolados. A determinacao de pontos isolados atraves dos algoritmos deperseguicao nao tem solucao conhecida, ate porque se assume que a curva esempre contınua e 1-dimensional, como e inerente a ideia de determinacao doproximo ponto atraves do ponto actual da curva.

Os algoritmos de conversao baseiam-se tambem na mesma ideia, uma vez quetoda a curva parametrica pode ser vista como a imagem do intervalo [0, 1].

Na generalidade, os software comerciais acima mencionados nao conseguem de-terminar pontos isolados. Por exemplo, o ponto isolado da curva (x3 + y)((x−2)2 + (y + 1)2)((x + 3)2 + (y − 1.5)2 − 0.5)2 ilustrada na Figura 3.11(c) nao ecomputavel por estes softwares.

Em geral, os algoritmos de subdivisao do espaco tambem nao determinam pontosisolados. As tecnicas baseadas na aritmetica intervalar conseguem detectar aexistencia de qualquer componente da curva num dado subespaco, incluindopontos isolados. Mas a determinacao de um ponto isolado requer a subdivisaorecursiva do subespaco que o contem ate que o subespaco final seja suficiente-mente pequeno. Nestas circunstancias, pode considerar-se que o ponto isolado e obaricentro dos vertices do subespaco. Pelo contrario, o algoritmo BSP determinatodos os pontos isolados (veja-se o exemplo da Figura 3.11(c)). No algoritmoBSP, qualquer ponto isolado e detectado e determinado pelo algoritmo TGFPatraves do calculo de uma sequencia de mınimos, sem ser necessario efectuarqualquer subdivisao do subespaco terminal que o contem, o que resulta numganho de desempenho.

- Componentes 1-dimensionais sem variacao de sinal. Os algoritmos de persegui-


cao, assim como os algoritmos baseados na conversao da representacao, precisamde um primeiro ponto em cada componente da curva de modo a poder representa-la.

Os algoritmos baseados na subdivisao espacial nao precisam de um primeiroponto em cada componente da curva porque, ao subdividir o espaco, a probabil-idade de encontrar uma componente da curva e bastante elevada, com excepcaode pontos isolados, componentes 1-dimensionais muito pequenas e componentes1-dimensionais sem variacao de sinal.

Como se referiu anteriormente, a aritmetica intervalar permite localizar com-ponentes da curva, mas so trabalha com subdivisoes regulares. Alem disso, esempre necessario usar um algoritmo numerico para determinar os pontos dacurva apos a subdivisao (veja-se, por exemplo, a subdivisao quadtree em [29]).

O problema e que, com a excepcao dos algoritmos numericos GFP, TGFPe NGFP descritos nesta tese, nao ha nenhum algoritmo numerico capaz dedeterminar os pontos que resultam da interseccao de um segmento com umacomponente da curva sem variacao de sinal. Quer dizer, o facto de a funcaoser sempre positiva (respectivamente, negativa) aquem e alem da curva implicaque os metodos numericos classicos nao conseguem determinar os pontos dumacomponente sem variacao se sinal, uma vez que todos eles se baseiam no Corolariode Bolzano. Isto explica a razao porque a maioria dos algoritmos apresentadosna literatura e implementados em pacotes de software comerciais nao determi-nam componentes 1-dimensionais sem variacao de sinal, como acontece com ascircunferencias apresentadas na Figura 3.11(a) e (c) e Figura 3.12(a).

- Descontinuidades e controlo do domınio da funcao. O controlo sobre o domıniode funcao e um assunto importante na amostragem de curvas implıcitas, umavez que para determinar um ponto de uma curva e necessario avaliar variasvezes a funcao atraves de um metodo numerico. Sem o controlo efectivo sobreo domınio da funcao, o algoritmo pode quebrar a qualquer momento. Portanto,ha que impor restricoes apropriadas sobre o domınio da funcao. Por exemplo,o calculo de f(x, y) = y − √x = 0 requer que a restricao x ≥ 0 seja satisfeita.Outro exemplo ocorre com a funcao f(x, y) = y − 1

x2 = 0 Figura 1.5(b), em queo calculo de f esta restrito a x 6= 0. Para a subdivisao quadtree controlada poraritmetica intervalar, isto ainda e mais importante porque a assimptota verticalx = 0 desta funcao sera considerada como uma componente da curva. Istoacontece porque as operacoes aritmeticas intervalares baseiam-se no Teorema doValor Intermedio. Este facto, como se pode constatar na Figura 1.5(b), conduza subdivisoes desnecessarias do espaco perto da assimptota vertical.

3.7.2 Resultados experimentais

A comparacao entre algoritmos em termos de desempenho (tempo e espaco ocupadoem memoria) esta sempre sujeita a diferentes interpretacoes. Por forma a minimizar


∆ Nf t(empo) L (nıveis) N(os)BSP UBSP UQSP BSP UBSP UQSP BSP UBSP UQSP BSP UBSP UQSP

0.01 11.09 25.43 38.50 5.05 6.20 4.84 17 22 11 4139 16163 171330.02 11.65 26.41 41.56 2.64 3.48 2.92 14 20 10 2161 8085 84650.03 12.37 26.32 41.56 1.79 3.44 2.39 13 19 10 1391 5831 84650.04 13.34 27.84 46.83 1.12 1.70 1.21 15 18 9 919 4029 41090.05 13.54 27.65 46.83 0.91 1.72 1.25 12 17 9 715 2905 41090.06 13.49 27.65 46.83 0.86 1.70 1.19 15 17 9 667 2905 41090.07 14.17 28.09 53.34 0.71 0.88 0.64 13 16 8 541 2009 19250.08 14.03 28.09 53.34 0.65 0.88 0.64 13 16 8 527 2009 19250.09 14.33 29.21 53.34 0.53 0.88 0.66 12 15 8 435 1445 19250.1 15.30 29.21 53.34 0.44 0.88 0.67 15 15 8 343 1445 1925

Tabela 3.1: y3 − x2 + 2xy = 0.

a subjectividade da comparacao entre o algoritmo BSP e os restantes, procedeu-se doseguinte modo:

- A comparacao sera feita entre tres algoritmos da mesma famılia: (i) subdivisaobintree nao-uniforme (BSP); (ii) subdivisao bintree uniforme (UBSP); (iii) sub-divisao quadtree uniforme (UQSP) controlada pela aritmetica intervalar.

- A amostragem da curva atraves dos tres algoritmos e levada a cabo com o metodoTGFP.

- A resolucao de amostragem ∆ (distancia maxima admissıvel entre pontos vizin-hos) e a mesma para os tres algoritmos. Mas, no caso da subdivisao quadtree, aresolucao ∆ refere-se a distancia entre os vertices mais afastados dum quadrante.

- O domınio geometrico das curvas esta restrito as curvas algebricas e as curvasracionais definidas por funcoes sem assimptotas porque —tanto quanto o autortem conhecimento— nao sao conhecidas operacoes aritmeticas intervalares paracurvas definidas por funcoes transcendentais.

Os testes de desempenho foram realizados num PC equipado com um processadorIntel Pentium a 800MHz e com 128MB RAM, sendo o Windows XP o sistema opera-tivo. Para estes testes utilizou-se as curvas representadas na Figura 3.14 e considerou-se o espaco do domınio Ω = [−4, 4]× [−4, 4] em R2. Considerou-se ainda dez valoresdiferentes para a resolucao de amostragem, nomeadamente ∆ ∈ 0.01, 0.02, ..., 0.1,como se ilustra nas Tabelas 3.1, 3.2, 3.3 e 3.4. No entanto, as curvas da Figura 3.14estao representadas com resolucao ∆ = 0.1, ate porque, como mostra o grafico (b)da Figura 3.15, o tempo medio de processamento dos tres algoritmos consideradose aproximadamente constante para uma resolucao igual ou superior a ∆ = 0.07. Econveniente precisar que o tempo de processamento t e o tempo total necessario paragerar a arvore e amostrar a curva (passos 1 e 2 do algoritmo ImplicitCurve).


(a) y3 − x2 + 2xy = 0.

(b) y3 − x2 + 2xy − x = 0.

(c) (x− 1y2+1 )(y − 1

x2+1 ) = 0.

(d)0.004+0.11x−0.177y−0.174x2+0.224xy − 0.303y2 − 0.168x3 +

0.327x2y − 0.087xy2 − 0.013y3 +0.235x4−0.667x3y+0.745x2y2−

0.029xy3 + 0.072y4 = 0.

Figura 3.14: Esquemas de subdivisao de curvas racionais sem assimptotas: (esquerda)BSP nao-uniforme, (centro) BSP uniforme, (direita) quadtree uniforme com AI.



0.01 11.84 27.23 42.18 6.58 8.71 6.56 15 20 11 5001 9915 243610.02 12.64 28.52 45.27 2.82 4.22 3.25 15 18 10 2279 4947 120370.03 12.85 28.52 45.27 2.31 4.21 3.25 14 18 10 1761 4947 120370.04 13.11 29.87 48.69 1.45 2.09 1.64 14 16 9 1165 2465 59210.05 13.68 29.87 48.69 1.26 2.09 1.63 12 16 9 999 2465 59210.06 14.63 29.87 48.69 0.98 2.09 1.64 16 16 9 749 2465 59210.07 14.55 32.38 58.30 0.86 1.05 0.88 13 14 8 641 1215 27890.08 14.40 32.38 58.30 0.71 1.05 0.87 11 14 8 569 1215 27890.09 14.20 32.38 58.30 0.70 1.05 0.87 13 14 8 577 1215 27890.1 14.81 32.38 58.30 0.57 1.05 0.87 12 14 8 453 1215 2789

Tabela 3.2: y3 − x2 + 2xy − x = 0.


0.01 8.99 21.54 21.11 4.45 7.44 2.90 17 20 11 4301 8835 101050.02 9.47 22.21 21.68 2.73 3.60 1.43 15 18 10 2273 4401 50370.03 9.43 22.21 21.68 2.03 3.61 1.43 15 18 10 1633 4401 50370.04 10.05 23.52 22.47 1.28 1.76 0.71 16 16 9 1085 2179 24970.05 9.70 23.52 22.47 1.28 1.76 0.73 13 16 9 1039 2179 24970.06 10.09 23.52 22.47 1.07 1.76 0.71 13 16 9 853 2179 24970.07 10.27 24.57 23.44 0.71 0.87 0.36 13 14 8 607 1061 12210.08 10.31 24.57 23.44 0.68 0.87 0.36 12 14 8 563 1061 12210.09 10.41 24.57 23.44 0.70 0.87 0.36 13 14 8 557 1061 12210.1 10.17 24.57 23.44 0.62 0.87 0.36 13 14 8 485 1061 1221

Tabela 3.3: (x− 1y2+1

)(y − 1x2+1

) = 0.

Para todos os valores de ∆, usou-se uma precisao numerica fixa ε = 10−6 nocalculo ou amostragem de cada ponto da curva atraves do metodo TGFP. Note-se quecada ponto amostrado requer que a funcao seja calculada varias vezes. O numero Nf

de vezes que a funcao e calculada por cada ponto da curva depende nao so do metodonumerico usado para aproximar o referido ponto num dos segmentos da grelha departicao, mas tambem do proprio esquema de subdivisao.

Uma analise breve as Tabelas 3.1, 3.2, 3.3 e 3.4 permite concluir o seguinte:

- Numero de vezes que a funcao e calculada versus resolucao. Qualquer que seja oesquema de subdivisao, o numero de vezes que a funcao e calculada aumenta amedida que a resolucao diminui. Isto acontece porque o espaco de procura (cadasegmento XiXi+1) e maior quando a resolucao e mais pequena. No entanto, emmedia, o algoritmo BSP nao-uniforme e claramente melhor que os outros dois,como se mostra no grafico (a) da Figura 3.15.

- Tempo de processamento versus resolucao. Qualquer que seja o esquema desubdivisao, o tempo de processamento t do algoritmo e o numero de nos da



0.01 11.59 24.39 118.55 4.39 6.32 20.39 16 20 7 3395 6667 606850.02 12.90 26.56 139.72 2.13 3.22 10.12 15 18 7 1705 3307 282650.03 14.28 26.56 139.72 1.49 3.18 10.05 13 18 7 1113 3307 282650.04 14.20 29.36 154.14 1.09 1.63 4.82 13 16 7 865 1631 126130.05 15.92 29.36 154.14 0.94 1.64 4.78 12 16 7 721 1631 126130.06 16.53 29.36 154.14 0.75 1.63 4.78 13 16 7 551 1631 126130.07 16.68 36.09 160.46 0.63 0.88 2.28 13 14 7 465 795 56970.08 17.97 36.09 160.46 0.53 0.89 2.36 13 14 7 415 795 56970.09 17.39 36.09 160.46 0.51 0.88 2.31 13 14 7 393 795 56970.1 18.92 36.09 160.46 0.49 0.89 2.31 11 14 7 363 795 5697

Tabela 3.4: 0.004+0.11x−0.177y−0.174x2+0.224xy−0.303y2−0.168x3+0.327x2y−0.087xy2 − 0.013y3 + 0.235x4 − 0.667x3y + 0.745x2y2 − 0.029xy3 + 0.072y4 = 0

arvore diminuem a medida que a resolucao diminui. Repare-se, no entanto,que o tempo de processamento medio do algoritmo BSP nao-uniforme e menorque os outros dois, como se confirma pelo grafico (b) da Figura 3.15. Mas,ocasionalmente, pode acontecer que o tempo de processamento do algoritmo BSPnao-uniforme seja superior, como acontece no caso da curva (c) da Figura 3.14(veja-se Tabela 3.3). Isto deve-se ao facto de, por vezes, um ou mais quadrantesdo primeiro nıvel da arvore nao serem mais subdivididos, como acontece coma subdivisao quadtree da curva referida acima. Mas, tambem, pode acontecerexactamente o oposto. Isto e, o tempo de processamento da subdivisao quadtreeatraves da AI pode ser consideravel, pois para cada no e necessario efectuarvarias operacoes de aritmetica intervalar que requerem sempre o calculo intensivoda funcao. Como mostra a Figura 3.14(d), pode ate acontecer que existamsubdivisoes quadtree onde nao existe qualquer segmento da curva, o que se deveao facto de, em certas zonas do subespaco, o potencial da funcao ser proximo dezero, o que induz o algoritmo a considerar que existe algum segmento da curvanaquela ou naquelas zonas.

- Altura da arvore versus resolucao. Qualquer que seja o esquema de subdivisao,a altura da arvore (i.e. o numero L de nıveis da arvore) e menor a medida quea resolucao diminui. Isto acontece porque sao necessarias menos subdivisoes doespaco. O grafico (c) da Figura 3.15 confirma que, em media, assim e, e quea subdivisao BSP gera mais nıveis na arvore que a subdivisao quadtree, o quenao e surpreendente pois a subdivisao BSP gera dois subespacos ao passo que asubdivisao quadtree gera quatro em cada subdivisao.

- Numero de nos versus resolucao. Qualquer que seja o esquema de subdivisao,o numero N de nos diminui a medida que a resolucao diminui. Isto aconteceporque sao necessarias menos subdivisoes do espaco. Contudo, em media, eapesar da altura da arvore BSP ser maior que a da arvore quadtree, o numerode nos gerado pelo algoritmo BSP e menor que o numero de nos gerado pelo

3.8. SUMARIO 73

0

10

20

30

40

50

60

70

80

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10

∆∆∆∆

Nf BSP

UBSP

UQSP

0

1

2

3

4

5

6

7

8

9

10

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10

∆∆∆∆

t(em

po) BSP

UBSP

UQSP

(a) (b)

0

5

10

15

20

25

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10

∆∆∆∆

L (n

ívei

s) d

a ár

vore

BSP

UBSP

UQSP

0

5000

10000

15000

20000

25000

30000

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10

∆∆∆∆

N(ó

s) BSP

UBSP

UQSP

(c) (d)

Figura 3.15: Comparacao dos tres metodos: BSP, BSU, Quadtree com AI.

algoritmo quadtree, como se confirma pelo grafico (d) da Figura 3.15.

Em resumo, os resultados experimentais obtidos pelos tres algoritmos parecemindicar que o algoritmo BSP, nao-uniforme tem, claramente, melhor desempenho queos outros dois.

3.8 Sumario

Neste capıtulo apresentou-se um algoritmo BSP para a representacao de curvas implı-citas com preservacao topologica. No desenvolvimento do algoritmo BSP, houve queenfrentar e ultrapassar os seguintes problemas:

- Singularidades. Uma curva implıcita em R2 podera ter varias singularidades


(e.g. apices e auto-interseccoes) onde f nao e diferenciavel, o que impossibilitaa utilizacao de metodos numericos baseados em derivadas (e.g. metodo deNewton). Caso contrario, corre-se o risco de o programa quebrar.

- Componentes falsas com variacao de sinal. Todos os algoritmos numericos saobaseados no Corolario de Bolzano. Ou seja, se f variar de sinal ao longo de umsegmento, entao existe um zero de f algures no segmento. Mas, nem sempreisso e verdade. Por exemplo, a curva y − 1/x = 0 tem uma assimptota verticalem x = 0, onde a direita e a esquerda da mesma a funcao tem sinais contrarios.Contudo, a assimptota x = 0 nao e uma componente com variacao de sinal.

- Componentes sem variacao de sinal. Todos os algoritmos numericos conhecidospara a determinacao de zeros sao baseados no Corolario de Bolzano. Tornou-se, pois, necessario encontrar um algoritmo capaz de amostrar componentes comvariacao de sinal, componentes sem variacao de sinal e componentes mistas dumacurva implıcita. Na verdade, poder-se-ia usar algum algoritmo de calculo deextremos para fazer isso, mas estes algoritmos nao sao suficientemente fiaveispara serem combinados com algum algoritmo numerico de calculo de zeros,porque o procedimento de analisar quando uma funcao muda de sinal pode seriludido por uma funcao que tem um extremo (veja-se [40] para mais detalhes).

- Pontos isolados. Tanto quanto se conhece, nao existe nenhum algoritmo paradeterminar pontos isolados de uma curva. Consequentemente, muitos algoritmosnao podem representar curvas degeneradas e nao-homogeneas.

Estes problemas estiveram sempre presentes no desenvolvimento dos algoritmosde representacao de curvas e superfıcies implıcitas, ate porque sofisticados sistemasde calculo matematico (por exemplo, Maple, MathLab e Mathematica) tem mostradodeficiencias na preservacao da forma topologica de grande numero de curvas e su-perfıcies com auto-interseccoes, componentes sem variacao de sinal e pontos isolados.

Capıtulo 4

Superfıcies ImplıcitasNao-Homogeneas

Este capıtulo descreve um algoritmo baseado na subdivisao do espaco para a rep-resentacao de superfıcies implıcitas nao-homogeneas em R3. Para isso, subdivide-seum subespaco Ω ⊂ R3 em celulas cubicas, sendo depois determinada a interseccao detodas estas celulas com a superfıcie.

Ao contrario de outros algoritmos, o algoritmo aqui proposto pode ser usado parapoligonizar nao so superfıcies homogeneas, como tambem superfıcies dimensionalmentenao-homogeneas. Uma superfıcie nao-homogenea contem fragmentos de dimensao 0(pontos isolados), fragmentos de dimensao 1 (segmentos de linha) e fragmentos dedimensao 2. Alem disso, o algoritmo e capaz de resolver auto-interseccoes dumasuperfıcie e, por conseguinte, interseccoes entre superfıcies distintas.

A amostragem dos pontos da superfıcie e levada a cabo atraves do algoritmoBSP (veja-se Capıtulo 3) que determina os pontos da curva que resulta da interseccaode cada face de uma celula com a superfıcie.

4.1 Introducao

Os algoritmos de poligonizacao de objectos implıcitos podem ser classificados emduas grandes classes: algoritmos baseados na subdivisao do espaco e algoritmos devisualizacao directa da superfıcie.

4.1.1 Algoritmos baseados na subdivisao do espaco

Os algoritmos desta classe caracterizam-se por efectuar uma subdivisao espacial deΩ para, numa primeira fase, identificar as celulas que intersectam a superfıcie. Cada

75

76 CAPITULO 4. SUPERFICIES IMPLICITAS NAO-HOMOGENEAS

celula e, depois, classificada com base no sinal do valor da funcao nos vertices dacelula, sendo entao efectuada uma interpolacao linear ao longo das arestas da celulapor forma a calcular pontos da superfıcie.

O tipo de celulas mais utilizadas e o das celulas cubicas, uma vez que possuemvarias vantagens: simetria, facilidade de indexacao, etc. O algoritmo marching cubes[30] e um exemplo dum algoritmo que utiliza celulas cubicas. Mas, por vezes, tambeme usada uma subdecomposicao das celulas cubicas em tetraedros de modo a resolveralguns problemas de ambiguidade. Por exemplo, a deteccao de zeros multiplos dafuncao numa aresta de uma celula cubica, assim como a deteccao de um contornofechado que resulta da interseccao da superfıcie com o interior de uma face, geraproblemas de ambiguidade que podem ser resolvidos —na maior parte dos casos—atraves da decomposicao das celulas em tetraedros.

Os algoritmos baseados na subdivisao do espaco apresentam alguns problemas. Oprimeiro e o tempo dispendido na poligonizacao. Este tempo depende essencialmentede dois factores: o numero de celulas a serem processadas e o processo de calculo dafuncao nos vertices de cada celula.

O segundo problema advem da dificuldade em garantir a consistencia geometricada poligonizacao. De facto, todos os algoritmos conhecidos baseiam-se no princıpio deque so existe algum fragmento da superfıcie onde existe variacao de sinal da funcao nasua vizinhanca (Corolario de Bolzano). Ora, como se viu nos capıtulos anteriores, issonem sempre e verdade. Por exemplo, a superfıcie esferica representada na Figura 4.8(d)tem uma componente sem variacao de sinal, ou seja, todos os seus pontos foramdeterminados sem nunca ter sido satisfeito o Corolario de Bolzano.

O terceiro problema tem que ver com a preservacao da forma topologica dasuperfıcie, ou seja, da sua consistencia topologica. A interseccao entre uma celulacubica e a superfıcie e, como se vera, um subespaco da superfıcie que pode apresentardiversos tipos topologicos ou, se se quiser, subespacos que nao sao topologicamenteequivalentes. Este problema, assim como o anterior, nunca foi devidamente tratadopela comunidade cientıfica. No entanto, algumas tentativas de resolver estes problemasforam feitas por Bloomenthal [38] e Stander [50], ainda que com poucos resultadospraticos.

A subdivisao do espaco podera ser adaptativa ou nao. Em geral, na subdi-visao adaptativa utilizam-se octrees [33]. O controlo da subdivisao pode ser feitoutilizando aritmetica intervalar ou aritmetica afim [51, 19]. No caso da subdivisao nao-adaptativa, as celulas tem tamanho fixo pre-definido, o que nao impede a utilizacaodas aritmeticas anteriores na deteccao da superfıcie.

4.2. SUBDIVISAO NAO-ADAPTATIVA: ESTRUTURA DE DADOS 77

4.1.2 Algoritmos de visualizacao directa da superfıcie

A poligonizacao por amostragem directa da superfıcie baseia-se na tecnica de visu-alizacao por sistemas de partıculas. Quando o sistema de partıculas atinge a suacondicao de equilıbrio, e efectuada uma triangulacao (por exemplo, triangulacao deDelaunay) dessas partıculas sobre a superfıcie.

No entanto, este tipo de algoritmos so garante a consistencia topologica parasuperfıcies homogeneamente 2-dimensionais. Por exemplo, o trabalho de Stander [50]conjuga o uso de um sistema de partıculas com a teoria de Morse para a determinacaode pontos crıticos. A determinacao dos pontos crıticos e feita atraves da matrizHessiana da funcao implıcita. A partir da analise destes pontos crıticos, a consistenciatopologica pode ser mantida a nıvel da forma global da superfıcie, no sentido de queas componentes e os furos da superfıcie sao detectados e representados.

No entanto, a consistencia topologica a nıvel da forma local de superfıcies nao-homogeneas nao esta garantida atraves destes algoritmos. Por exemplo, os pontosisolados so por acaso podem ser encontrados. Alem disso, a utilizacao da Hessianarestringe a utilizacao das tecnicas de partıculas a superfıcies descritas por funcoes difer-enciaveis. Consequentemente, superfıcies com auto-interseccoes e apices nao podemser representadas.

4.2 Subdivisao nao-adaptativa: estrutura de dados

A subdivisao utilizada pelo algoritmo descrito neste capıtulo e uma subdivisao nao-adaptativa. Mas, como se vera adiante, a grande diferenca entre a abordagem aquiapresentada e outras existentes na literatura reside no modo como e tratada cadacelula e na forma como sao resolvidos os problemas referidos na seccao anterior.

Na representacao de superfıcies implıcitas, considerou-se a seguinte estrutura dedados:

class Face

Point *v1, *v2, *v3, *v4;

List<Polyline> *SurfacePoints;

bool flag;

class CubeCell

Point *c;

Real l;

int i,j,k;

Face *ArrayOfFaces[6];


List<Polyline> *SurfacePointsInside;

List<Triangle> *TriangleList;

int CellType;

CubeCell*3DCubeCellGrid[N][N][N];

List<CubeCell> *CubeCellList;

A classe Face caracteriza cada face de uma celula cubica, sendo cada face definidapor quatro vertices v1, v2, v3 e v4. Cada face guarda tambem a informacao dos seuspontos de interseccao com a superfıcie em SurfacePoints.

A classe CubeCell descreve cada celula da decomposicao de Ω. De facto, consi-dera-se que o espaco Ω = [−L,L]× [−L,L]× [−L,L] e subdividido em celulas cubicasdo mesmo tamanho. Cada celula e definida pelo seu centro c e pelo comprimento l dassuas arestas. A variavel ArrayOfFaces[6] permite armazenar as faces da celula. Avariavel SurfacePointsInside permite guardar os pontos amostrados da superfıcieque se encontram no interior da celula. A variavel TriangleList permite guardaros triangulos resultantes da triangulacao da superfıcie. Por fim, a variavel CellTypepermite guardar o tipo de celula.

A estrutura de dados 3DCubeCellGrid, e uma matriz tridimensional que permitearmazenar N3 celulas de tamanho l, tal que N × l = 2L. Na fase inicial do algoritmo,todas as celulas da matriz 3DCubeCellGrid sao inicializadas a NULL. A determinacaoda posicao central (xc, yc, zc) de cada celula no espaco Ω e feita atraves da funcaoΠ : N3 → Ω definida por

xc = 2Li

N− L

yc = 2Lj

N− L

zc = 2Lk

N− L

(4.1)

Ou seja, dada uma posicao (i, j, k) na matriz 3DCubeCellGrid, e possıvel saberqual a posicao (xc, yc, zc) da celula em Ω. Portanto, cada uma das celulas cubicas emΩ pode ser vista como a imagem de um elemento daquela matriz tridimensional. Afuncao Π e importante para gerar a geometria de cada celula a partir do seu centro.

Repare-se que duas celulas adjacentes partilham sempre uma face. E por issoque cada celula tem um array de ponteiros para faces, e nao um array de faces, o quepossibilita a partilha de faces por celulas adjacentes. Isto e particularmente importanteporque evita a repeticao do calculo dos pontos da interseccao da face de cada celulacom a superfıcie.

4.3. O ALGORITMO IMPLICITSURFACE 79

Por fim, a variavel CubeCellList identifica uma lista que permite guardar ascelulas onde existe algum fragmento da superfıcie. Deste modo, como se explicaramais adiante, a visualizacao da superfıcie sera mais rapida do que atraves da matriz3DCubeCellGrid.

4.3 O algoritmo ImplicitSurface

O algoritmo ImplicitSurface baseia-se na subdivisao do espaco Ω. Este espaco eparticionado uniformemente em celulas cubicas, sendo de seguida efectuado o proces-samento de cada uma delas. Assim, tem-se:

Algoritmo 4.1 ImplicitSurface

INPUT:

f : funcao que descreve a superfıcie SΩ : subespaco de R3

Begin

foreach i = 0, . . . , N − 1

foreach j = 0, . . . , N − 1

foreach k = 0, . . . , N − 1

(1) C ← nova celula com centro em (xc, yc, zc) = Π(i, j, k) e lado l;

(2) foreach face fi de C ainda nao processada faz

(a) Φ ← subespaco definido pela face fi;

(b) fi.SurfacePoints ← BSP(f, Φ);

(3) CellPolygonizer(C, f)

(4) if (CellType 6= -1)

(a) 3DCubeCellGrid[i,j,k] ← C

(b) CubeCellList ← CubeCellList + C(5) else

(a) delete C

(b) 3DCubeCellGrid[i,j,k] ← NULL

End


O algoritmo ImplicitSurface atravessa todas as posicoes da matriz de celulas3DCubeCellGrid. Para cada posicao da matriz 3DCubeCellGrid e gerada uma celulaC (passo 1). Na criacao de uma celula (i, j, k) sao normalmente criadas tres novasfaces, pois as outras tres foram ja criadas aquando da criacao das tres celulas adjacentes(i − 1, j, k), (i, j − 1, k) e (i, j, k − 1), como se ilustra na Figura 4.1. As novas facestem a flag a false. No passo 2 determina-se a interseccao das novas faces com asuperfıcie, apos o que a flag de cada uma delas passa a true.

C(i,j,k)

(a)

C(i,j,k-1)

(b) C(i-1,j,k)

(c)

C(i,j-1,k)

(d)

i j

k

Figura 4.1: Criacao da celula C(i, j, k).

Como se sabe, a interseccao entre uma face e a superfıcie e uma curva. Assim,a determinacao desta curva de interseccao pode ser feita atraves do algoritmo BSP(passo 2 do Algoritmo 4.1). Recorde-se que o algoritmo BSP pode ser usado paradeterminar a interseccao de um plano qualquer em R3 com uma superfıcie.

Apos determinar a interseccao entre todas as faces de uma celula e a superfıcie,procede-se a poligonizacao da referida celula, para o que se invoca o subalgoritmoCellPolygonizer do passo 3.

So e guardada uma dada celula na matriz 3DCubeCellGrid se aquela contiveralgum fragmento da superfıcie (passo 4(a) e (b)). No caso em que nao existe qualquerfragmento da superfıcie na celula, procede-se a sua eliminacao (passo 5(a) e (b)).

4.4 Poligonizacao das celulas

Como se viu acima, a poligonizacao de cada celula ocorre no passo 3 do AlgoritmoImplicitSurface, e e descrita pelo seguinte algoritmo:

Algoritmo 4.2 CellPolygonizer

INPUT:

C : celula a tratar

f : funcao que descreve a superfıcie S

4.4. POLIGONIZACAO DAS CELULAS 81

Begin

(1) if (Fr(C) ∩ S 6= ∅)

(a) C.CellType ← DetermineCellType(C)

(b) switch(C.CellType)

1: 1-TypeCellPolygonizer(C, f) return;





(2) else

(a) if (I = 0-TypeCellPolygonizer(C, f))

C.SurfacePointsInside ← C.SurfacePointsInside ∪ IC.CellType ← 0

(b) else C.CellType ← −1

End

Como se mostra no algoritmo anterior, a poligonizacao de uma dada celuladepende do seu tipo. Ha seis tipos de celulas:

Tipo -1. As celulas deste tipo nao contem nenhum fragmento da superfıcie;

Tipo 0. As celulas deste tipo contem um ponto isolado (dimensao 0) no seuinterior (Figura 4.2(a));

Tipo 1. As celulas deste tipo sao atravessadas por um segmento 1-dimensional(dimensao 1), sendo pois caracterizadas pela existencia de dois pontos (extremosdo segmento) na sua fronteira (Figura 4.2(b));

Tipo 2. Este e o caso usual. Uma celula do tipo 2 e transversal a pelo menosum fragmento 2-dimensional da superfıcie que e topologicamente equivalente aum plano (Figura 4.2(c));

Tipo 3. Uma celula deste tipo e uma celula onde a superfıcie se auto-intersecta(Figura 4.2(d));

Tipo 4. Uma celula deste tipo e uma celula onde a superfıcie e topologicamenteequivalente a superfıcie de Cartan (Figura 4.2(e)).


(c) Tipo 2

(d) Tipo 3 (f) Tipo 5 (e) Tipo 4

A

B

(b) Tipo1 (a) Tipo 0

Figura 4.2: Varios tipos de celulas.

Tipo 5. Uma celula deste tipo e transversal a um fragmento 2-dimensional(topologicamente equivalente a um plano) e a um segmento 1-dimensional si-multaneamente (Figura 4.2(f)).

Repare-se que so ha poligonizacao de uma celula quando as suas faces intersectama superfıcie (passo 1 do Algoritmo 4.2). Quando nao existe interseccao das faces dacelula com a superfıcie, entao duas situacoes podem ocorrer (passo 2). Ou existe umponto isolado (passo 2(a)) ou nao existe qualquer fragmento da superfıcie (passo 2(b)),neste caso a celula tem o tipo -1.

P0

(a)

P1

(c)

P0

(b)

Figura 4.3: Processo de determinacao do ponto isolado.

4.4.1 Celulas do tipo 0: pontos isolados

Uma celula do tipo 0 e uma celula que contem um fragmento da superfıcie homeomorfoa R0, i.e. um ponto isolado (Figura 4.2(a)). Um ponto X e um ponto isolado deuma superfıcie se numa pequena vizinhanca de X nao existe qualquer fragmento da


superfıcie a nao ser o proprio ponto X. A determinacao dos pontos isolados e efectuadanas celulas onde nao existe interseccao da fronteira da celula com a superfıcie.

Assim, o processo de poligonizacao de celulas do tipo 0 esta ilustrado na Figura 4.3,o qual e descrito pelo seguinte algoritmo:

Algoritmo 4.3 0-TypeCellPolygonizer

INPUT:

(a) C: celula

(b) fi: face de C

(c) ρ : nıvel maximo de recursividade

(d) ε : erro numerico maximo admissıvel

OUTPUT:

(a) Pi: um ponto isolado

Begin

1. if (ρ == 0) return NULL

2. if (fi == NULL)

Pi ← NULL

nfi ← 0

while (nfi < 6) ∧ (Pi == NULL)

fi ← face numero nfi de C

Pi ← mınimo local no interior da face fi de C

nfi ← nfi + 1

else Pi ← mınimo local no interior da face fi

3. if (Pi == NULL) return NULL

4. if (|f(Pi)| ≤ ε) return Pi

5. fi+1 ← fragmento planar perpendicular a fi em Pi e paralelo a alguma das facesde C dentro de Ω

6. return 0-TypeCellPolygonizer(Ω, fi+1, ρ− 1, ε)

End


A poligonizacao de uma celula do tipo 0 resume-se a determinacao dum pontoisolado. Assim, e de acordo com o Algoritmo 4.3:

• Escolhe-se uma face fi onde exista um ponto Pi que e um mınimo absoluto dafuncao no sentido do interior de C (passo 2 do algoritmo 0-TypeCellPolygonizer)(Figura 4.3(a)). O calculo deste ponto mınimo em fi e feito atraves de umasequencia de mınimos, a semelhanca do algoritmo IsolatedPoint do Capıtulo3.

• Se |f(Pi)| < ε entao encontrou-se o ponto isolado a menos de ε (passo 4 doalgoritmo 0-TypeCellPolygonizer).

• fi+1 ← fragmento planar perpendicular a fi em Pi dentro de C (passo 5 doalgoritmo 0-TypeCellPolygonizer). Esta face tem de ser paralela a algumadas faces de C (Figura 4.3(b) e (c)).

• Invoca-se novamente o metodo para a face fi+1 (passo 6).

4.4.2 Celulas do tipo 1

Uma celula do tipo 1 e uma celula que contem um fragmento da superfıcie home-omorfo a R1, ou seja, um segmento 1-dimensional da superfıcie transversal a celula(Figura 4.4). Este tipo de celula caracteriza-se pelo facto de existirem duas faces comum ponto cada uma (Figura 4.4(a)).

A

B

(a)

A

B

(c)

A

B

(b)

Figura 4.4: Celulas do tipo 1.

Tres situacoes podem acontecer. Na primeira situacao nao existe nenhum seg-mento 1-dimensional entre os pontos A e B (Figura 4.4(a)), o que significa queA e B sao dois pontos isolados. Pode acontecer tambem que nao exista nenhumsegmento entre A e B, mas pelo menos algum deles e um ponto tangente a superfıcie(Figura 4.4(b)). Na terceira situacao existe um segmento 1-dimensional entre os pontosA e B (Figura 4.4(c)).

No caso em que existe um segmento 1-dimensional entre A e B, e necessariocalcular pontos intermedios da superfıcie entre A e B, de forma a permitir efectuaruma poligonizacao correcta da superfıcie. O algoritmo e o seguinte:


A

B

A

BX0

X1

A

B

X0

X1

A

B

X0

X1

Figura 4.5: Determinacao dos pontos de um segmento 1-dimensional de uma superfıcie.

Algoritmo 4.4 1-TypeCellPolygonizer

INPUT:

C : celula

f : funcao que descreve a superfıcie

Begin

(1) l ← ∅;

(2) −→v ← −→BA;

(3) −→u ← −→v10

;

(4) X0 ← B;

(5) repeat

(a) X1 ← X0 +−→u ;

(b) Ωi ← subespaco com centro em X1 e ⊥ a −→u dentro de C;

(c) X1 ← interseccao do subespaco Ωi com a superfıcie, usando o algoritmoBSP;

(d) if (X1 == NULL) return;

(e) l ← l + X1(f) −→u ← X1 −X0

(g) X0 ← X1

until (X1 /∈ Ω);


(6) se (l 6= NULL)

(a) SurfacePointsInside← SurfacePointsInside + Polyline(l)

else

(b) SurfacePointsInside← SurfacePointsInside + Polyline(A)+ Polyline(B)

End

Como se ilustra na Figura 4.5, cada ponto intermedio X1 do segmento 1-dimen-sional entre A e B e determinado pela interseccao do segmento 1-dimensional comum plano que lhe e tendencialmente perpendicular (passo 5(b) e (c)). Este plano eperpendicular ao vector −→u que e aproximadamente tangente a curva em X0. Istosignifica que o plano acompanha as variacoes de curvatura do segmento 1-dimensionalpertencente a superfıcie. Se nao existirem pontos da superfıcie entre B e A (passo 5(d))entao e porque se esta perante os casos ilustrados na Figura 4.4(a) e (b). Se for possıvelir de B a A, entao e construıda uma polilinha com os pontos determinados (passo 6(a)).No caso em que nao existe essa polilinha, cada um dos pontos A e B e guardado numapolilinha singular (passo 6(b)).


Uma celula do tipo 2 e uma celula que contem um fragmento da superfıcie homeomorfoa R2 (Figura 4.2(c)). Recorde-se que antes da poligonizacao de cada celula ja foramcalculados os pontos de interseccao da superfıcie com cada uma das faces de uma celula.Os pontos de interseccao da superfıcie com cada face de uma celula sao determinadosatraves do Algoritmo BSP descrito no capıtulo anterior. Cada ponto e comum asfronteiras de dois subespacos adjacentes, left e right, da subdivisao BSP de cadaface da celula, conforme expresso na classe Point. Deste modo, e ao contrario doutrosalgoritmos existentes na literatura, e trivial formar uma ou mais polilinhas abertas deinterseccao entre a superfıcie e cada face da celula. Noutras palavras, este processo eindependente do numero de polilinhas de interseccao entre a celula e a superfıcie.

Numa celula do tipo 2, o unico problema que subsiste e o de estabelecer a conexaoentre as polilinhas abertas das varias faces de uma celula de modo a formar umaou mais polilinhas fechadas, como se ilustra na Figura 4.2(c). A ligacao entre duaspolilinhas abertas de faces adjacentes faz-se atraves da identificacao do ultimo pontoda primeira com o primeiro ponto da segunda, os quais tem de estar a uma distanciamaxima de ε que e a precisao numerica usada no Algoritmo BSP. Nesta altura, epor forma a evitar lacunas na poligonizacao da superfıcie, elimina-se um dos pontosficando o outro como o ponto comum ou de ligacao de duas polilinhas abertas de facesadjacentes.



Uma celula do tipo 3 e uma celula que contem um fragmento da superfıcie homeomorfoa dois planos cruzados (Figura 4.2(d)). Este caso caracteriza-se pelo facto de haver doispontos A e B de auto-interseccao em duas faces distintas de uma celula. Portanto, haque determinar pontos do segmento 1-dimensional de auto-interseccao da superfıciedentro da celula. Estes pontos sao determinados por um processo semelhante aoilustrado na Figura 4.5. A diferenca e que neste caso a sequencia de pontos adeterminar sao pontos de auto-interseccao que resultam da interseccao de um planocom a superfıcie dentro da celula.

Apos determinar a polilinha de auto-interseccao, ha que sequencia-la com asoutras polilinhas existentes nas faces da celula. A conexao entre duas polilinhas faz-sede forma analoga a das celulas do tipo 2: o ultimo ponto de uma polilinha tem decoincidir com o primeiro da outra polilinha.


Antes da poligonizacao, este tipo de celula caracteriza-se pelo facto de a sua interseccaocom a superfıcie dar origem a uma ou mais polilinhas fechadas e a um ponto isoladonuma das faces da celula (Figura 4.2(e)).

Uma celula do tipo 4 e uma celula que contem um fragmento da superfıciehomeomorfo a superfıcie de Cartan. Portanto, a semelhanca das celulas do tipo 3,as celulas do tipo 4 tambem contem um fragmento 2-dimensional da superfıcie quese auto-intersecta, assim como um segmento 1-dimensional conectado ao fragmento2-dimensional. A determinacao dos pontos de auto-interseccao e efectuado por umprocesso semelhante ao das celulas tipo 3. A determinacao dos pontos do segmento1-dimensional da superfıcie dentro da celula sao determinados por um processo semel-hante ao das celulas tipo 1.

A unica questao que resta resolver e determinar o ponto de uniao Xt entre asduas partes. A determinacao de Xt esta ilustrada na Figura 4.6. Na Figura 4.6(a),o ponto X0 e um ponto do segmento 1-dimensional ao passo que X1 e um ponto deauto-interseccao. O ponto Xt e um ponto que se situa algures entre X0 e X1. Tendoem conta que X0 e X1 estao muito proximos um do outro, toma-se Xt como sendo oponto medio dos pontos P e Q, onde P = AX0 ∩ S e Q = BX1 ∩ S, como se ilustrana Figura 4.6(b).


Uma celula do tipo 5 e semelhante a uma celula do tipo 4, pois contem um segmento1-dimensional conectado ao um fragmento 2-dimensional da superfıcie. A diferenca e


X0

X1

(a)

X0

X1

(b)

Xt

B A

Figura 4.6: Ponto de uniao entre um fragmento 2-dimensional e um segmento 1-dimensional de uma superfıcie numa celula do tipo 4.

que o fragmento 2-dimensional nao se auto-intersecta.

Portanto, uma celula do tipo 5 caracteriza-se por ter somente uma polilinhafechada sem pontos de auto-interseccao e um ponto isolado numa das faces da celula(Figura 4.2(f)). A determinacao do ponto de uniao entre o fragmento 2-dimensionale o segmento 1-dimensional da superfıcie e precisamente igual ao descrito no caso dacelula do tipo 4.

4.5 Triangulacao

Apos a amostragem dos pontos da superfıcie implıcita, torna-se necessario proceder asua poligonizacao. A poligonizacao de uma superfıcie e normalmente levada a cabopor triangulacao. Ha diversos algoritmos de triangulacao de superfıcies homogeneasna literatura [20, 24, 25, 26]. Contudo, a triangulacao de superfıcies nao-homogenease ainda um problema em aberto, ate porque nao sao conhecidos algoritmos parasuperfıcies deste genero.

Assim, no ambito desta tese, e por forma a poder visualizar as superfıcies geradaspelo algoritmo aqui apresentado, houve que desenvolver um algoritmo —ainda quegrosseiro— de representacao grafica de superfıcies nao-homogeneas.

A visualizacao de pontos isolados e segmentos 1-dimensionais da superfıcie re-sume-se a representacao de polilinhas com um ponto (caso dos pontos isolados) ede polilinhas abertas (no caso de segmentos 1-dimensionais), o que e equivalente arepresentar o conteudo de SurfacePointsInside de cada celula.

Na representacao dos fragmentos 2-dimensionais da superfıcie ha que gerar umatriangulacao daqueles nas celulas respectivas, sendo os triangulos armazenados emTriangleList. Esta triangulacao e gerada a partir da polilinha SurfacePoints decada uma das seis faces da celula. Recorde-se que a uniao destas polilinhas das facesda celula formam uma ou mais polilinhas fechadas.

No algoritmo adoptado nesta tese, a triangulacao de cada fragmento 2-dimensio-

nal de uma superfıcie requer o calculo preliminar do centro geometrico PC =Pn−1

k=0 Pk

nde

4.5. TRIANGULACAO 89

cada polilinha fechada que resulta da interseccao entre uma dada celula e a superfıcie.Dado que PC nao pertence normalmente a superfıcie, e preciso calcular a sua projeccaosobre a superfıcie. Desta projeccao resulta um novo ponto amostrado na superfıcie.Este novo ponto resulta da interseccao da superfıcie com uma recta que passa por PC

e que tem a direccao de um vector −→w tendencialmente perpendicular a superfıcie. Oponto PC passa entao a ser o novo ponto amostrado. Uma condicao importante quedeve ser satisfeita e que o novo ponto da superfıcie pertenca ao interior da celula.

A determinacao do vector −→w e feita atraves do produto externo de dois vectores−→u e −→v . Para uma polilinha P0, P1, P2, ..., Pn−1 com n pontos, os vectores −→u =

−−−−→P0Pn/3

e −→v =−−−−→P0P2n/3 sao definidos pelos pontos P0, Pn/3 e P2n/3. Estes pontos formam um

plano que aproxima de alguma maneira a polilinha.

Portanto, o processo de triangulacao de uma polilinha P0, P1, P2, ..., Pn−1 podeser descrito atraves do seguinte algoritmo:

(1) Determinacao do ponto medio PC da polilinha (Figura 4.7(b));

(2) Ajustar o ponto PC a superfıcie;

(3) Determinacao dos pontos Qi = P2i+PC

2(Figura 4.7(c));

(4) Ajustar os novos pontos Qi a superfıcie;

(5) Efectuar a triangulacao (Figura 4.7(d)).

(d) (a)

P0 P1

P2

(b)

PC

(c)

Q0 Q1

P0 P1

P2

Q2

Figura 4.7: Triangulacao de um fragmento 2-dimensional de uma superfıcie.

O ajustamento de PC a superfıcie e essencial para que a triangulacao representeo mais fielmente possıvel a superfıcie (passo 2). No passo 3 calcula-se os pontos Qi quesao os pontos medios de P2i e PC . O passo 4 consiste em ajustar os varios pontos Qi

a superfıcie. O processo de ajustamento dos pontos Qi e de algum modo semelhanteao efectuado aquando do ajustamento do ponto PC . No entanto, o vector director−→w = −→u × −→v de ajustamento de cada ponto Qi e calculado duma forma um pouco

diferente, pois os vectores parcelares sao agora dados por −→u =−−−→P2iPC e −→v =

−−−−−→P2iP2i+1.

Apos o ajustamento dos pontos PC e Qi e efectuada a triangulacao como se mostrana Figura 4.7(d). Este metodo de triangulacao e aplicado somente as celulas do tipo2,3,4,5, pois so estas contem segmentos 2-dimensionais da superfıcie.


4.6 Resultados experimentais

O algoritmo ImplicitSurface foi testado para varias superfıcies. As superfıciesrepresentadas nas Figuras 4.8 e 4.9 sao alguns exemplos de superfıcies algebricas eracionais que foram testadas.

(a) f1(x, y, z) = 44z2 − (1− (x6 )2 − ( y

3.5 )2) (b) f2(x, y, z) = (z − (x2 + y2))((x− 3.9)2 + y2 − 1.22) (x2 + y2 + (z + 2)2) = 0

((x + 3.9)2 + y2 − 1.22) = 0

(c) f3(x, y, z) = x2

7 + y2

8 + z2

9 − 1 = 0 (d) f4(x, y, z) = (x2 + y2 + z2 − 2)2 = 0

Figura 4.8: Superfıcies algebricas.

4.6.1 Apreciacao geral

Como se mostra nas Figuras 4.8 e 4.9, o algoritmo ImplicitSurface preserva a formatopologica das superfıcies implıcitas. De facto, ao inves da generalidade dos algoritmos,o algoritmo permite gerar uma poligonizacao que processa e representa correctamente


os furos, os pontos isolados, as componentes sem variacao de sinal, os segmentos 1-dimensionais remanescentes e as auto-interseccoes duma superfıcie implıcita.

Na Figura 4.8(a) esta representada uma superfıcie algebrica com dois furos. NaFigura 4.8(b) apresenta-se um paraboloide com um ponto isolado. Na Figura 4.8(c)esta representado um elipsoide que nao esta totalmente incluıdo no subespaco ambienteΩ. A Figura 4.8(d) mostra uma superfıcie esferica onde nao existe variacao de sinalda funcao, pois a funcao e sempre positiva em qualquer ponto do espaco, a nao ser napropria superfıcie onde f(x, y, z) = 0.

(a) f5(x, y, z) = (x2 + y2 − z) (b) f6(x, y, z) = x2 − zy2 = 0(x2 + y2 + z − 1) = 0

(c) f7(x, y, z) = x2 + y2

6− 1 = 0 (d) f8(x, y, z) = z − 1

x2+y2 = 0

Figura 4.9: Superfıcies algebricas e racionais.

A superfıcie da Figura 4.9(a) apresenta uma linha de auto-interseccao, pois ea reuniao de dois paraboloides. Na Figura 4.9(b) esta representada a superfıcie deCartan, a qual contem um segmento 1-dimensional (eixo negativo dos ZZ) e umalinha de auto-interseccao (eixo positivo dos ZZ). Note-se que o facto de Ω ser um


subespaco de R3 nao obriga a expressao da funcao a ter tres variaveis. Por exemplo,a superfıcie elıptica da Figura 4.9(c) e descrita por uma funcao em que nao e usada acoordenada z. Como acontece na representacao de curvas atraves do algoritmo BSP,tambem e possıvel representar superfıcies em que o domınio da funcao nao coincidecom Ω, ou seja, pode acontecer que existam pontos de Ω onde o calculo da funcao naoe valido. Por exemplo, na Figura 4.9(d), os pontos (0, 0, z) de Ω nao pertencem aodomınio da funcao.

Apos a realizacao de varios testes de representacao de superfıcies atraves desoftwares comerciais (por exemplo, Maple, MathLab e Mathematica), verificou-se quetais softwares apresentavam falhas notorias na preservacao da forma topologica desuperfıcies implıcitas. Por exemplo, as superfıcies representadas na Figura 4.8(b) e(d) e na Figura 4.9(a) e (b) nao aparecem desenhadas correctamente nestes sistemasmatematicos, quer ao nıvel global (furos e componentes, incluindo pontos isolados),quer ao nıvel local (pontos de auto-interseccao e segmentos 1-dimensionais remanes-centes) da forma topologica.

4.6.2 Desempenho do algoritmo

A Tabela 4.1 apresenta alguns resultados interessantes relativos as superfıcies dasFiguras 4.8 e 4.9 em funcao do comprimento l e do numero de celulas N em queo espaco ambiente Ω foi subdividido. Estas N celulas de Ω estao codificadas ouassociadas a estrutura 3DCubeCellGrid. Por outro lado, K e o numero de celulas ondeexiste algum fragmento da superfıcie, pelo que estao tambem associadas a estruturaCubeCellList.

Os resultados da Tabela 4.1 referem-se ao tempo ∆p de processamento e ao tempototal ∆p + ∆v de processamento e de visualizacao. O tempo de processamento e otempo que leva a processar todas as N celulas desde a amostragem ate ao termino datriangulacao da superfıcie. O tempo de visualizacao ∆v e o tempo dispendido para vi-sualizar a superfıcie triangulada em 3DCubeCellGrid (∆v(Grid)) ou em CubeCellList

(∆v(List)).

Superfıcie l N K ∆p ∆p + ∆v(List) ∆p + ∆v(Grid)f1 0.4 253 506 35.750 36.135 43.321f2 0.8 123 55 7.110 7.390 8.640f3 0.7 103 65 6.350 6.462 7.110f4 0.6 143 80 9.102 9.310 10.780f5 0.5 123 126 8.559 8.660 9.812f6 0.4 123 580 12.327 12.760 13.120f7 0.5 113 484 9.023 9.145 9.687f8 0.8 73 62 5.017 5.141 5.301

Tabela 4.1: Tempos de processamento e de visualizacao de superfıcies (em segundos).

4.7. SUMARIO 93

Pela analise da Tabela 4.1, constata-se que a visualizacao utilizando a estruturaCubeCellList e sempre mais rapida que utilizando a estrutura 3DCubeCellGrid. Istodeve-se ao facto de que em CubeCellList estao somente as celulas onde existe algumfragmento da superfıcie. Na fase da visualizacao, a estrutura 3DCubeCellGrid tambemso aloja celulas que contem algum fragmento da superfıcie. No entanto, testar todasas posicoes de 3DCubeCellGrid vai consumir tempo de processamento. Este facto ebem notorio no caso da funcao f1, pois o numero de celulas que nao contem qualquerfragmento da curva e muito elevado em relacao ao numero de celulas que contemalgum fragmento da superfıcie. Consequentemente, operacoes de natureza interactivacomo as transformacoes geometricas 3D tornam-se bastante lentas e frustrantes.

4.7 Sumario

Os algoritmos existentes na literatura para a representacao de superfıcies impoem umconjunto de restricoes que impossibilitam a representacao de um numero consideravelde superfıcies implıcitas. A maior parte dos algoritmos so serve para representarsuperfıcies dimensionalmente homogeneas e sem auto-interseccoes, embora alguns de-les consigam preservar a forma topologica global (componentes e furos) daquelassuperfıcies.

Ao inves, este capıtulo descreve um algoritmo mais geral de superfıcies implıcitas,pois quer superfıcies homogeneas, quer superfıcies nao-homogeneas podem ser repre-sentadas com preservacao da sua forma topologica local e global. Para conseguiratingir este objectivo fez-se uso do Algoritmo BSP descrito no Capıtulo 3, com o quale possıvel determinar a curva de interseccao entre a superfıcie e cada face das celulasdo espaco Ω.

E e neste ponto que reside a diferenca essencial entre os algoritmos existentesna literatura e este que e aqui proposto. De facto, aqueles algoritmos so determi-nam a interseccao entre a superfıcie e cada uma das arestas das celulas. Daı queestes algoritmos sejam tao limitativos na representacao da forma local e global desuperfıcies implıcitas, pois dificilmente podem determinar as suas singularidades (porexemplo, auto-interseccoes). Por outro, o algoritmo aqui proposto permite inferir assingularidades da superfıcie a partir das singularidades da curva de interseccao entrea superfıcie e cada uma das faces das celulas.

Capıtulo 5

Conclusoes, Contribuicoes eTrabalho Futuro

5.1 Conclusoes

Tanto quanto se julga saber, este e o primeiro trabalho que aborda o problema darepresentacao de curvas e superfıcies dimensionalmente nao-homogeneas. Ao invesde outros trabalhos, os algoritmos propostos nesta tese foram desenvolvidos pararesolver nao so o problema das curvas e superfıcies dimensionalmente degeneradas,mas tambem o problema da preservacao da forma topologica local e global daquelesobjectos geometricos.

Assim, e por forma a atingir estes objectivos principais, no Capıtulo 1 introduziu-se os conceitos fundamentais relativos a funcoes reais de variaveis reais e sua utilizacaona descricao de objectos geometricos, em particular curvas e superfıcies implıcitas.Ainda neste capıtulo, fez-se o ponto da situacao actual relativamente as varias famıliasde algoritmos existentes para a representacao de curvas e superfıcies implıcitas, porforma a contextualizar o trabalho de investigacao que foi desenvolvido ao longo dostres anos do programa de doutoramento.

O Capıtulo 2 surgiu da constatacao das limitacoes dos algoritmos numericosactuais na determinacao de zeros. A principal limitacao destes algoritmos e a condicaode variacao de sinal subjacente ao Corolario de Bolzano, o que impossibilita a deter-minacao de componentes sem variacao de sinal, incluindo pontos isolados, de curvase superfıcies. Daı ter sido proposto um novo metodo numerico, designado por GFP,que tem a vantagem de calcular zeros e extremos de uma funcao atraves de umaunica formula de interpolacao. Foram ainda introduzidas duas variantes deste metodonumerico —metodos TGFP e NGFP— por forma a aumentar a razao de convergenciae desempenho global na procura dos zeros.

No Capıtulo 3 introduziu-se um algoritmo de representacao de curvas implıcitas

95

96 CAPITULO 5. CONCLUSOES, CONTRIBUICOES E TRABALHO FUTURO

nao-homogeneas, em que a amostragem dos pontos e feita em funcao da subdivisaoBSP do espaco ambiente Ω. A determinacao dum ponto da curva resulta sempre dainterseccao da curva com um segmento de bisseccao de Ω, o que e feito com recursoa formula de interpolacao do metodo numerico GFP ou das suas variantes. Alias, aprocura de pontos da curva (ou zeros da funcao descritora) ao longo de cada segmentode bisseccao e exaustiva e feita por recursao ate que alguma condicao de paragem sejasatisfeita. A amostragem dos pontos e realizada independentemente das condicoes decontinuidade e diferenciabilidade da funcao descritora, assim como da forma topologicada curva. Curiosamente, os pontos isolados da curva podem tambem ser determinadosatraves de uma sequencia de mınimos absolutos usando o metodo GFP.

O modo generico como o algoritmo BSP de curvas implıcitas foi desenvolvidopermite determinar a interseccao entre um plano qualquer em R3 e uma dada su-perfıcie. Esta capacidade do algoritmo BSP permitiu a sua utilizacao na amostragemde superfıcies implıcitas, como se viu no Capıtulo 4. Recorde-se que a partir das curvasde interseccao entre a superfıcie e as faces de cada celula torna-se possıvel inferir aforma da superfıcie. Interessante foi notar que as singularidades da superfıcie podiamser inferidas a partir das singularidades das suas curvas de interseccao com as facesdas celulas, o que permitiu tipificar as celulas da decomposicao do espaco ambiente.Isto e de algum modo diferente relativamente aos algoritmos actualmente conhecidos,pois estes o que fazem e determinar a interseccao entre a superfıcie e as arestas dascelulas. Em parte, isto explica os problemas de preservacao da forma topologica dosalgoritmos actuais.

5.2 Contribuicoes

As principais contribuicoes desta tese sao, pois, as seguintes:

• Algoritmo de representacao de curvas implıcitas nao-homogeneas com preserva-cao topologica. Isto significa que componentes 0-dimensionais, componentes 1-dimensionais, assim como pontos de auto-interseccao, sao detectados e amostra-dos, independentemente da satisfacao do Corolario de Bolzano e das pre-condicoesde continuidade e diferenciabilidade.

• Algoritmo de representacao de superfıcies implıcitas nao-homogeneas com preser-vacao topologica. Este algoritmo tem a capacidade de determinar componentes0-dimensionais, componentes 1-dimensionais, componentes 2-dimensionais, as-sim como singularidades, incluindo auto-interseccoes, fragmentos 1-dimensionaise pontos isolados.

No decorrer deste trabalho de investigacao, outras contribuicoes foram apresen-tadas, das quais se destacam:

5.3. TRABALHO FUTURO 97

• Algoritmo(s) numerico(s) para a determinacao de zeros e extremos (maximose mınimos) atraves duma unica formula de interpolacao. Alem disso, este(s)algoritmo(s) numerico(s) nao necessitam de satisfazer qualquer das pre-condicoesusuais dos metodos numericos classicos para funcoes reais.

• Um algoritmo de subdivisao BSP nao-uniforme para a representacao de curvasdefinidas implicitamente por funcoes reais de duas variaveis em R2. Com estealgoritmo, todas as componentes da curva (componentes com variacao de sinal,componentes sem variacao de sinal e componentes mistas) podem ser detectadas,amostradas e representadas. Embora as falsas componentes sejam detectadas,nao sao amostradas nem representadas.

• O proprio algoritmo BSP pode ser visto como um algoritmo de determinacao dospontos de auto-interseccao duma curva implıcita, pois, como ficou demonstradomatematicamente, a subdivisao BSP converge para tais singularidades.

• Algoritmo de determinacao de pontos isolados de curvas e superfıcies implıcitasatraves dum processo de convergencia de mınimos absolutos, sendo cada mınimoabsoluto calculado por algum dos algoritmos numericos do Capıtulo 2.

• Um algoritmo para a representacao de superfıcies definidas implicitamente porfuncoes reais de tres variaveis em R3. Com este algoritmo, todas as componentesda superfıcie (componentes com variacao de sinal, componentes sem variacao desinal e componentes mistas) podem ser detectadas, amostradas e representadas.

• Algoritmo de poligonizacao de superfıcies nao-homogeneas.

• Todos os algoritmos propostos aplicam-se a qualquer funcao real de duas ou tresvariaveis reais, nomeadamente: funcoes algebricas, funcoes racionais e funcoestranscendentais.

5.3 Trabalho futuro

Durante o trabalho realizado, varias questoes foram levantadas que irao mereceratencao no futuro, nomeadamente:

• Determinar matematicamente a razao de convergencia dos metodos numericosTGFP e NGFP.

• Estudar o comportamento do algoritmo BSP em curvas de oscilacao rapida, comopor exemplo, f(x, y) = sin( 1

x). Este estudo ja foi de algum modo empreendido

no artigo intitulado ”A new 2D implicit curve-tracking algorithm”(veja-se [36]),no qual se resolve este problema atraves dum algoritmo de perseguicao.

98 CAPITULO 5. CONCLUSOES, CONTRIBUICOES E TRABALHO FUTURO

• Desenvolver algoritmos de interseccao SSI (Surface-Surface Intersection) combase no algoritmo aqui desenvolvido para superfıcies implıcitas. Isto e particu-larmente importante para determinar a interseccao de solidos em modeladoresgeometricos e sistemas de CAD, ate porque os operadores Booleanos de solidospodem ser facilmente construıdos se se souber a partida qual ou quais sao aslinhas de interseccao entre as suas fronteiras.

• Desenvolver um algoritmo de superfıcies baseado na subdivisao BSP do espacoambiente Ω, e nao na subdivisao regular em celulas cubicas. Deste modo, espera-se obter maior desempenho, pois, em princıpio, o numero de celulas poliedricassera menor do que o numero de celulas cubicas.

• Desenvolver algoritmos de triangulacao e poligonizacao de superfıcies nao-homo-geneas.

• Desenvolver algoritmos interactivos de deformacao de curvas e superfıcies implı-citas.

Bibliografia

[1] S. Abhyankar and C. Bajaj. Automatic parameterization of rational curvesand surfaces I: conics and conicoids. Computer Aided Design, 19(1):11–14,January/February 1987.

[2] S. Abhyankar and C. Bajaj. Automatic parameterization of rational curvesand surfaces II: cubics and cubicoids. Computer Aided Design, 19(9):499–502,November 1987.

[3] S. Abhyankar and C. Bajaj. Automatic parameterization of rational curves andsurfaces III: algebraic plane curves. Computer Aided Geometric Design, 5(4):309–321, November 1988.

[4] E. Allgower and K. George. Numerical Continuation Methods: An Introduction,volume 13 of Springer Series In Computational Mathematics. Springer-Verlag,New York, USA, 1990.

[5] A. Appel. Some techniques for shading machine renderings of solids. In AFIPS1968 Spring Joint Computer Conference, volume 32, pages 37–45, 1968.

[6] D. Arnon. Topologically reliable display of algebraic curves. ACM SIGGRAPHComputer Graphics, 17(3):219–227, July 1983.

[7] A. Barr. Superquadrics and angle-preserving transformations. IEEE ComputerGraphics & Applications, 1(1):11–23, 1981.

[8] P. Baxandall. Vector Calculus. Applied Mathematics and Computing ScienceSeries. Oxford University Press, 1986.

[9] J. Bloomenthal. Polygonization of implicit surfaces. Computer Aided GeometricDesign, 5(4):341–355, November 1988.

[10] J. Bloomenthal, Chandrajit B., J. Blinn, M. Gascuel, A. Rockwood, B. Wyvill,and G. Wyvill. Introduction to Implicit Surfaces. Morgan Publishers Inc, 1997.

[11] R.P. Brent. An algorithm with guaranteed convergence for finding a zero of afunction. The Computer Journal, 14(4):422–425, 1971.

99

100 BIBLIOGRAFIA

[12] P. Bezier and R. Riesenfeld. Mathematical and practical possibilities ofUNISURF. Computer Aided Geometric Design, pages 127–152, 1974.

[13] R. Chandler. A tracking algorithm for implicitly defined curves. IEEE ComputerGraphics & Applications, 8(2):83–89, 1988.

[14] S. Conte and C. de Boor. Elementary Numerical Analysis: An AlgorithmicApproach. McGraw-Hill Book Co., 3rd edition, 1981.

[15] D. Cox, J. Little, and D. O’Shea. Ideals, Varieties, and Algorithms: AnIntroduction to Computational Algebraic Geometry and Commutative Algebra.Undergraduate Texts in Mathematics. Springer-Verlag, New York, USA, 1992.

[16] C. de Boor. A Practical Guide to Splines. Springer-Verlag, 1978.

[17] M. Desbrun, N. Tsingos, and M. Gascuel. Adaptive sampling of implicit surfacesfor interactive modelling and animation. Computer Graphics Forum, 15(5):319–325, 1996.

[18] L. Figueiredo, J. Gomes, D. Terzopoulos, and L. Velho. Physically-based methodsfor polygonization of implicit surfaces. In Proceedings of Graphics Interface(GI’92), pages 250–257, September 1992.

[19] L. Figueiredo and J. Stolfi. Adaptive enumeration of implicit surfaces with affinearithmetic. Computer Graphics Forum, 15(5):287–296, 1996.

[20] S. Fortune. Voronoi Diagrams and Delaunay Triangulations, volume 1 of LectureNotes Series on Computing, pages 193–230. World Scientific, Singapore, 1992.

[21] H. Fuchs, Z. Kedem, and B. Naylor. On visible surface generation by a priori treestructures. Computer Graphics, 14(3):124–133, July 1980.

[22] M. Gascuel and M. Desbrun. Animation of deformable models using implicitsurfaces. IEEE Transactions on Visualization and Computer Graphics, 3(1):39–50, 1997.

[23] J. Hart. Ray tracing implicit surfaces. Siggraph Notes, pages 13.1–13.15, 1993.

[24] E. Hartmann. A marching method for the triangulation of surfaces. The VisualComputer, 14(3):95–108, 1998.

[25] A. Hilton and J. Illingworth. Marching triangles: Delaunay implicit surfacetriangulation. Technical Report CVSSP 01, University of Surrey, Surrey, England,January 1997.

[26] T. Hlavaty and V. Skala. The brute-force generator of triangulations with requiredproperties. In Proceedings of the International Conference on Computer Visionand Graphics (ICCVG’2002), pages 325–336, Zakopane, Poland, September 2002.

BIBLIOGRAFIA 101

[27] J.D. Hobby. Rasterization of nonparametric curves. ACM Transactions onGraphics, 9(3):262–277, July 1990.

[28] C. Hoffmann. A dimensionality paradigm for surface interrogations. ComputerAided Geometric Design, 7(6):517–532, 1990.

[29] H. Lopes, J. Oliveira, and L. Figueiredo. Robust adaptive polygonal approx-imation of implicit curves. In Proceedings of the 14th Brazilian Symposiumon Computer Graphics and Image Processing (SibGrapi’2001). IEEE ComputerSociety, 2001.

[30] W. Lorensen and H. Cline. Marching cubes: a high resolution 3d surfaceconstruction algorithm. Computer Graphics, 21(4):163–169, July 1987.

[31] D. Manocha and J. Canny. Algorithm for implicitizing rational parametricsurfaces. Computer Aided Geometric Design, 9(1):25–51, 1992.

[32] R. Martin, H. Shou, I. Voiculescu, A. Bowyer, and G. Wang. Comparison ofinterval methods for plotting algebraic curves. Computer Aided Geometric Design,19(7):553–587, 2002.

[33] D. Meagher. Geometric modeling using octree encoding. Computer Graphics andImage Processing, 19(2):129–147, 1982.

[34] T. Moller and R. Yagel. Efficient rasterization of implicit functions.OSU-CISRC-11/95-TR49, Ohio State University, 1995.

[35] R. Moore. Interval analysis. Prentice-Hall, 1966.

[36] F. Morgado and A. Gomes. A new 2D implicit curve-tracking algorithm. InProceedings of the International Conference on Computer Vision and Graphics(ICCVG’2002), pages 562–567, Zakopane, Poland, September 2002.

[37] F. Morgado and A. Gomes. A non-uniform binary space partition algorithm for2D implicit curves. In V. Kumar, M. Gavrilova, C. Tan, and P. L’Ecuyer, editors,Proceedings of the International Conference on Computational Science and ItsApplications (ICCSA’2003), volume 3 of Lectures Notes in Computer Science,pages 418–427. Montreal, Canada, May 2003.

[38] P. Ning and J. Bloomenthal. An evaluation of implicit surfaces tilers. IEEEComputer Graphics & Applications, 13(6):33–41, 1993.

[39] A. Opalach and S. Maddock. Disney effect using implicit surfaces. In Proceedingsof the 5th Eurographics Workshop on Animation and Simulation, Oslo, Norway,September 1994.

[40] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling. Numerical Recipes in C:The Art of Scientific Computing. Cambridge University Press, 1992.

102 BIBLIOGRAFIA

[41] A. Ralston and P. Rabinowitz. A First Course in Numerical Analysis. McGraw-Hill Book Co., 1978.

[42] A. Requicha. Representation for rigid solids: theory, methods and systems. ACMComputing Surveys, 12(4):437–464, 1980.

[43] A. Requicha and J. Rossignac. Solid modeling and beyond. IEEE ComputerGraphics & Applications, 12(5):31–44, September 1992.

[44] C. Ridders. A new algorithm for computing a single root of a real continuousfunction. IEEE Transactions on Circuits and Systems, 26:979–980, 1979.

[45] K. Russel. Implicit surfaces for interactive animated characters. Master’s thesis,Massachussets Institute of Technology, The Media Lab, MIT, May 1999.

[46] F. Sederberg, T. Chen. Implicitization using moving curves and surfaces. ACMSIGGRAPH Computer Graphics, 29(3):301–308, August 1995.

[47] C. Shene and J. Johnstone. On the planar intersection of natural quadrics.In Proceedings of the 1st ACM Symposium on Solid Modeling Foundations andCAD/CAM Applications, pages 233–242, Austin, Texas, USA, 1991. ACM Press.

[48] J. Snyder. Interval arithmetic for computer graphics. ACM SIGGRAPHComputer Graphics, 22(2):121–130, July 1992.

[49] B. Stander and J. Hart. Guaranteeing the topology of an implicit surfacepolygonization for interactive modeling. ACM SIGGRAPH Computer Graphics,31(3):279–286, August 1997.

[50] T. Stander. Polygonization implicit surfaces with guaranteed topology. Phd thesis,Washington State University, Washington, USA, 1997.

[51] N. Stolte and A. Kaufman. Parallel spatial enumeration of implicit surfacesusing interval arithmetic for octree generation and its direct visualization. InProccedings of the 3rd International Workshop on Implicit Surfaces (IS’98), pages81–87, Seattle, USA, July 1998.

[52] K. Subramanian and B. Naylor. Converting discrete images to partitioning trees.IEEE Transactions on Visualization and Computer Graphics, 3(3):273–288, July1997.

[53] K. Sung and P. Shirley. Ray tracing with the BSP tree, pages 271–274. AcademicPress Inc., San Diego, California, USA, 1992.

[54] W. Thibault and B. Naylor. Set operations on polyhedra using binary spacepartition trees. ACM SIGGRAPH Computer Graphics, 21(4):153–162, July 1987.

BIBLIOGRAFIA 103

[55] F. Triquet. Fast polygonization of implicit surfaces. In Proceedings of 9thInternational Conference in Central Europe on Computer Graphics, Visualizationand Computer Vision (WSCG’2001), volume 2, pages 283–290, Plzen, CzechRepublic, 2001.

[56] J. Van Aken and M. Novak. Curve-drawing algorithms for raster displays. ACMSIGGRAPH Computer Graphics, 4(2):147–169, April 1985.

[57] L. Velho, L. Figueiredo, and J. Gomes. A unified approach for hierarchicaladaptative tesselation of surfaces. ACM Transactions on Graphics, 18(4):329–360, 1999.

[58] A. Witkin and P. Heckbert. Using particles to samples and control implicitsurfaces. ACM SIGGRAPH Computer Graphics — Proceedings of InternationalConference on Computer Graphics and Interactive Techniques (SIGGRAPH’94),28:269–277, July 1994.

[59] G. Wyvill and A. Trotman. Ray tracing soft objects. CG International ’90:Computer Graphics Around the World, pages 469–476, 1990.

Documents

Universidade da Beira Interior Departamento de Inform ...agomes/PhDthesis-FranciscoMorgado.pdf · ii. Tese realizada sob a ... 2004, Lecture Notes in Computer Science, Vol.3039, Springer-Verlag,