ahd: alternate hierarchical decomposition towards lod ... · d i s s e r t a t i o n ahd: alternate...
Post on 20-Jun-2020
2 Views
Preview:
TRANSCRIPT
AHD: Alternate Hierarchical Decomposition
Towards LoD Based Dimension Independent Geometric
Modeling
Final
January 27, 2011
D I S S E R T A T I O N
AHD: Alternate Hierarchical Decomposition
Towards LoD Based Dimension Independent Geometric Modeling
ausgefuhrt zum Zwecke der Erlangung des akademischen Grades eines Doktors
der technischen Wissenschaften unter der Leitung von
Univ. Prof. Dipl. Ing. Dr. techn. Andrew U. Frank
E127
Institut fur Geoinformation und Kartographie
eingereicht an der Technischen Universitat Wien
Fakultat fur Mathematik und Geoinformation
von
MS. Rizwan Bulbul
Matr. Nr. 0727884
Gymnasiumstrasse 85
1190 Wien
Wien, am 27.01.2011
D I S S E R T A T I O N
AHD: Alternate Hierarchical DecompositionTowards LoD Based Dimension Independent Geometric Modeling
A thesis submitted in partial fulfillment of the requirements for the degree of
Doctor of Technical Sciences
submitted to the Vienna University of Technology
Institute of Geoinformation and Cartography
Faculty of Mathematics and Geoinformation
Advisory Committee:
Univ.Prof.Dr. Andrew U. Frank
Institute for Geoinformation and Cartography
Vienna University of Technology
Univ.Prof.Dr. Walter G. Kropatsch
Institute of Computer Graphics and Algorithms
Pattern Recognition and Image Processing Group
Vienna University of Technology
Submitted by
MS. Rizwan Bulbul
Gymnasiumstrasse 85
1190 Vienna
Vienna, 27.01.2011
3
i
In the name of ALLAH, the most Gracious, the most Merciful
Read! in the name of your Lord, Who created.
Created man from a clot of congealed blood.
Read! and your Lord is most Bountiful.
He Who taught knowledge by the pen.
Taught man that knowledge which he knew not.
Al-Qur’an (Chapter 96, Verses 1-5)
Dedicated to my beloved parents. Dear Aajee and Buba it was possible all
because of your prayers, your efforts spent for educating me, and your
encouragement.
Acknowledgments
First I owe sincere and earnest gratitude to my supervisor Prof. Andrew Frank
for his guidance, advice, and support. He was always available whenever I needed
him, helping, motivating, and encouraging me during the whole duration of my
PhD studies. I also thank Prof. Walter Kropatsch for his insightful comments
and guidance in writing of the dissertation.
I wish to thank my wife, my brothers and sisters for their patience, for mo-
tivating me and being always around me. I feel obliged to all of my colleagues
in the department who supported me at all levels, especially Farid Karimipour,
Gehard Navratil, Gwen Wilke, Ivana Wechselberger, Amin Abdalla, Eva-Maria
Holy, and Christian Gruber.
I cannot forget to thank my friends in Vienna for their support and encourage-
ment especially Aamir Habib, Dr. Basanta Raj Adhikari, Syed Khuram Shahzad,
Guido Petillo, Abbas Chang, M. Hanif, and Shariq Bashir. I will always remem-
ber the time we spent together.
I thank all my friends in Pakistan especially Muhammad Umer, Ejaz Hus-
sain, Dr. Khaliq Aman, Munawar Hussain, Kashif Ali Baig, Shahrzad Khattak,
Shahzad Shigri, Ajmal Hussain, Imtiaz Shah, Abdeen Ali Zia, Yousaf Ali, Ma-
sood Wali, and Meraj Ali. Friends in UK especially Arif Iqbal and Ali Abbas.
Your good wishes and prayers always served as motivation for my achievements
in life.
Finally I wish to thank Higher Education Commission of Pakistan for fund-
ing my PhD studies. This dissertation would not have been possible without
their financial support. I also acknowledge the support and administrative assis-
tance provided by the OeAD office in Vienna (Austrian agency for international
mobility and cooperation in education, science and research) during my stay in
Vienna.
iii
Abstract
The thesis shows that the separation of metric and topological processing for GIS
geometry is possible and opens the doors for better geometric data structures.
The separation leads to the novel combination of homogeneous coordinates with
big integers and convex polytopes. Firstly, the research shows that a consistent
metric processing for geometry of straight lines is possible with homogeneous co-
ordinates stored as arbitrary precision integers (so called big integers). Secondly,
the geometric model called Alternate Hierarchical Decomposition (AHD), is pro-
posed that is based on the convex decomposition of arbitrary (with or without
holes) regions into their convex components. The convex components are stored
in a hierarchical tree data structure, called convex hull tree (CHT), each node of
which contains a convex hull. A region is then composed by alternately subtract-
ing and adding children convex hulls in lower levels from the convex hull at the
current parent node. The solution fulfills following requirements:
� Provides robustness in geometric computations by using arbitrary precision
big integers.
� Supports fast Boolean operations like intersection, union and symmetric
difference etc. Supports level of detail based processing.
� Supports dimension independence, i.e. AHD is extendable to n-dimensions
(n ≥ 2).
The solution is tested with three real datasets having large number of points.
The tests confirm the expected results and show that the performance of AHD
operations is acceptable. The complexity of AHD based Boolean operation is
near optimal with the advantage that all operations consume and produce the
same CHT data structure.
v
Keywords
Convex decomposition, geometric robustness, geometric modeling, geometric rep-
resentation, spatial data modeling, dimension independence, level of detail, Boolean
operations
vii
Contents
Acknowledgements iii
Table of Contents ix
List of Tables xv
List of Figures xvii
1 INTRODUCTION 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Research Question and Approach . . . . . . . . . . . . . . . . . . . 3
1.3 Objective and Hypothesis . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Thesis Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.5 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2 GEOMETRIC DATA REPRESENTATION: SURVEY 7
2.1 Geometric Data Model . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1 Elementary representation schemes . . . . . . . . . . . . . . 9
2.1.2 Hierarchical representation schemes . . . . . . . . . . . . . 10
2.1.2.1 Partition based schemes . . . . . . . . . . . . . . . 10
2.1.2.2 Decomposition based schemes . . . . . . . . . . . 11
2.1.2.3 Constructive solid geometry . . . . . . . . . . . . 12
2.1.2.4 Topological representation schemes . . . . . . . . 13
2.2 The Issue of Robustness in Geometric Computations . . . . . . . . 14
2.2.1 Examples of robustness issues in geometric computations . 15
2.2.1.1 Convex hull computation . . . . . . . . . . . . . . 15
2.2.1.2 Intersection computation of two polyhedra: . . . . 18
2.2.2 Approaches to address the issue of robustness . . . . . . . . 19
ix
x
2.2.2.1 Exact approach . . . . . . . . . . . . . . . . . . . 19
2.2.2.2 Inexact approach or approximate geometry . . . . 22
2.3 An Overview of Geometric Operations . . . . . . . . . . . . . . . . 24
2.4 Survey Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3 MODEL REQUIREMENTS 29
3.1 The Requirement for Improved Robustness . . . . . . . . . . . . . 29
3.1.1 Improved robustness: The requirement . . . . . . . . . . . . 29
3.1.2 Proposed solution: Arbitrary precision arithmetic with ho-
mogeneous coordinates . . . . . . . . . . . . . . . . . . . . . 30
3.1.3 Floating point numbers . . . . . . . . . . . . . . . . . . . . 30
3.1.4 Integers with special arithmetic . . . . . . . . . . . . . . . . 31
3.1.5 Big integers . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.1.5.1 Rational numbers with big integers. . . . . . . . 32
3.1.5.2 Homogeneous coordinates with big integers. . . . 32
3.2 The Requirement for Geometric Primitive . . . . . . . . . . . . . . 33
3.2.1 Geometric Primitive: The Requirement . . . . . . . . . . . 33
3.2.2 Proposed solution: Convex polytope . . . . . . . . . . . . . 33
3.3 The Requirement for Closed Geometric Operations . . . . . . . . . 34
3.3.1 Geometric Operations: The Requirement . . . . . . . . . . 34
3.3.2 Proposed solution: Intersection of convex polytopes . . . . 34
3.4 The Requirement for Level of Detail . . . . . . . . . . . . . . . . . 35
3.4.1 Level of detail: The Requirement . . . . . . . . . . . . . . . 35
3.4.2 Proposed solution: Hierarchical tree data structure . . . . . 35
3.5 The Requirement for Dimension Independence . . . . . . . . . . . 36
3.5.1 Dimension independence: The requirement . . . . . . . . . 36
3.5.2 Proposed solution: Convex polytope . . . . . . . . . . . . . 36
3.6 Requirements Summary . . . . . . . . . . . . . . . . . . . . . . . . 37
4 SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSI-
TION 39
4.1 Homogeneous Coordinates and Big Integers . . . . . . . . . . . . . 39
4.1.1 Experiment for testing the viability of big integers . . . . . 40
4.1.1.1 Ease of programming . . . . . . . . . . . . . . . 40
4.1.1.2 Performance – Approach I . . . . . . . . . . . . 40
4.1.1.3 Size of representation – Approach I . . . . . . . 42
xi
4.1.1.4 Performance – Approach II . . . . . . . . . . . . 43
4.1.1.5 Size of representation – Approach II . . . . . . 43
4.1.2 Discussion on the viability of big integers for GIS . . . . . . 45
4.2 Alternate Hierarchical Decomposition . . . . . . . . . . . . . . . . 46
4.2.1 AHD: An Overview of the Approach . . . . . . . . . . . . . 46
4.2.2 AHD: An Example . . . . . . . . . . . . . . . . . . . . . . . 46
4.2.3 AHD: Convex Hull Tree . . . . . . . . . . . . . . . . . . . . 47
4.2.4 AHD: Functions . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2.4.1 Build function . . . . . . . . . . . . . . . . . . . . 48
4.2.4.2 Build example . . . . . . . . . . . . . . . . . . . . 48
4.2.4.3 Build algorithm . . . . . . . . . . . . . . . . . . . 48
4.2.4.4 Build complexity . . . . . . . . . . . . . . . . . . . 48
4.2.4.5 Eval function . . . . . . . . . . . . . . . . . . . . . 48
4.2.4.6 Eval example . . . . . . . . . . . . . . . . . . . . . 52
4.2.4.7 Eval algorithm . . . . . . . . . . . . . . . . . . . . 52
4.2.4.8 Eval complexity . . . . . . . . . . . . . . . . . . . 54
4.2.5 AHD: Characteristics . . . . . . . . . . . . . . . . . . . . . 54
4.3 Dimension Independent AHD . . . . . . . . . . . . . . . . . . . . . 56
4.3.1 Pseudocode: n-Dimensional AHD . . . . . . . . . . . . . . . 57
4.3.1.1 Convex hull Computation . . . . . . . . . . . . . 58
4.3.1.2 n-Dimensional AHD . . . . . . . . . . . . . . . . 59
4.4 AHD: Level of Detail . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.5 Solution Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.5.1 Robustness . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.5.2 Closed Boolean operations . . . . . . . . . . . . . . . . . . . 62
4.5.3 Dimension independence . . . . . . . . . . . . . . . . . . . . 62
4.5.4 Level of detail . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.5.5 Number of components . . . . . . . . . . . . . . . . . . . . 63
5 AHD BASED GEOMETRIC OPERATIONS 65
5.1 Intersection Operation . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.1.1 AHD based intersection computation . . . . . . . . . . . . . 66
5.1.1.1 Basic intersection operation: Convex-convex in-
tersection . . . . . . . . . . . . . . . . . . . . . . . 67
5.1.1.2 Convex-nonconvex intersection . . . . . . . . . . . 68
5.1.1.3 Generic intersection operation . . . . . . . . . . . 71
xii
5.1.2 Intersection special cases . . . . . . . . . . . . . . . . . . . 74
5.1.3 Intersection solution characteristics . . . . . . . . . . . . . 75
5.2 Complement Computation . . . . . . . . . . . . . . . . . . . . . . . 78
5.3 Union Computation . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.4 Difference and Symmetric Difference Computation . . . . . . . . . 80
5.5 Line-Region Intersection Computation . . . . . . . . . . . . . . . . 81
5.6 Point in Polygon Test . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.6.1 The PIP problem . . . . . . . . . . . . . . . . . . . . . . . . 84
5.6.2 AHD based PIP computation . . . . . . . . . . . . . . . . . 84
5.7 Operations Summary . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6 AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL 87
6.1 Input Representation . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.2 Implementation of AHD . . . . . . . . . . . . . . . . . . . . . . . . 90
6.2.1 2D implementation of AHD . . . . . . . . . . . . . . . . . . 90
6.2.1.1 Geometric data types . . . . . . . . . . . . . . . . 90
6.2.1.2 Geometric data structure . . . . . . . . . . . . . . 90
6.2.1.3 AHD Algebra . . . . . . . . . . . . . . . . . . . . 91
6.2.1.4 Instances for AHD algebra . . . . . . . . . . . . . 93
6.2.2 Dimension independent implementation of AHD . . . . . . 96
6.3 Implementation of AHD based Geometric Operations . . . . . . . . 98
6.3.1 AHD operations algebra . . . . . . . . . . . . . . . . . . . . 98
6.3.2 Implementation of intersection operation . . . . . . . . . . . 99
6.3.3 Implementation of complement operation . . . . . . . . . . 100
6.3.4 Implementation of union operation . . . . . . . . . . . . . . 101
6.3.5 Implementation of difference and symmetric difference op-
erations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.3.6 Implementation of point in polygon . . . . . . . . . . . . . 102
6.4 Implementation of Generalized LoD Data Structure . . . . . . . . 102
6.4.1 Implementation of segment-region intersection . . . . . . . 103
6.4.2 Implementation of strip tree . . . . . . . . . . . . . . . . . 105
7 CONCLUSION AND FUTURE WORK 109
7.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
7.1.1 Is the use of big integers practical? How slow are big integers?110
xiii
7.1.2 Is the combination of convex polytope and big integers at-
tractive? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
7.1.3 Does the model support Boolean operations? . . . . . . . . 111
7.1.4 Support for level of detail . . . . . . . . . . . . . . . . . . . 111
7.1.5 Dimension independence . . . . . . . . . . . . . . . . . . . . 112
7.1.6 Performance analysis . . . . . . . . . . . . . . . . . . . . . . 112
7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.2.1 Support for non simple polygons . . . . . . . . . . . . . . . 117
7.2.2 Implemention of AHD in a database . . . . . . . . . . . . . 117
7.2.3 AHD based geometric modeling of buildings . . . . . . . . . 117
Bibliography 119
Appendix 129
A A Brief Introduction to Haskell 129
A.1 High Order Functions-Functors . . . . . . . . . . . . . . . . . . . . 130
A.2 Haskell Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
A.2.1 Generic data types . . . . . . . . . . . . . . . . . . . . . . . 131
A.2.2 Algebraic data types . . . . . . . . . . . . . . . . . . . . . . 131
A.3 Classes in Haskell . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
B Prototype Application 135
C Haskell Code 137
Biography of the Author 169
List of Tables
2.1 Summary of algorithms for Boolean operations . . . . . . . . . . . 27
4.1 Test program approaches used in experiment . . . . . . . . . . . . 41
4.2 Characterization of the decomposition technique . . . . . . . . . . 55
5.1 Special case treatment . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.1 Functions of Input and Prim classes . . . . . . . . . . . . . . . . . 93
6.2 Operations on simplexes . . . . . . . . . . . . . . . . . . . . . . . . 97
7.1 Description of the three test data sets . . . . . . . . . . . . . . . . 112
7.2 Region statistics where r1, r2, and r3 are input regions: c1, c2
and c3 are the clip regions intersected with regions r1, r2, and r3
respectively: i1=(r1∩ c1), i2=(r2∩ c2) and i3 =(r3∩ c3) . . . . . 115
7.3 Time for AHD and intersection computation . . . . . . . . . . . . 115
B.1 Description of AHD prototype GUI . . . . . . . . . . . . . . . . . . 135
xv
List of Figures
1.1 Separation of metric and topological processing; possible approaches
combined . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1 Properties of a geometric data model . . . . . . . . . . . . . . . . . 9
2.2 A region with the resulting quadtree (from wikipedia) . . . . . . . 10
2.3 CSG tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Geometric computation components . . . . . . . . . . . . . . . . . 15
2.5 Example to show robustness issue for 3D polyhedral intersection
[Hoffmann, 2001] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.6 Fixed grid: some intersecting points are not on the grid points . . 20
2.7 Adding a straight line to a simplicial complex . . . . . . . . . . . . 23
2.8 Line with an envelope . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.1 Two lines intersect, but the computed point is not on either line . . . . 30
3.2 Projection of a homogeneous point to w = 1 plane [Bloomenthal
and Rokne, 1994] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.1 Average time per intersection computation in µsec with related
operations-approach I . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.2 Number of digits per coordinate value – approach I . . . . . . . . 43
4.3 Average time per intersection in µsec with related operations-
approach II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.4 Number of digits per coordinate value – approach II . . . . . . . . 44
4.5 Input polytope and its convex decomposition . . . . . . . . . . . . 47
4.6 CHT for demonstration example of Figure . . . . . . . . . . . . . 47
4.7 AHD steps. L represents the level, dr represents the delta region
and ch represents the convex hull . . . . . . . . . . . . . . . . . . . 49
4.8 Pseudocode of the algorithm build that populates the CHT . . . . 50
xvii
xviii
4.9 Pseudocode of QuickHull algorithm and its helping function qHull,
which computes the sub hulls of partitions recursively . . . . . . . 51
4.10 AHD “eval” using merging of convex hull regions . . . . . . . . . . 53
4.11 Pseudocode of the algorithm eval that processes the CHT to ex-
tract represented region . . . . . . . . . . . . . . . . . . . . . . . . 54
4.12 AHD characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.13 Approaches for dimension independent AHD . . . . . . . . . . . . 56
4.14 n-dimensional AHD- An example for 3D case . . . . . . . . . . . . 57
4.15 Simplexes of dimensions 0 to 3: node, edge, triangle, and tetrahedron 57
4.16 Pseudocode of dimension independent incremental convex hull al-
gorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.17 n-Dimensional AHD . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.18 AHD example for a nonconvex polygon . . . . . . . . . . . . . . . 60
4.19 The concept of LoD in AHD for example polygon in Fig. 4.18. The
increasing level of detail from LoD0 to LoD2 increases the quality
of approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.1 Different polygon intersection scenarios . . . . . . . . . . . . . . . 66
5.2 Intersection computation approach . . . . . . . . . . . . . . . . . . 67
5.3 Convex-convex intersection computation . . . . . . . . . . . . . . . 67
5.4 Pseudocode of convex-convex intersection . . . . . . . . . . . . . . 68
5.5 Convex-convex intersection example . . . . . . . . . . . . . . . . . 69
5.6 Convex-nonconvex intersection computation . . . . . . . . . . . . . 70
5.7 Pseudocode of convex-nonconvex intersection algorithm . . . . . . 71
5.8 Convex-nonconvex intersection example . . . . . . . . . . . . . . . 72
5.9 Generic intersection operation . . . . . . . . . . . . . . . . . . . . . 73
5.10 Pseudocode of generic intersection algorithm . . . . . . . . . . . . 75
5.11 Nonconvex-nonconvex intersection example . . . . . . . . . . . . . 76
5.12 Some special cases for intersection computation . . . . . . . . . . . 77
5.13 Complement R of a given region R . . . . . . . . . . . . . . . . . . 78
5.14 A4B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.15 Difference computation (A\B example) . . . . . . . . . . . . . . . 81
5.16 Line-region intersection . . . . . . . . . . . . . . . . . . . . . . . . 83
5.17 Line-region intersection steps . . . . . . . . . . . . . . . . . . . . . 83
5.18 AHD based point in a polygon test . . . . . . . . . . . . . . . . . . 84
xix
6.1 Hierarchy of geometric data types for region representation . . . . 89
6.2 A Region R represented by two polygons . . . . . . . . . . . . . . . 89
6.3 Relationship between input representation and AHD . . . . . . . . 92
6.4 The function f(x) with three approximations . . . . . . . . . . . . 103
6.5 One piece of a linear function approximating the given function . . 103
7.1 Test datasets and clip regions . . . . . . . . . . . . . . . . . . . . . 113
7.2 Intersection regions . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
7.3 Relationship between time per point (in µs) and max. no. of nodes
for AHD operations build and eval . . . . . . . . . . . . . . . . . . 116
7.4 Relationship between time per point (in µs) on logarithmic scale
and max. no. of nodes for AHD intersection operation . . . . . . . 116
7.5 AHD based geometric modeling of buildings . . . . . . . . . . . . . 118
A.1 Hierarchy of Haskell classes . . . . . . . . . . . . . . . . . . . . . . 132
B.1 Screen-shot of AHD prototype . . . . . . . . . . . . . . . . . . . . 136
Chapter 1
INTRODUCTION
1.1 Motivation
GIS manage the geometry of identified geographic objects together with descrip-
tive data. Research in the field of representation of geometry in GIS was popular
in the 1980s and 1990s, but recently somewhat less attention has been paid.
The current solutions, both in commercial and research software, do not separate
metric and topological processing.
If the two, metric and topological processing, are separated as shown in Fig-
ure 1.1 and an interface between metric and topology established, then we can
investigate independently (a) approaches for metric operations and coordinate
representations and (b) the representation of the topology of geometric objects.
This has been proved advantageous for the construction of libraries for compu-
tational geometry (for example, Computational Geometry Algorithms Library,
CGAL1 (computational geometry algorithms library) and LEDA (Library of Ef-
ficient Data Types and Algorithms)[Mehlhorn and Naher, 1999]) and is likely
useful for GIS.
Figure 1.1 shows 16 possible combinations shown by linking lines between the
approaches for metric and topological processing. The attractive combinations
are;
� the combination of floating point numbers with polytopes; not attractive
because rounding errors will be difficult to control,
� the combination of big numbers with non-topological objects relations;
1http://www.cgal.org/
1
Chapter 1 - INTRODUCTION 2
Figure 1.1: Separation of metric and topological processing; possible approachescombined
topology can be deduced consistently and equality of points checked,
� and big integers combined with polytopes. The combination of big inte-
gers with polytopes (shown with bold link in Figure 1.1) leads to a novel
approach, radically different from the approaches advocated in current re-
search e.g. [Penninga, 2008, Thompson, 2007]. This approach is attractive,
because straight lines remain straight, independent of how many intersec-
tion points are computed (see problem in Figure 3.1) and all results are
consistent.
The current approaches for geometric modeling have limitations; the accuracy of
computational results using current approaches is approximate. The problem of
dealing with approximations is further aggravated by the pressing requirement of
inclusion of 3D (e.g., city models) and temporal data (e.g., movements of people
or cars) causing a re-evaluation of current approaches to geometric modeling.
Revisiting the topic starting with fundamental question is thus warranted:
� The current systems use satisfactory solutions, but the code is complex.
Extension of such solutions, for handling 3D and temporal data, appears
difficult.
� The solutions for geometric computations in CAD (Computer Aided De-
sign), e.g., CGAL, may be helpful but are not necessarily the best solutions
for the GIS domain: For GIS we may restrict geometry to straight lines
and exclude arcs of circle; this is the first distinction between GIS and
CAD, where complex curves are necessary. The second distinction is that
in GIS most points have to be inserted and the new points are very seldom
constructed from previously computed ones.
Chapter 1 - INTRODUCTION 3
1.2 Research Question and Approach
The objective questions of the research are;
� Is the use of big integers practical? How slow are computations with big
numbers?
� Is the combination of convex polytope and big integers attractive?
� Does the model support all operations (like intersection, union and sym-
metric difference etc.)?
The research is planned in following steps;
1. Test the viability of big integers
2. Design the data model
3. Explore data model support for basic Boolean operations
4. Model verification with real data set
5. Test the model for offering additional support for
� Level of detail
� Dimension independence
1.3 Objective and Hypothesis
The objective of this research is to develop a geometric model for the represen-
tation of arbitrary (convex or nonconvex, with or without holes) objects. The
solution that will be based on the separation of metric and topological processing
must fulfill these requirements,
� improved robustness in geometric computations (a requirement for metric
part)
� simple geometric primitive
� efficient geometric operations e.g. intersection, union, symmetric difference
and other (similar) operations
Chapter 1 - INTRODUCTION 4
� generalization to n dimensions (n≥2)
� level of detail oriented processing
Hypothesis: The separation of metric and topological processing is possible if
the issues of consistency in metric processing are addressed.
1.4 Thesis Contribution
In this thesis I have showed that the separation of metric and topological pro-
cessing for GIS geometry is possible and opens the doors for better geometric
data structures. The separation leads to the novel combination of homogeneous
coordinates with big integers and convex polytopes which is shown to be practical
for the representation of GIS geometry that is guaranteed to be consistent.
Firstly I showed that a consistent metric processing for geometry of straight
lines is possible with homogeneous coordinates stored as arbitrary precision inte-
gers (so called big integers). An experiments is performed to test the performance
and storage requirements of big integers for cascaded intersection computation
and to confirm the expected result that big integers (performance wise) are not
expensive for doing GIS geometry.
Secondly I proposed a geometric model, called Alternate Hierarchical Decom-
position (AHD), that is based on the convex decomposition of arbitrary (with or
without holes) regions into their convex components. The convex components are
stored in a hierarchical tree data structure, called convex hull tree (CHT), each
node of which contains a convex hull. A region is then composed by alternately
subtracting and adding children convex hulls in lower levels from the convex hull
at the current parent node.
The modularization in metric processing and topological treatment gives more
freedom to design a data structure which supports the most important operations
in a GIS. The important requirements for GIS today are listed as:
� Support for robustness in geometric computations
� Simple geometric primitive
� Fast processing of intersection, union, symmetric difference and other (sim-
ilar) operations
� Generalization to n dimensions (n≥2).
Chapter 1 - INTRODUCTION 5
� Support for level of detail oriented processing.
The proposed approach, AHD, fulfills theses requirements. Robustness in geo-
metric computations is achieved by using arbitrary precision big integers. AHD
is defined following an algebraic approach and uses convex polytope as a basic
building block (geometric primitive) together with its related operations. This
makes possible to achieve dimension independence (as convex hull computation
is the main operation). The models hierarchical data structure makes the com-
putations efficient especially if implemented in functional languages like Haskell
that support “lazy evaluation”. Operations on arbitrary regions lead to nested
traversal of the trees representing the regions and applying the corresponding op-
eration on the convex building blocks. The algebraic approach suggest to add an
operation complement to the Boolean operation union and intersection - which
is unusual in GIS geometric processing. Implementation of union follows then
using duality from A∪B = (A∩B). AHD has level of detail properties in the sense
that all lower levels in the tree are included in and smaller than the current node.
Tree traversals are thus reduced to the relevant path.
Finally, the solution is tested with three real datasets having large number
of points (in thousands). The tests confirm the expected results and show that
the performance of AHD operations is acceptable. The complexity of AHD based
Boolean operation is near optimal with the advantage that all operations consume
and produce the same CHT data structure.
1.5 Thesis Organization
The thesis is organized as follows;
Chapter 2: This chapter will provide a comprehensive survey of geometric
data modeling schemes. The issue of robustness in geometric computations will
be discussed with some examples and solutions available in literature to deal with
the issue. The last section will discuss the techniques for performing geometric
operations.
Chapter 3: This chapter will list the requirements for the geometric data mod-
eling. This chapter will highlight the basic requirements of the solution for the
metric and the topological processing, the requirement of support for Boolean and
Chapter 1 - INTRODUCTION 6
other related geometric operations, and additional requirements for supporting
dimension independence and level of detail.
Chapter4: This chapter will discuss the solution that fulfills the requirements
listed in previous chapter. First experiment will be performed to show that the use
of big integers for GIS geometry is viable from performance point of view. Next
it will be showm that AHD as a geometric model fulfills the requirements listed
in Chapter 3. AHD based on the convex decomposition of nonconvex polytopes
using homogeneous big integer coordinates for metric processing is needed for
performing robust linear geometry using projective geometry. AHD uses convex
polytope as the basic geometric primitive that provides support for achieving
dimension independence, closed operations and level of detail.
Chapter 5: This chapter will provide the details of using AHD for performing
closed geometric operations. The Boolean intersection is important as the inter-
section of two convex hulls is always convex. In addition, it is the most expensive
operation computationally and other operations like union and difference (which
are not closed over convexity) use it. Intersection and complement operations
will be used to compute union and other Boolean operations using AHD. The
pseudocode of algorithms is provided with appropriate examples that will help in
understanding the algorithmic details.
Chapter 6: This chapter will provide implementation details of all the algo-
rithms presented in Chapter 5. The algorithms will be implemented in Haskell
which is attractive because of its support for lazy evaluation, higher order func-
tions, big integers and compact code.
Chapter 7: In this chapter I will give the conclusion of my work with some
results of AHD for three real datasets. Suggestions for future improvements will
follow.
Chapter 2
GEOMETRIC DATA
REPRESENTATION:
SURVEY
Spatial data models form abstractions from reality which include discretization of
reality which is needed for the reality to be implementable on a computer system
[Schneider, 1997a]. The digital representation of spatial or geometric data 1
describing reality is a key issue [Frank, 1992] and consists of the representation
of geometry of the spatial objects and their associated non-geometric attribute
values. The representation of geometry requires a suitable geometric data model.
A geometric data model provides at an abstract level the geometric primitives
and geometric operations supported by those primitives. The implementation of
a geometric model involves mapping of geometric primitives and operations to
geometric data structures and geometric algorithms respectively. For GIS the
two major models are raster and vector (topological) [Frank, 1992, De Floriani
et al., 1999].
The implementation of a geometric data model is characterized by the follow-
ing properties and shown in Figure 2.1 (motivated by tutorial [Schneider, 2002],
page 11).
1. Closure of geometric operations: The geometric primitives should be
algebraically closed under geometric operations like intersection, union, dif-
ference and so forth. The results can then be used as inputs for further
1The terms spatial and geometric are used as synonyms as for example by Schneider [1997b]
7
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 8
operations. For example, the intersection operation is closed as the inter-
section of two convex hulls is always convex. On the other hand, the union
of two convex hulls is not a closed operation as the union of two convex
hulls may be nonconvex.
2. Efficient data structure: Geometric primitives describe the structure of
the spatial objects and geometric data structures define the storage struc-
ture and organization of primitives facilitating geometric operations to be
performed efficiently and effectively. Geometric data structures together
with geometric algorithms provide implementations details for achieving
the desired behaviour of operations specified in the data model. Geomet-
ric algorithms realizing the spatial operations should be processed as fast
as possible using efficient data structures and methods of computational
geometry [Schneider, 1997a].
3. Robustness of geometric computations: The representational model
should take into account the fixed precision arithmetic available in comput-
ers. Current systems use much improved software for geometric computa-
tions, where all known“holes”were plugged. Will extending this software to
3D or temporal data require the same long period of gradual improvement?
Industry standards like International Standard Organization’s, ISO19107,
is silent on the issue of final digital representation and considering the issue
as an “implementation issue”, that is left to individual programmer to deal
with. The same is true for Open Geospatial Consortium, OGC, specifi-
cations which are based on ISO standards [Alexandroff and Aleksandrov,
1961, Thompson, 2007]. This work reports details about the robustness
issue and survey of approaches to resolve the issue in section 2.2.
4. Dimension Independence: Most of the approaches in literature are di-
mension specific, mostly 2D, 3D (e.g. [Penninga, 2008]) or some are for
specific higher dimensions (e.g. 4D [Breunig, 2001], 5D [van Oosterom and
Stoter, 2010]). Although some dimenion independent methods already exist
e.g. [Gunther, 1988a, Paoluzzi et al., 1993, Feito and Rivero, 1998], these
approaches are not truly dimension independent from an implementation
point of view because the underlying algorithms need to be coded sepa-
rately for each dimension. The task of seprately coding for every dimension
is laborious and needs lot of programming effort if an existing implementa-
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 9
Figure 2.1: Properties of a geometric data model
tion needs to be redone from scratch for adding more dimensions. The need
arises for dimension independent models in which a single implementation
of which works for all dimensions e.g. [Karimipour et al., 2008].
2.1 Geometric Data Model
A geometric data model is used to describe a formalized abstract set of spatial
object classes (geometric primitives e.g. point, line, face, region, volume, simpli-
cial complex, regular polytope, convex polytope etc.) and operations performed
on them. A geometric data structure is the specific implementation of a geomet-
ric data model which fixes the storage structure (e.g. a list, a graph or a tree
etc) and set of algorithms which perform geometric computations [Frank, 1992]
(which have two parts, numerical and topological or combinatorial computation
see Figure 2.4) on the objects stored in the associated structure.
In this section I will give a brief description of related research in the field
of representation of geometry of spatial objects. Gunther [1988a] categorizes the
representation schemes in two, namely elementary representation schemes and
hierarchical representation schemes.
2.1.1 Elementary representation schemes
Elementary schemes do not represent an object by some combination of simpler
objects of the same dimension, but from primitives of lower dimension. Elemen-
tary representation schemes include boundary based representations, sweep rep-
resentation and skeleton representation. Boundary based representation schemes
include techniques for 2D objects, vertex list for general polygons and Fourier
descriptors [Zahn and Roskies, 1972, Persoon and Fu, 1977] for planar curves
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 10
Figure 2.2: A region with the resulting quadtree (from wikipedia)
and for 3D objects including BRep2 and wireform schemes. Sweep representa-
tion schemes are used for the representation of both 2D and 3D objects using
a space curve which acts as the spine or axis of the object. Skeleton schemes
represent a 2D or 3D object by means of a graph. They are not a general pur-
pose representation scheme and useful for giving a rough, short description of an
object.
2.1.2 Hierarchical representation schemes
In these schemes the objects are represented by some combination of simpler
objects of the same dimension. The most common hierarchical schemes include
partition based or decomposition based schemes, constructive solid geometry and
topology based schemes.
2.1.2.1 Partition based schemes
In partition based schemes the object space is divided into non-overlapping parti-
tions. Partition based schemes include region quadtree Samet [1984, 2006] (Figure
2.2), octree3, strip tree [Ballard, 1981] etc. Partition based schemes are not used
as a major representation scheme because they represent an approximation of
the actual object using provided primitives and are sensitive to translation and
rotation operators.
2http://en.wikipedia.org/wiki/Boundary representation3http://en.wikipedia.org/wiki/Octree
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 11
2.1.2.2 Decomposition based schemes
In decomposition based approaches, the object is divided into parts that may
overlap. The spatial data representation model by Ai et al. [2005], uses hierarchi-
cal convex based polygon decomposition approach for multi scale organization of
vector data on spatial data servers. Their approach has a run time4 of O(nlog2n)
and is complex separately dealing holes from the rest of the polygon and no im-
plementation details provided. The work by Keil [2000] focuses on decomposition
of orthogonal polygons. Another hierarchical strategy by Lien and Amato [2006]
decomposes a polygon with or without holes into its approximate convex compo-
nents. Their algorithm results smaller set of approximately convex components
efficiently in time5 O(nr). Chazelle and Palios [1990] have proposed a triangu-
lation based algorithm for the problem of partitioning a polytope in R3. The
algorithm requires O(n + r2) space and runs in O((n + r2)logr) time. A similar
work in the pattern recognition domain has been done by Badawy and Kamel
[2005], presenting a convex hull based algorithm for concavity tree extraction
from pixel grid of 2-D shapes.
Convex decomposition: The problem of decomposing a nonconvex object
into its convex components is known as “convex decomposition” and has appli-
cations in diverse domains ranging from pattern recognition, motion planning
and robotics, computer graphics and image processing, etc. [Chazelle and Palios,
1990, Lien and Amato, 2004, 2006, Xu, 2007, Liu et al., 2008]. The applica-
tion of decomposition techniques in the spatial data modeling domain is seldom
treated in literature, only few address the issue in the spatial database domain
[Kriegel et al., 1991, Ai et al., 2005]. The object representations based on sim-
pler geometric structures are more easily handled [Fernandez et al., 2000] than
complex structures. The algorithms for convex objects are simple and run faster
than for arbitrary objects [Schachter, 1978, Kriegel et al., 1991, Bajaj and Dey,
1992, Palios, 1992] , for example, the decomposition of complex spatial objects
into simple components considerably increases the query processing efficiency
[Kriegel et al., 1991].
Classification criteria for decomposition appraoches: The decomposi-
tion approaches are classified based on the following criteria [Kriegel et al., 1991,
Keil, 2000, Lien and Amato, 2006];
4 n is the number of vertexes5 r is the number of notches or delta regions or concavities
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 12
A. Type of sub-polygons (components)
Depending on application, different decompositions can result in variety of com-
ponents, e.g, convex, star-shaped, or fixed orthogonal shapes.
B. Type of polygon being decomposed
Decompositions can also be classified on the basis of input polygons, e.g., closed
or open, simple or connected, with or without holes, etc.
C. Relation of sub-polygons
Decompositions that result in components, which do not overlap except at bound-
aries, are classified as partitioning. Other decompositions allowing component
overlaps are classified as covering.
D. Objective function
Most of the decompositions are classified based on the objective function. The two
broad categories are decompositions that minimize the number of convex compo-
nents and decompositions that minimize computational complexity in terms of
execution time.
E. Dimension specific or dimension independent
Decomposition approaches can also be classified depending on the input poly-
gon dimension. The majority of the decomposition techniques in literature are
dimension dependent and mostly applicable in 2-D and 3-D [Palios, 1992] cases
only.
2.1.2.3 Constructive solid geometry6
Constructive solid geometry (CSG) together with boundary representation is an
important representation scheme in current CAD and CAM. In constructive solid
geometry a 3D object is represented by 3D volumetric primitives and a set of ge-
ometric operators. Typically the primitives are simple shapes like cuboids, cylin-
ders, prisms, pyramids, spheres, cones and operator are the Boolean set operators
intersection, union and difference. The real advantage of constructive solid ge-
ometry is that basic primitives can be parametrized (thus scaling, translating
6http://en.wikipedia.org/wiki/Constructive solid geometry
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 13
Figure 2.3: CSG tree
and rotating the primitives) and therefore less storage space is needed [Penninga,
2008]. The data structure used is a binary tree, CSG tree (Figure7 2.3), whose
nodes correspond to the operators while the leaf nodes correspond to primitives.
2.1.2.4 Topological representation schemes
The most common approach is the polyhedron approach in which the boundary
of a 3D object is represented by polygonal faces. The polyhedron approach may
have implicit topology or explicit topology. The concept of polyhedral chains
by Gunther [1988a] using convex chains (cells) is a representation scheme for
polyhedra in arbitrary dimensions. Cells are represented as intersection of half
spaces, encoded in a vector.
Winged scheme [Paoluzzi et al., 1993] based on simplicial decompositions is a
representation scheme for piecewise-linear polyhedra of any dimension and curved
polyhedra which are approximated by simplicial maps. The objects are collected
into structures which allowto combine bothother structures and primitive objects
with affine transformations and higher order operators.
An irregular subdivision of space contains nonconvex regions. A region is
defined as a set of polygons [Rigaux et al., 2002], most of them may be nonconvex.
The representation of geometry of convex regions as a collection of inequalities as
y y ≥ ax + b has been proposed in [Gunther, 1988a, Gunther and Wong, 1989].
Penninga [2008] presents a new topological approach for 3D data modeling
based on a tetrahedral network (called TIN or TEN) defined using simplicial
complexes. Hui et al. [2006] proposed a representation method for non-manifold
3D shapes described as 3D simplicial complexes through a decomposition based
approach.
7Source: http://en.wikipedia.org/wiki/Constructive solid geometry
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 14
2.2 The Issue of Robustness in Geometric Computa-
tions
The issue of robustness in geometric computations is a well known and challeng-
ing problem in the implementation of geometric algorithms [Hoffmann, 1989, Li
et al., 2005]. In general, the geometric algorithms are designed based on two
assumptions. First, the inputs are in “general position” meaning that degenerate
input is precluded [Hanniel and Wein, 2007]. Second, geometric algorithms are
designed under a computational model that assumes exact computation with so-
called real RAM model in which the quantities are allowed to be arbitrary real
numbers8. However, these assumptions are not reliable in practice. The input
data may be erroneous and quality of the data may be questionable. Also, in
real implementations of geometric algorithms, exact arithmetic is approximated
by finite precision floating point arithmetic making the implementation notori-
ously difficult [Egenhofer et al., 1990, Burnikel et al., 1999]. The implementa-
tions based on approximations work well for majority of the geometric problem
instances, but occasionally fail as the rounding errors due to fixed precision can
cause catastrophic errors [Schirra, 1998] including program crashes, infinite loops,
wrong and inconsistent results.
The issue of robustness in geometric algorithms arises from the special nature
of geometric computation as it not only consists numerical part but also the
combinatorial part [Li et al., 2005]. The numerical part is further divided into
two parts termed predicate and constructor [Shewchuk, 2009]. Numerical part,
is therefore, used in the computation of geometric predicates (e.g. is a given
point on a line?) and construction of new geometric objects (e.g. computation
of the intersection point of two intersecting lines). The numerical computation,
especially the predicate component is used for determining the combinatorial
relations among geometric objects, shown by a bold arrow in the Figure 2.4. The
fixed precision arithmetic provided by machines can produce numerical errors
leading to inconsistencies, which in turn may lead to inconsistent combinatorial
structures or inconsistent program states [Li et al., 2005].
The importance of the robustness in geometric algorithms is evident by the
wide range of their application domains ranging from computer aided design
(CAD), computer aided manufacturing (CAM), computer aided engineering (CAE),
8http://www.cs.berkeley.edu/˜jrs/meshpapers/robnotes.pdf
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 15
Figure 2.4: Geometric computation components
image processing, computer graphics, robotics and GIS. In all these domains, the
geometric operations based on fixed precision arithmetic may lead to inconsistent
results or complete algorithmic crashes because of rounding errors in computa-
tions. Early GIS research focused—among other things—on overlay calculations
(e.g., [Chrisman et al., 1992]). In the late 1980s the results on testing overlay
operations in the then available commercial GIS were disastrous. The result
was disastrous: none could perform the (intentionally difficult) problem with-
out errors; many of the program’s operations stopped unexpectedly and others
produced substantially wrong answers.
The issue of robustness makes the process of algorithm development costly and
hinders full automation [Smith, 2009] and may lead to inconsistent combinatorial
structures or inconsistent program states [Li et al., 2005].
First I will discuss some geometric algorithms from literature to illustrate the
issue of robustness in geometric computations and then discuss the methods used
to resolve the issue.
2.2.1 Examples of robustness issues in geometric computations
2.2.1.1 Convex hull computation
Kettner et al. [2008] have shown an example of the failure of an algorithm for
the computation of incremental convex hull for a finite set of points in 2D. The
incremental convex algorithm starts by taking any three non-collinear points and
then adding points one by one. Whenever a point is added, it is checked whether
to be inside or outside the previously computed convex hull. If it is inside the
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 16
convex hull, nothing changes, otherwise, the convex hull is recomputed with the
new point to get an updated hull.
The mathematical predicates needed for implementing the algorithm are;
Orientation Predicate: Given three points p = (px, py), q = (qx, qy), and
r = (rx, ry) the predicate tells whether the points are oriented clockwise, coun-
terclockwise or collinear. The arithmetic value that determines the value of the
predicate is the signed area function defined as;
twicearea(p, q, r) = (qx=px)(ry=py)=(qy=py)(rx=px) (2.1)
Equation 2.1 computes the twice the area of the triangle made by three points
p, q and r. The sign of the area in equation 2.1 provides the orientation state i.e.
orientation(p, q, r) = sign(twicearea(p, q, r)) (2.2)
> 0 counterclockwise
< 0 clockwise
= 0 collinear
Positive sign in equation 2.2 denotes the counterclockwise orientation, while
negative denotes the clockwise and zero value shows that the points are collinear.
When the orientation predicate is implemented using floating point arithmetic,
the three potential ways in which the result could differ from the correct orien-
tation are mentioned by Kettner et al. [2008]:
� rounding to zero: ’+’ or ’=’ signs are miss-classified as 0;
� perturbed zero: 0 is miss-classified as either ’+’ or ’=’;
� sign inversion: ’+’ is miss-classified as ’=’ or vice-verse.
Visibility Predicate: Given an edge e = (pi, pi+1) and a point q, the predicate
determines whether the point q is visible to edge e or not. The visibility predicate
is deduced from the orientation predicate.
visible((pi, pi+1), q) = orientation(pi, pi+1, q) (2.3)
Negative orientation i.e. (twicearea(pi, pi+1, q)) < 0 means that the point q is
visible to edge e = (pi, pi+1), and weakly visible if (twicearea(pi, pi+1, q))≤ 0.
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 17
Kettner et al. [2008] have identified two key properties for the operation of
the algorithm;
Property A: A point p is outside a convex hull, ch, iff p is visible to any
edge of the ch.
Property B: For a point outside ch, the weakly visible edges of ch make a
nonempty subset of consecutive edges of ch. Likewise the edges not weakly visible
make a nonempty subset of consecutive edges of ch.
The algorithm is based on these two properties. Suppose that (vi, vi+1), ..., (v j=1, v j)
is the sub-sequence of weakly visible edges of ch, then the updated hull is obtained
by replacing the sub-sequence (vi+1, ...,v j=1) by p.
Thus the key algorithmic steps are;
1. Determination of any visible edge from p: This is done iteratively
using the visibility predicate over the list of convex hull edges L. If no visible
edge, then p is neglected. Otherwise, continue with a any visible edge e to
next step.
2. Determination of weakly visible edges from p: This is also done by
using visibility predicate. Starting with any visible edge e, iterate over L in
counterclockwise direction until a non weakly visible edge is encountered.
The step gives the subset sequence (vi, vi+1), ..., (v j=1, v j) of L.
3. Updating of Convex hull: The edges of subset sequence determined in
previous step (vi, vi+1), ..., (v j=1, v j) are removed from L and two new edges
(vi, p) and (p, v j) are added to L.
The properties A and B may not, however, may not hold for computations done
with floating point arithmetic.Kettner et al. [2008] have identified four logical
ways to negate the properties A and B:
Failure 1: There is not any visible edge for a point outside the current hull.
Failure 2: There is a visible edge for a point inside the current hull.
Failure 3: All the edges of the convex hull are visible for a point outside the
current hull.
Failure 4: There is a non-contiguous set of visible edges for a point out the
current hull.
Kettner et al. [2008] have provided concrete examples to demonstrate these
failures. The effects of these failures could be as shown by Kettner et al. [2008]
are;
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 18
(a) Intersection of a cube witha tetrahedron
(b) Possible inconsistent facepartitions
Figure 2.5: Example to show robustness issue for 3D polyhedral intersection[Hoffmann, 2001]
� The algorithm computes a convex polygon, but misses some of the extreme
points.
� The algorithm crashes or does not terminate.
� The algorithm computes a non-convex polygon
2.2.1.2 Intersection computation of two polyhedra:
Another example of failure of geometric computation drawn from polyhedral
intersection is given by Hoffmann [2001] as shown in Figure 2.5a. The intersection
is done by computing the face-face intersection points i.e. by computing the
intersection points of edges of one solid with the faces of the other solid.
Assuming that the the front edge of the tetrahedron, eFT , cuts the top face of
the cube, fTC, at a steep angle and the front face of the cube, fFC, at a shallow
angle, it is possible that the intersection eFT∩ fTC deemed to be in the interior
of the face as shown by A in Figure 2.5b while the intersection eFT∩ fFC may be
deemed to lie on eFT as shown by B in Figure 2.5b. The two partitions resulted
from incidence decisions are logically inconsistent , because the front face of the
tetrahedron, fFT , now intersects twice, one with fTC and one with fFC at A and
B respectively as shown in Figure 2.5b. This situation leads to the failure of the
intersection algorithm.
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 19
2.2.2 Approaches to address the issue of robustness
The state of the art is best reflected by libraries like CGAL and LEDA[Mehlhorn
and Naher, 1999]. It has been pointed out that one way out of the rounding
problem is to compute each value only once [Knuth, 1992]; contradiction is then
logically impossible (also described as “Do not ask stupid questions” [Fortune,
1989]). This had been the guiding idea for Frank and Kuhn [1986] earlier. Hoff-
mann [2001] identifies three categories of strategies to address the robustness
problem;
1. Using exact arithmetic (integer arithmetic, extended precision arith-
metic, nonstandard or symbolic arithmetic).
2. Using symbolic reasoning (reuse previously computed results)
3. Using reliable calculations (interval arithmetic)
These approaches can broadly be divided into two categories based on either
using inexact or exact arithmetic[Li et al., 2005, Smith, 2009];
1. Inexact approaches: Approaches based on inexact arithmetic - making
fixed precision robust.
2. Exact approaches: Approaches based on exact arithmetic - making exact
approach viable.
2.2.2.1 Exact approach
This approach uses the exact rational arithmetic to avoid the problems associated
with fixed precision floating point arithmetic. The computation is done with a
precision that is sufficient to keep the theoretical correctness of an algorithm,
thus allowing exact implementation of algorithms developed for real arithmetic
without modifications. However, there is criticism for using an exact approach
as it is thought to be slow computationally, especially for cascaded operations
where the output of one operation is used as input of another operation [Schirra,
1998, Li et al., 2005].
The two problems that make an exact approach complicated are [Hoffmann,
2001];
1. Proliferation: If the input representation of a geometric computation has
a precision of k digits, the output may require a multiple of that precision to
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 20
Figure 2.6: Fixed grid: some intersecting points are not on the grid points
be representable e.g. consider the intersection computation of line segments
(whose point coordinates are represented by fixed precision) in a fixed grid
as shown in Figure 2.6. It is seen that, while some of the intersection points
are representable being on the grid points, others are not. That is, the
input points have a precision of k digits, the output computed points need
a multiple of input precision to be representable. For cascaded operations,
this may lead to an exponential growth in the required precision.
2. Irrationality: Some computations may result in output which is not rep-
resentable in the given fixed precision e.g. consider rotating a unit square,
whose vertex coordinates are represented by fixed precision, by 45◦. The
new vertex coordinates of the rotated square now involve radicals and are
irrational. Franklin [1984] has also shown that scaling or rotating an object
can cause a point to move to the other side of a line, thus failing point
inclusion tests.
Many geometric algorithms may involve computations involving integers and or
rationals. Every rational value can be represented by a pair of integer values
indicating the numerator and denominator. Infinite precision rational numbers
are also suggested by [Franklin, 1984] has the advantage of satisfying the mathe-
matical field axioms [Thompson, 2007]. The issue with overflow in such computa-
tions can be avoided using arbitrary precision arithmetic provided by big number
software packages, arbitrary precision integer (for integers aka big integers) and
arbitrary precis on rational (aka big rational represented by a pair of big integers).
Many software libraries and programming languages provide support for big
number packages e.g. LEDA [Mehlhorn and Naher, 1999]. The functional pro-
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 21
gramming language Haskell [Jones and Jones, 2003] also provides support for
arbitrary precision integers, also known as big integers. However, the use of
big integers is not encouraged in the literature accounting for slowing down the
computation especially the cascaded computations. Karasick et al. [1991] have
reported that the implementation of 2D Delaunay triangulation using arbitrary
precision arithmetic is 10,000 times slower than floating point implementation.
Computing exactly may not always be needed. Thus, the computational effi-
ciency of geometric algorithms involving exact computation using arbitrary preci-
sion arithmetic can be improved by technique of adaptive evaluation (also known
as lazy evaluation). Adaptive evaluation computes exactly only if it is necessary
to do so, otherwise it uses the speed of approximate arithmetic (floating point
arithmetic). The simplest form of adaptive evaluation is to use filters e.g. floating
point filter, which filters computations where floating point arithmetic gives the
correct result. A filter compares the absolute value of the computed numerical
value of the expression to the error bound. If the error bound is smaller, the
computed approximation and the exact value have similar sign. Otherwise, the
expression is reevaluated using exact arithmetic. Depending on the computation
of error bounds, filters are classified as static or dynamic. In case of so called
static filters [Fortune and Van Wyk, 1996], error bounds are computed a priori
if the specific information on the input data is available. Static filter is efficient
because it makes use of information known at compile time [Smith, 2009]. On
the other hand, dynamic filter computes error bounds on the fly using the run
time information and so is less efficient.
Franklin [1984] discusses the idea of using interval arithmetic in which real
numbers are represented by intervals, whose endpoints are floating point numbers.
The interval representing the exact value of an operation is computed by floating
point operations on the endpoints of the interval representing the operands. The
lower and upper bounds of a value are thus rounded to nearest or the rounding
mode is changed to rounding toward −∞ and +∞ respectively [Schirra, 1998].
Fortune [1995] has showed that boundary based polyhedral modeling support-
ing regularized set operations can be implemented using exact arithmetic with
minimal performance overhead compared to floating point arithmetic.
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 22
2.2.2.2 Inexact approach or approximate geometry
As discussed, the algorithms based on fixed precision arithmetic may fail oc-
casionally because of the rounding off errors. Exact computation is thought as
being computationally expensive, so inexact methods which can tackle the robust-
ness issue are used. Avoiding inconsistencies among decisions is a primary goal
in achieving robustness in implementations with imprecise arithmetic [Schirra,
1998]. With approximate geometry, the method of representation and the con-
straints that apply are selected to maintain certain properties and ensure correct
behavior of the algorithm [Smith, 2009]. Some design principles for robustness
under computation with imprecision as listed by Schirra [1998] follow.
Representation and model approach introduced by Hoffmann et al. [1988] for-
malizes the “compute the correct solution for an input” idea [Schirra, 1998]. The
model distinguishes between the models (real mathematical objects), and their
computer representations. An algorithm is correct if the model corresponding to
the computed output representation is the solution to the problem for the model
corresponding to the input representation. For proving the robustness of an algo-
rithm, it is needed to show that there is always a model for which the algorithm
takes correct decisions. But this is a highly non-trivial task [Schirra, 1998].
Epsilon geometry introduced by Salesin et al. [1989] is a theoretical framework
in which the predicate returns a real number instead of a Boolean. The technique
computes an exact solution for a perturbed version of the input and returns a
bound on the size of this implicit perturbation. Epsilon geometry has been applied
to few basic geometric primitives and the framework seems difficult [Schirra,
1998].
Topology oriented approach focuses on the topological and combinatorial cor-
rectness than numerical correctness (e.g. simulation of simplicity [Edelsbrunner
and Mucke, 1990], Realms [Guting and Schneider, 1993], ROSE algebra [Guting
and Schneider, 1995]). The decisions based on numerical computations that may
lead to topological inconsistencies or violations are not taken and are replaced
by topology conforming decisions. This is achieved by providing set of rules,
following these lead to topological consistency [Schirra, 1998]. In a most recent
work, Smith [2009] has provided two algorithms for performing Boolean opera-
tions on piecewise linear shapes. The algorithms are guaranteed to generate a
topologically valid result from topologically valid input.
Simplex based approaches fall into the same class of topology based appraoch
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 23
Figure 2.7: Adding a straight line to a simplicial complex
and they define the world with points, line segments and tetrahedrons. Insert-
ing all geometry into a simplicial complex was proposed by Frank and Kuhn
[1986] where presentations stress formal definitions building multi-sorted alge-
bras. This means that every point—including any intersection point is inserted
in a triangulation [Egenhofer et al., 1990]; the arithmetic precision of this process
is not important (it was described as “oracle”). Important is that all topologi-
cal relations are deduced from the topology of the triangulation—no repetition
of finite precision computation is permitted, to avoid inconsistencies by “asking
stupid questions” i.e. question to which the answers are already known. The
disadvantages of this solution are; (1) processes to insert new points or lines
establish topology with relation to all nearby objects, which may never be of in-
terest—resulting in large amounts of unnecessary computations and (2) straight
lines are broken in small pieces when they are inserted into the simplicial complex
and their straightness may be lost due to rounding (Figure 2.7).
The construction of the simplicial complex produces a large number of smaller
objects [Egenhofer et al., 1990]. A commercial implementation (TIGRIS Inter-
graph, [Herring, 1991]) avoided the division to triangles and used a generalization
of simplicial complexes, namely cell complexes, but lost also the advantage of the
fully defined topology. The approach relates to quad edge algebra for partitions
[Guibas et al., 1983].
Numerous investigations in computing geometry with integers, dubbed “ge-
ometry on a grid”. Such approaches often consider the envelope of a line as the
set of grid points and the intersection points of this line with other lines may
belong (see Figure 2.8) to this set of grid points [Greene and Yao, 1986]. The
method of [Guting and Schneider, 1993] has been built on a special geometric
engine to compute geometry in which certain types of errors are avoided.
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 24
Figure 2.8: Line with an envelope
2.3 An Overview of Geometric Operations
Boolean operations on geometric data have been studied in diverse domains
ranging from image processing, computer graphics [Weiler and Atherton, 1977,
Thibault and Naylor, 1987, Greiner and Hormann, 1998], robotics, CAD and
CAM [Lauther, 1981, Paoluzzi et al., 1993]. Here I provide a brief survey of al-
gorithms for performing intersection, union and symmetric difference operations
and Table 2.1 provides a summary of the surveyed literature for Boolean opera-
tions. In GIS Boolean set operations are essential for the manipulation of spatial
data needed for performing spatial queries e.g for overlay operations, clipping
etc. Traditional algorithms may fail because of the robustness issues in geometric
computations and (or) may need special preprocessing of degenerate cases.
Many algorithms for Boolean operations on polygons have been reported in
the literature. The preliminary work by Shamos and Hoey [1976] provided a basis
for geometric intersection problems. They have shown that the intersection of two
simple plane n-gons can be detected in O(n logn). The intersection of two convex
n-gons and two nonconvex n-gons can be computed in O(n) and O(n2) respec-
tively. The Wieler-Atherton algorithm [Weiler and Atherton, 1977] for intersec-
tion of two simple polygons with vertexes n and m of any shape with or without
holes has the worst time complexity of O(nm). Bentley and Ottmann [1979] gave
the classical sweep line algorithm for counting and reporting all intersections by
extending the work by Shamos and Hoey. They provided an algorithm for report-
ing all k intersections between two general planar objects in O(n logn + k logn).
The work by Bentley and Ottmann was further extended by Lauther [1981], and
the reported sweep line algorithm for intersection computation has the expected
time complexity of O(n logn). Another O(n) time algorithm was presented by
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 25
O’Rourke etal [1998]. The algorithm is simple but is limited to convex polygons
only. The two plane sweep algorithms by Nievergelt and Preparata [1982] com-
pute the geometric intersection of two nonconvex polygons in O((n+k) logn) and
two convex polygons in O(n logn + k). The polygons can have self intersecting
edges but degenerate cases are not tackled. The data structure is complex and
implementation details are not given.
The work by Chazelle and Dobkin [1987] provides lower bounds on algorithms
for intersection of convex objects in two and three dimensions. Their work is based
on the assumption that the intersecting objects are available in random access
memory, eliminating reliance on linear input reading time. The time bounds for
two convex polygons in 2D and two convex polyhedra in 3D cases are O(logn)
and O(log3 n) respectively (in 3D case an additional multiplicative factor of logn
for data structure preprocessing for standardization). The work by Margalit and
Knott [1989] presents an algorithm for set operations on polygon pairs having
worst time complexity of O(n2). They give partial correctness proof of their
solution and implementation is discussed but is still complex and not easily un-
derstandable. The approach by Rappoport [1991] for the representation of n-
dimensional polyhedra called extended convex difference tree (ECDT) extends
the convex difference tree approach. The approach can deal only two polygons
and can not deal with holes and islands. It also needs special topological process-
ing to avoid the robustness issues.
The state of the art for finding all intersections among segments is given by
Chazelle and Edelsbrunner [1992]. The algorithm by Vatti [1992] is for clipping
arbitrary polygons against arbitrary polygons. The polygons may be convex,
concave or self intersecting. However, self intersecting polygon is converted to a
non-intersecting polygon by inserting the points of intersection during the clip-
ping process. The algorithm also supports polygon decomposition by allowing
the output in the form of trapezoids if required. The solution is complex and
implementation is not easy. Its performance has not been proved asymptoti-
cally, rather a comparison with traditional clipping methods is provided. The
algorithm proposed by Greiner and Hormann [1998] also deals clipping arbitrary
polygons like Vatti’s algorithm. However, it is a simple algorithm based on the
boundary segment manipulation and performs better than Vatti’s algorithm over
randomly generated general polygons. The data structure is a doubly link list or
lists in case of multiple polygons. Only few degenerate cases are mentioned and
it can not treat overlap degeneracies. The solution is limited for 2D polygons
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 26
and robustness issues related to the fixed precision floating point arithmetic are
not considered. The algorithm by Rivero and Feito [2000] to calculate Boolean
operations for general planar polygons (manifold and non manifold, with and
without holes) is based on simplicial chains and their operations. The strategy
has been demonstrated for 2D and claimed to be valid for 3D polyhedra. The
algorithm does not need special treatment of degenerate cases, and has same time
complexity as the Greiner’s algorithm. The slight modifications in the work of
Rivero and Feito by Peng et al. [2005] resulted in an algorithm which has been
shown to be more efficient (execution time less than one third of that by Rivero
and Feito). CGAL [Fogel et al., 2006] provides Boolean operations for polytopes
in 2-dimensional Euclidean space. Robustness is ensured through the use of exact
arithmetic. The regularized operations are provided for two simple polygons with
or without holes. The time complexity is O(n2) for simple polygons. The algo-
rithm by Kui Liu et al. [2007] is for clipping arbitrary polygons with holes. The
algorithm is based on segment manipulation and works by classification of inter-
section points into entry or exit points. Unlike the solutions by Vatti and Greiner,
this algorithm uses a single linked list data structure and performs better than
Vatti’s solution for smaller number of input points. The solution is limited to 2D
polygons and modifications are needed for dealing multiple polygons and holes.
The degenerate cases are specially treated and the methods are demanding hav-
ing difficult implementation. The algorithm by Martınez et al. [2009] is based on
classical plane sweep algorithm for computing intersections performing in time
O((n + k)logn). They claim the solution works for general polygons, although
its working for polygons with holes is not demonstrated. Algorithmic details
are given but the implementation issues are not discussed. The implementation
seems difficult and edge overlaps are specially treated. Dimension independence
and robustness issues are not discussed.
2.4 Survey Summary
In this chapter I provided a literature survey of geometric modeling schemes.
The implementation of a geometric data model is should satisfy the following
requirements; (1) closed geometric computations, (2) robustness in geometric
computations, (3) dimension independence, (4) efficient data structure, and (5)
support for level of detail. The surveyed literature shows that the geometric
modeling approaches to date satisfy some of these requirements and none of
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 27
Table 2.1: Summary of algorithms for Boolean operations
Algorithm Operations Scheme Dim Robust Holes Complexity
[Hoffmannet al., 1989]
∩∗, ∪∗, −∗ star-edge 3D yes yes —
[Rivero andFeito, 2000]
∩, ∪, − simplicialchains
2D no yes —
[Martınez et al.,2009]
∩, ∪, − plane sweep 2D no yes O((n + k) logn)
[Karinthi et al.,1994]
∩∗, ∪∗, −∗ —– 2D no yes O(max(2logm, log2 n))
[Milenkovic,1995]
∩∗, ∪∗, −∗ linesegment ar-rangement
2D yes no claimed to be linear butnot proved
[Lauther, 1981] ∩, ∪, − vertical linesweep
2D no yes O(nlogn)
[Gunther,1988b]
∩∗ polyhedralchains
3Dandmore
no yes poly-logarithmic
[Margalit andKnott, 1989]
∩∗, ∪∗, −∗ BRep 2D yes yes O((n + m + k)∗ log(n +m + k))
[Paoluzzi et al.,1993]
∩, ∪, − BSP nD yes no —
[Peng et al.,2005]
∩∗, ∪∗, −∗ Simplex 2D no yes —-
[Smith andDodgson, 2007]
∩, ∪, − BRep 3D yes yes O(nlogn)
[Kui Liu et al.,2007]
∩∗, ∪∗, −∗ Edges 2D no yes ——
[Greiner andHormann, 1998]
∩ Vertex orEdge
2D no yes O(n2)
[O’Rourke,1998]pages:266-268
∩ Sweep line 2D no no O(n logn + k)
[Bentley andOttmann, 1979]
∩ Sweep line 2D no no O(nlogn + k logn)
[Shamos andHoey, 1976]
∩ Sweep line 2D no no O(n2)
[Nievergelt andPreparata,1982]
∩ plane sweep 2D no no O((n + k) logn)
[Vatti, 1992] ∩ scan beam(classicplanesweep)
2D no no ——
[Sabau, 2008] ∩ sweep line 2D no yes ——-
[Weiler andAtherton, 1977]
∩ Vertex orEdge
2D no yes O(n2)
* means regularized operationsn and m stand for number of vertexesk stands for number of intersections
Chapter 2 - GEOMETRIC DATA REPRESENTATION: SURVEY 28
these satisfy all of the aforementioned requirements. Therefore, the need arises
for a geometric data model that provides all the properties, all in one, and that
provides ease of implementation and understanding.
Chapter 3
MODEL REQUIREMENTS
The proposed model for representation of geometry of objects is based on the
basic idea of separation of metric and topological processing as discussed section
1.1 and shown in Figure 1.1 on page 2. The model has been designed keeping in
mind the requirements the model should fulfill. In this chapter I will discuss the
model requirements and proposed solutions to fulfill these requirements.
3.1 The Requirement for Improved Robustness
3.1.1 Improved robustness: The requirement
Current commercial approaches in GIS use floating point numbers and hope that
the resolution is enough to avoid problems due to finite precision arithmetic
used. Computations with floating numbers may produce contradictory results
difficult to deal within numeric processing and particularly difficult for geometric
operations. Geometric calculations in finite precision arithmetic can yield results
that contradict geometric reasoning. Years ago, Franklin in a landmark paper
[Franklin, 1984], has given instructive examples.
Line intersections are most important in GIS (e.g., in overlay computations).
Two lines AB and CD intersect at point P as shown in Figure 3.1. If test the
coordinates of point P , we may find that P is neither on line AB nor on line
CD.
Such cases occur and are caused by the rounding of results to finite precision.
Modern floating point implementations of the IEEE 754 standard, found on most
29
Chapter 3 - MODEL REQUIREMENTS 30
Figure 3.1: Two lines intersect, but the computed point is not on either line
current hardware, have reduced glaring shortcomings. Nevertheless, the numeric
example shown in Figure 3.1 which uses single precision floating point numbers,
demonstrates that the problem occurs because of the rounding errors caused by
fixed precision.
Figure 3.1 shows that testing whether the intersection point P is on the two
lines comes out negative always. Increasing the length of the floating point num-
bers, which increases mostly the precision of arithmetic operations, reduces but
not completely eliminates such vexing problems. Still such problems occur which
are not really found and fixed during debugging. This is especially annoying
concerning computation of overlay operations in GIS [Chrisman et al., 1992].
3.1.2 Proposed solution: Arbitrary precision arithmetic with ho-
mogeneous coordinates
The possible approaches for metric processing are;
3.1.3 Floating point numbers
Floating point numbers approximate real numbers with two restrictions: (1) they
are bounded in the sense that there is a largest and smallest floating point number
of a given representation and (2) computations give approximate results because
resolution is limited. Especially annoying is that results of computations are not
invariant under translation. Results may differ when a configuration is shifted,
the computation performed and the result shifted back compared with the com-
putation directly done.
Chapter 3 - MODEL REQUIREMENTS 31
3.1.4 Integers with special arithmetic
Using integers gives results, which are invariant under shifting; integers are
equally spaced, whereas floating point numbers are concentrated around 0 and
sparser the further away from 0. Care is necessary to remain within the bounds
of standard integers (typically 2.1×109). Intersections of lines do not necessarily
result in an integer but a rational. The proposals for geometry on a grid avoid
the problem shown in Figure 3.1 but they are complex and not easily convincing
that they cover all possible unfortunate inputs.
3.1.5 Big integers
Integers of varying size give arbitrary-precision arithmetic (sometimes imprecisely
called “infinite precision arithmetic” because the number of digits is both finite
and bounded in practice). The precision and the size of these integers are only
limited by the size of the available memory; for this reason they are called “big
integers”. Modern computer software includes routines for arithmetic with big
numbers (packages like big integer, big floating point and big rational supporting
arbitrary precision arithmetic), as this alleviates programs from errors with over-
flow considerations, which if not properly handled, will lead to erroneous results
[Goldenhuber, 1997] or crash.
Computation with big numbers is said to be very slow and complicated [For-
tune and Van Wyk, 1996, Yap, 1997, Li et al., 2005]. Three issues are raised
against the use of big number arithmetic for geometry in GIS:
� difficult implementation,
� slow execution,
� memory intensive
The viability of using big integers depends on the effective performance of geo-
metric operations typical for GIS with big integers. In the next chapter I will
perform experiments to test the viability of using big integers for doing GIS ge-
ometry. I will present the results of an experiment to determine whether solutions
using big numbers are realistic for GIS.
Chapter 3 - MODEL REQUIREMENTS 32
3.1.5.1 Rational numbers with big integers.
The geometric constructions with straight lines can be performed with rational
numbers; the topological relations between regions can be derived precisely from
coordinates. Floating point numbers are necessary if circles, rotations, etc. enter
the picture—but this is seldom for GIS geometry constructions.
A pair of integers forms a rational number; packages to implement regular
arithmetic are available or can be coded easily. With such rational numbers all
geometric computations for geometry with straight lines are possible except for
limitations such as integers that are bounded (there is a smallest and largest
integer for a given size).
Rational numbers represented as pair of big numbers give a domain in which
all of straight line geometry can be done precisely and is limited only by the size
of the available memory.
3.1.5.2 Homogeneous coordinates with big integers.
In case of homogeneous coordinates an n dimensional point is represented by
(n+1 ) – tuple of big integers. Like rational numbers with big integers, all straight
line geometry can be performed precisely with homogeneous coordinates with big
integers, and the precision is limited by the size of the available memory. Ho-
mogeneous coordinates with big integers is a cheaper form of using rationals and
the proposed model will use homogeneous big integer numbers for representing
the coordinates of a point. Homogeneous coordinates form a basis for projective
geometry [Bloomenthal and Rokne, 1994].
The transformations between Cartesian and homogeneous coordinates are:
pc (x, y) → ph(1, x, y)
ph(w, x, y) → pc ( xw ,
yw)
For a point (x,y) in Euclidean plane, its corresponding homogeneous repre-
sentation is (w,x,y). The ordinary plane augmented with points at infinity is
known as the projective plane which can not be represented in Euclidean co-
ordinate system. This is a fundamental reason for using homogeneous coordi-
nates for representing projective planes in projective geometry [Bloomenthal and
Rokne, 1994]. Figure shows the projection of a homogeneous point (w,2,3) for
w = {4,2,1,0.5,0}. As w approaches to 0, the projected Euclidean point moves
away from origin in the direction of (2,3) and at w=0, the point is at infinity.
Chapter 3 - MODEL REQUIREMENTS 33
Figure 3.2: Projection of a homogeneous point to w = 1 plane [Bloomenthal andRokne, 1994]
3.2 The Requirement for Geometric Primitive
3.2.1 Geometric Primitive: The Requirement
The model should use a geometric primitive that offers simple, closed and consis-
tent framework for the representation of geometry of objects. Geometric primi-
tives are the elementary construct used to approximate and represent the spatial
objects digitally. According to “ISO 1907”, a geometric primitive is a geometric
object that is not decomposed further into other primitives in the system. So the
geometric representation of any spatial feature is done by a collection of these
geometric primitives.
In order to represent the geometry of nonconvex regions, a simple and con-
sistent geometric primitive is needed that offers;
� ease of implementation and
� supports consistent topological operations
3.2.2 Proposed solution: Convex polytope
The proposed solution uses convex polytope as the basic geometric primitive,
which is the convex hull in the dual space. Convexity is an important property
of geometric objects and thus the use of convex polytope brings into the model
the benefits associated with. However, most of the objects in real world are
nonconvex, so the model should be able to deal the nonconvex objects. The
proposed solution will tackle the nonconvex objects by decomposing them into
their convex components i.e. convex decomposition and storing the objects in a
Chapter 3 - MODEL REQUIREMENTS 34
hierarchical tree data structure, the convex hull tree. An arbitrary region is then
composed by alternately subtracting convex regions stored in the convex hull tree.
3.3 The Requirement for Closed Geometric Opera-
tions
3.3.1 Geometric Operations: The Requirement
The model and the associated data structure must facilitate the basic geometric
operations. In the GIS domain, the geometric operations that form the basis for
most of the spatial operations are intersection, union, symmetric difference, and
test for point in a polygon etc. Among these operations intersection is the most
important as it is used for overlay operations, the equivalent of join operation
in databases and other operations (e.g. union and symmetric difference) use it.
The model should offer easily understandable and implementable algorithms for
performing the aforementioned operations.
3.3.2 Proposed solution: Intersection of convex polytopes
Intersection of two convex objects is always convex and optimal algorithms for
intersection computation are known. Using the duality principle, the complement
of a region is easily computed. The algorithm for intersection and complement
computation is then used for computing the union and other Boolean operations.
In the next chapter I will analyze the proposed model for supporting the
following closed operations;
� Intersection
� Complement
� Union
� Difference and Symmetric difference
� Point in a polygon
Chapter 3 - MODEL REQUIREMENTS 35
3.4 The Requirement for Level of Detail
3.4.1 Level of detail: The Requirement
Level of detail (LoD) is a fundamental concept in geometric and graphical de-
scription of figures and objects. Level of detail indications are crucial for sharing
geographic information, because only data structures, which can produce just the
right amount of detail that is needed for an application, are capable to serve di-
verse applications. Keeping in view the importance of LoD, the proposed model
should support level of detail.
To construct an LoD data structure and the corresponding progressive algo-
rithm, the following issues must be fixed:
- a description of the original data, with a concept of gradual approximations.
- a simpler description, which is capable of a piecewise gradual approximation
to the original (the primitive).
- a tree in which the primitives are stored.
The selection of the primitive is the major design decision, which influences
the design of algorithms substantially. There are two well known types of dividing
data over LoD:
� by replacement: the next level is a complete and improved description
of what is described above. Examples: linear functions approximated by
linear pieces, strip tree. (in computer vision this is called lowpass pyramid
Perlin and Velho [1995]),
� by delta: the next level is just an addition to what is described in the
higher levels and the improved approximation is achieved by combining the
values from this levels with the sum from above. Example: region quadtree
[Samet, 1984],(in computer vision, this is called bandpass pyramid).
3.4.2 Proposed solution: Hierarchical tree data structure
The model having a hierarchical tree data structure provides support for level of
detail generalization, thus enabling to write progressive algorithms for geometric
computations. The convex hull tree is an LoD structure in which the data is
divided over LoD by delta i.e. the convex hulls in the lower level nodes are
included and smaller than the convex hull at their respective parent node. The
tree evaluation is done progressively by recursive traversal of children trees of
Chapter 3 - MODEL REQUIREMENTS 36
a node only if the children trees are relevant otherwise the children trees are
pruned and not evaluated further. This in combination with lazy evaluation
provides opportunity to design generic progressive algorithms that are efficient
as the computations are done up to the level needed.
3.5 The Requirement for Dimension Independence
3.5.1 Dimension independence: The requirement
The model should be dimension independent i.e. it should model objects of any
dimension. The current models for the representation of geometry are dimension
specific mostly 2D or 3D i.e. they have been implemented to support objects of
a specific dimension, mostly 2D and few for 3D. Even if a model supporting both
2D and 3D would have separate implementations for each dimension. However, a
dimension independent model has a single implementation that works for objects
of any dimension. Dimension independence is important because in some cases we
may need 4D or more models e.g. 3D + time and (or) movement. The proposed
model should be able to deal with objects of any dimension.
3.5.2 Proposed solution: Convex polytope
I will provide a single implementation of proposed model using simplexes that
works for n-dimensions.
Thus the two possible modeling approaches are;
� The model is implemented separately for each dimension e.g. 2D or 3D etc.
� A single implementation that works for any dimension.
The second approach is more attractive, because its provides less code which
means less errors and less maintenance issues. So a single implementation of
representation model is needed that works for n-dimensions.
The necessary topological relations in GIS are deducible from point coor-
dinates by three computations. The dimension independent implementation is
achieved by lifting these basic computations for 2D so that they generalize im-
mediately to 3D and more dimensions.
� Distance between points to decide if two points have the same location.
Chapter 3 - MODEL REQUIREMENTS 37
� Signed areas of a triangle, used to decide if a point is left, in, or right of a
line[Knuth, 1992].
� Testing if a point is inside the circle defined by three other points.
3.6 Requirements Summary
In this chapter I discussed the requirements the model should support and the
proposed solutions to fulfill these requirements. The separation of metric and
topological processing is the key to fulfill the geometric model requirements of
the proposed solution. I discussed the model requirements and proposed solutions
that fulfills the specific requirement. In the next chapter, I will discuss how the
proposed solution, Alternate Hierarchical Decomposition (AHD) achieves these
requirements by integrating the proposed solutions for each requirement.
Chapter 4
SOLUTION: ALTERNATE
HIERARCHICAL
DECOMPOSITION
The solution consists of two parts motivated by the separation of metric and
topological processing as discussed in section 1.1.
1. Homogeneous coordinates and big integers
2. Alternate hierarchical decomposition
4.1 Homogeneous Coordinates and Big Integers
The solution uses homogeneous coordinates stored as arbitrary precision big inte-
gers. A consistent metric processing for geometry of straight lines is then possible
with homogeneous coordinates stored as big integers. The use of arbitrary preci-
sion arithmetic ensures that the (nonrobustness) issues related to fixed precision
arithmetic in geometric computations are avoided. However, the use of big in-
tegers is discouraged in literature as is thought to be computationally expensive
especially for cascaded operations in which case the output of one operation is
the input of next operation. Are big integers really slow? In order to answer
this question and to see the viability of using big integers in the representation
model, an experiment is performed using a cascaded operation of computing the
line intersections.
39
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION40
4.1.1 Experiment for testing the viability of big integers
Line intersection is the primary method to construct new points (besides import-
ing new points from“outside”. Given two lines AB and CD, the intersection point
P is P = (A×B)× (C×D) where × stands for ordinary vector cross on 3-vectors
(note: 2D points are represented as homogeneous 3-vectors).
The behavior of geometry is tested with repeatedly intersecting lines and con-
necting the new points to create new lines, which are then intersected again in
the next iteration. This cascaded intersection computation resembles the genera-
tions of geometrically constructed points. In order to avoid artifacts two different
test programs are produced. Table 4.1 gives a brief of the two test program
approaches.
The assessment included the observation of:
� the ease of programming with big integers,
� execution times compared to the same computation with floating point
numbers,
� and the growth of size of the big integers.
4.1.1.1 Ease of programming
Most modern languages (including the one, Haskell [Jones and Jones, 2003], used
here) offer ready to use packages for big numbers. In Haskell big integers are
implemented as data type Integer, and all ordinary meaningful operations are
overloaded to work on this type.
For both test approaches, the same program — written as polymorph func-
tions — is used in the test for computations with big numbers (Integer) and
floating point (Float or Double). Programming with big numbers is easier than
when working with Float or Integer because no consideration to check for over-
flow is necessary and no hidden bugs may emerge later because the programmer
did not think of the required test, nor are there tolerances to guess.
4.1.1.2 Performance – Approach I
In order to assess the performance penalty caused by the use of big numbers, the
intersection computation is performed for ten cascaded iterations and results com-
pared for the same program executing with big integers and with double precision
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION41
Table 4.1: Test program approaches used in experiment
Approach - I Approach - II
Initial PointGeneration
Random Random
Line Construction(in a generation)
Pairing of points inorder of their positionin the data structure
Pairing of points in order of theirposition in the
Point Construction(in a generation)
Taking lines one by oneand calculating theirintersection with allother succeeding linesin the data structure
Pairing succeeding lines in thedata structure and calculatingtheir intersection point. Eachpair of lines thus gives a singleintersection point to be used inthe next generation
Point Selection (ina generation)
As the number of pointsis increasing per gen-eration, a filter box isused to select only thosepoints for the next gen-eration which fall withinthe boundaries of filterbox
As the number of points is de-creasing per generation, no filteris used. All intersection pointscalculated in one generation areused as input points to the newgeneration
Stopping Criteria Explicitly stopped byspecifying the number ofrounds to compute
Implicitly stopped when thenumber of input points for a gen-eration becomes less than 4.
Test Environment Compiler: Haskell ghc6.8.2 [22] Processor &Memory: 1.2 MHz, 2GBRAM Operating Sys-tem: Ubuntu 7.10
Compiler: Haskell ghc 6.6.1[22]Processor & Memory: P4CPU2.4 GHz, 503MB RAMOperating System: Ubuntu 7.10
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION42
Figure 4.1: Average time per intersection computation in µsec with relatedoperations-approach I
floating point numbers (there is usually no performance difference with modern
computers and hardware between single and double precision float). The differ-
ence was less than expected; total time necessary to compute the intersections of
one iteration was divided by the number of intersections computed. Figure 4.1
shows that the average time for floating point intersections is 1.4×10−4 seconds,
a constant for every iteration. However, computation with big integers takes the
nearly same time for the first 4 iterations and then increasing from iteration 5 up
till > 10−2 seconds for iteration 9.
4.1.1.3 Size of representation – Approach I
In order to answer the question, how quickly the size of big integers grows, the
average length and maximum length of big integers for ten cascaded iterations is
plotted in Figure 4.2.
The result is unfortunately less encouraging: the average length of the rep-
resentation doubles with every following iteration of the cascaded intersection
computation of intersecting lines. From iteration 3 upwards, the size is larger
than the (less precise) floating point numbers. As long as only a small percentage
of points are from constructions, the increase in size is certainly negligible.
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION43
Figure 4.2: Number of digits per coordinate value – approach I
4.1.1.4 Performance – Approach II
Although the implementation using approach II was simpler than approach I,
the results in both cases are comparable as shown in Figure 4.3 and Figure 4.4
(plotted for 700,000 initial random points both big integers and double floating
point numbers).
As shown in Figure 4.3, the average time for double is gradually decreasing,
because the number of intersections to be calculated in subsequent iterations are
decreasing while the average length is constant. Also, the graph shows that the
difference in average time for double and big integers is small in initial iterations
and it increases rapidly in higher iterations. This is because of the fact that the
average length of intersection point coordinates using big integers is increasing
rapidly in higher iterations.
4.1.1.5 Size of representation – Approach II
The behavior of average length of big integers, as shown in Figure 4.4, is same
as for approach I. The average length of double is obviously constant throughout
generations.
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION44
Figure 4.3: Average time per intersection in µsec with related operations-approach II
Figure 4.4: Number of digits per coordinate value – approach II
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION45
4.1.2 Discussion on the viability of big integers for GIS
Programming is easier and straightforward than writing code using conventional
Integer or Float data types. The advantage of metric calculations with big integer
homogeneous coordinates is that they are accurate, not approximates and all
topological relations derived are consistent. Intersection points calculated are
exactly on the line, not only “very near”, etc.
Coarsening of big number homogeneous coordinates is possible, but then the
guarantee is lost. Such schemes for producing point coordinates that are approx-
imate are common in a GIS where points are observed with finite precision. For
example, points constructed in reality as forming a straight line, then measured
and the deviation from a straight line computed yields never the result “exactly
on the line”.
The performance penalty is relatively small and affects only the part of GIS
operations computing metric relations, when points are constructed geometrically
using cascaded intersection of lines. For applications using cascaded construc-
tions, the growth in the size of the big numbers representing the point coordinates
may be problematic. The average size doubles with each cascaded iteration of the
computation. Typical GIS operations do not construct long chains of cascaded
intersections.
The effect on GIS operations is difficult to quantify, but it is most likely small:
� Cascaded constructions of points that go for several iterations are rare in
GIS. Time penalties for such situations are up to 4 or 5 iterations of points—
based on the test results—negligible. The speed reduction depends therefore
on the percentage of “higher iteration points”.
� Only a small amount of the execution time of GIS operations is spent in
geometric operations; assume 5%. If the time for geometric operations
increases by a factor of 10, the total slowing of the program is only 50%.
This slowdown is compensated by the increase in the processing speed of
the processor in about one year (Moore’s law). It is expected that the effect
will be less drastic as all necessary tests for overflow, which use up execution
time, especially in modern pipelined processors, etc., are inessential.
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION46
4.2 Alternate Hierarchical Decomposition
Alternate Hierarchical Decomposition, AHD is an algebraic approach that uses
convex polytope as the basic building block for the geometric representation of an
arbitrary (convex or nonconvex, with or without holes, islands or without islands)
region. The arbitrary region is decomposed into convex hull component regions
(convex decomposition) which are stored in a tree. The arbitrary region is then
formed by alternately subtracting and adding convex hull component regions from
the convex region at the root of the tree. AHD supports operations on regions
which includes the nested traversal trees representing the regions and applying
the corresponding operation on the component convex regions.
4.2.1 AHD: An Overview of the Approach
AHD uses convex polytope as the basic geometric primitive. The use of convex
polytope as the fundamental geometric construct provides simple, straight forward
and consistent topological processing. Convex polytope is defined by the intersec-
tion of finite set of half planes (half spaces in case of 3D or more), however, the
representation of a convex polytope in dual space is used, which is represented by
a set of points and lines, depicting lines and the connection points of half planes
in the Euclidean space respectively. In dual space, a convex polytope is defined
as the convex hull of the dual points [De Berg et al., 2008].
The use of convex polytope as the geometric primitive requires that the non-
convex objects are decomposed into their convex hull components (convex de-
composition). AHD is therefore, based on convex decomposition in which the
nonconvex region is decomposed into convex hull components arranged hierar-
chically in a Convex Hull Tree, CHT (related concepts are convex deficiency tree
[O’Rourke, 1998] and concavity tree [Badawy and Kamel, 2005]), the alternate
levels of which represent convex regions to be added or subtracted depending on
even or odd level respectively. The hierarchical tree data structure is useful for
achieving level of detail.
4.2.2 AHD: An Example
In order to explain the AHD approach, an example that demonstrates the decom-
position steps and shows the output CHT is presented. Suppose a 2-D nonconvex
polytope with a hole (H ) (as shown in Figure 4.5a) is given as input. The ap-
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION47
(a) Input nonconvex polytope with hole H (b) Convex decomposition of input
Figure 4.5: Input polytope and its convex decomposition
Figure 4.6: CHT for demonstration example of Figure
plication of AHD results in the decomposition of the input into its component
convex hulls (Figure 4.5b).
4.2.3 AHD: Convex Hull Tree
The data structure, CHT, is a tree in which every node represents a convex
hull. The root node of the CHT contains the convex hull of the input polytope,
non-leaf nodes contain convex hulls of the nonconvex delta regions (notches) of
the nonconvex structure whose convex hull is represented by parent node. Leaf
nodes contain the convex hulls of the convex delta regions. The positive (+) and
negative (-) signs at alternate levels (even levels are + and odd levels are -) of
the CHT show the convex regions to be added and subtracted for forming the
arbitrary region represented by the CHT (see Figure 4.6).
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION48
4.2.4 AHD: Functions
Two basic functions build and eval are used for populating the data structure
and evaluating the data structure to retrieve the original object respectively.
4.2.4.1 Build function
The AHD build consists of two steps: (1) convex hull computation, and (2) delta
region or regions extraction (using symmetric difference). The process is then
recursively applied to the delta regions until the input region is convex or up to
a specified level of detail ( see section 4.4) specified.
4.2.4.2 Build example
In order to demonstrate the AHD build function, an example is shown to in
Figure 4.7.
4.2.4.3 Build algorithm
The pseudocode for build function is given in Figure 4.8. The pseudocode of the
algorithm for computing the QuickHull [O’Rourke, 1998] is given in Figure 4.9.
4.2.4.4 Build complexity
Given a nonconvex polytope p with n points, the convex hull and delta regions of
p can be computed in O(nlogn) and in liner time, respectively. If r is the number
of delta regions of p then complexity becomes O(rnlogn). However, size of r is
negligibly smaller than n for higher values of n. Thus, the worst case scenario can
not be O(n2logn) and the build function of AHD has an average time complexity
of O(nlogn).
4.2.4.5 Eval function
Given a arbitrary object represented in CHT, the eval function of AHD evaluates
the tree to retrieve the object by recursive traversal of convex hulls stored in the
nodes of the tree up to the level of detail (see section 4.4) specified. The function
eval uses the concept of region merging i.e. given a convex hull parent and its
convex hull components, the components are merged into the parent convex hull
recursively.
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION49
Figure 4.7: AHD steps. L represents the level, dr represents the delta region andch represents the convex hull
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION50
Input: A list pe of polytope edges
1 t = Node [] [] // Initialize t as an empty convex hull tree
2 if pe is empty
3 then return t
4 else if pe contains only a single edge
5 then
6 ch = list of points containing unique endpoints of the
single edge
7 t = Node ch [] // a tree containing single node
8 return t
9 else
10 ch = convex hull of pe
11 che = list of convex hull edges of pe
12 dre = list of delta region edges which are not convex
hull edges
13 dhp = list of delta region points which are on hull
14 dr = list of closed delta regions // region is a list of
edges
15 for all delta regions in dr
16 ct = go to 1 //each region in dr results in a
children tree recursively
17 t = Node ch [ct]
18 return t
Output: A convex hull tree t
Figure 4.8: Pseudocode of the algorithm build that populates the CHT
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION51
Input: A list p of polytope points
1 Initialize ch as an empty list of convex hull points
2 x = point in p with max x-coordinate value
3 y = point in p with max y-coordinate value
4 s1 = list of points in p strictly above line xy
5 s2 = list of points in p strictly above line yx
6 ch = ( [x] + qHull (x, y s1) + [y] + qHull (y, x, s2))
7 return ch
Output: A list ch of convex hull points
Input: Two end points a and b of a line ab and a list S of points above that lineab
1 Initialize sh as an empty list of sub hull points
2 if S is empty
3 then return []
4 else
5 c = farthest point of S from line ab
6 partA = points in S strictly above line ac
7 partB = points in S strictly above line cb
8 chA = qHull (a, c, partA)
9 chB = qHull (c, b, partB)
10 sh = chA + [c] + chB
11 return sh
Output: A list of sub hull points sh
Figure 4.9: Pseudocode of QuickHull algorithm and its helping function qHull,which computes the sub hulls of partitions recursively
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION52
4.2.4.6 Eval example
An example is shown in Figure 4.10 to illustrate the eval function that using
region merging for forming the object stored in CHT. Four scenarios are possible
for merging a child hull into a parent hull:
1. An edge of child hull completely overlaps an edge of parent hull i.e. both
edges have same starting and ending points. For example, in Figure 4.10f
e12 of region represented by child hull 1 edge e03 of parent region overlap.
In this this case the common edge is removed from both child and parent
regions. The orientation of the remaining edges of the region represented
by the child tree is reversed (as shown in Figure 4.10g).
2. An edge of child hull overlaps an edge of parent hull in the middle i.e.
starting and ending points of the child edge lie on the edge of the parent
edge and are not equal to the starting and ending edges of the parent edge.
For example, in Figure 4.10f e24 of child region lies on the edge e01 of the
parent region. In this this case e24 is removed from child region while
the edge e01is split into two edges e011 and e012 . Also, the orientation of
the remaining edges of the region represented by this child tree is reversed
(Figure 4.10g).
3. An edge of the child hull overlaps an edge of the parent hull but both edges
have at least one common point. For example, edge e52 of child region
overlaps the edge e02 of the parent region and both have a common starting
point. In this case a new single edge e021is created excluding the overlapping
portion. The orientation of the remaining edges of the region represented
by this child tree is reversed (Figure 4.10g).
4. There is no any edges overlap i.e. all the points of the child hull are inside
the parent hull. This shows that the child hull represents a hole and the
merging is done just by reversing the edges of region represented by child
tree. For example, in Figure 4.10f child region 3 represents a hole region.
4.2.4.7 Eval algorithm
The pseudocode of the algorithm for complete evaluation of a CHT, which re-
trieves the original object, is given in Figure 4.11.
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION53
(a) An input nonconvex region with ahole
(b) Decomposed input
(c) CHT of input re-gion
(d) Merging child hull 1 into parent hull0-find common edge
(e) Merging child hull 1 into parent hull 0-remove common edge from both hulls andreverse orientation of remaining edges ofthe child hull 1
(f) Merging all children of parent hull 0 (g) Output of “eval”
Figure 4.10: AHD “eval” using merging of convex hull regions
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION54
Input: A convex hull tree t, Node x xts : x is convex hull at root node and xtsis list of children trees
1 Initialize pe, cr as empty list of edges
2 if t is empty
3 then eturn pe
4 else if xts is empty // a leaf node
5 then
6 pe = list containing edges of convex hull x
7 return pe
8 else
9 pe = list containing edges of convex hull x
10 for all child trees in xts
11 cr = cr + go to 1 // each children tree in xts results in a
component region recursively
12 mr = merge regions in cr into the region pe
13 pe = mr
14 return pe
Output: A list pe of polytope edges
Figure 4.11: Pseudocode of the algorithm eval that processes the CHT to extractrepresented region
4.2.4.8 Eval complexity
For a arbitrary simple polygon with n vertexes;
Max number of possible concavities or notches = (n-3)+1 ≈ n if n is too large
Max nodes of a CHT = (n-3)+1 ≈ n if n is too large
If the complexity of the algorithm for point in a convex hull is logn then the
worst case complexity of eval is n logn
4.2.5 AHD: Characteristics
AHD uses convex decomposition for representation of the object in CHT. In sec-
tion 2.1.2.2, the criteria for the classification decomposition techniques is men-
tioned. Most of the algorithms for decomposition are optimized for one of the
these criteria [Kriegel et al., 1991] depending on application needs. Table 4.2
shows a characterization of AHD based on these criteria.
The important characteristics of AHD as shown in Figure 4.12 are;
1. Closed operations: AHD supports Boolean operations which are closed
and implementation is efficient because of hierarchical data structure.
2. Dimension independence: AHD is dimension independent as a single
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION55
Table 4.2: Characterization of the decomposition technique
Criteria AHD
Type of components Convex
Type of input Polytope with or without holes
Relation of components Covering- as the component convex hulls at level l+1are
contained within the convex hull at level l
Objective function Does not fall to any category, but the objective is to
achieve generalization of data modeling framework using
convex polytopes
Input dimension Dimension independent
Figure 4.12: AHD characteristics
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION56
(a) (b)
Figure 4.13: Approaches for dimension independent AHD
implementation of it works for all dimensions by lifting basic geometric
primitives.
3. Level of detail: AHD supports level of detail because of hierarchical
data structure that facilitates generalization and design of progressive al-
gorithms.
4.3 Dimension Independent AHD
In this section the dimension independence of AHD is shown by providing a
dimension independent convex decomposition algorithm. The AHD approach
has two steps:
1. Convex hull computation
2. Delta region (s) or notch (s) extraction
For a dimension independent AHD, the two steps need to be dimension indepen-
dent. There are two possible approaches as shown in Figure 4.13. First approach
deals both of the aforementioned steps independent of each other (Figure 4.13a).
This means any of the existing dimension independent convex hull algorithm [Bar-
ber et al., 1996, Karimipour et al., 2008] can be used and the result is transformed
to the object representation model before extracting the delta regions. In second
approach, a data representation model is used that supports both operations to
be performed in sequence without conversions (fig. 3b). The proposed approach
for n-Dimensional AHD is based on the second approach which is better as it
avoids unnecessary conversions.
The approach is based on the simplex based data representation model which
provides a consistent implementation of the two basic steps of AHD i.e. convex
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION57
Figure 4.14: n-dimensional AHD- An example for 3D case
Figure 4.15: Simplexes of dimensions 0 to 3: node, edge, triangle, and tetrahedron
hull computation and delta region extraction. Figure shows the decomposition
example for a 3D case.
A n-simplex Sn (abbreviated as simplex henceforth) is defined as the smallest
convex set in the Euclidean space that contains n + 1 of n-dimensional vertexes
v0, . . . ,vn that do not lie in a hyperplane of dimension less than n [Alexandroff
and Aleksandrov, 1961]. A simplex is represented by the list of its vertexes as
Sn =< v0, ... ,vn >. Figure 4.15 shows simplexes of dimensions 0 to 3.
4.3.1 Pseudocode: n-Dimensional AHD
Dimension independent convex hull computation is the most important step for
n- dimensional AHD. In this section the pseudocode of the algorithms for Inc-
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION58
Input: A set of n−dimensional points LP= {p0, ..., pM}1 I = an n-simplex constructed by LP
2 E = boundary of I (which consits of n+1 of (n-1)-simplexes)
3 for all points p in {pN+1, ..., PM}
4 if p is outside the region delimited by all (n-1)-
simplexes in E
5 S = set of all (n-1)-simplexes in E that make a cw
turn with respect to point p
6 B = set of (n-2)-simplexes that are on the border of
S
7 N = set of (n-1)-simplexes constructed by adding p to
all (n-2)-simplexes b in E
8 E = E-S+N
9 return E
Output: A List E contains the (n-1)-simplexes that construct the convexhull ofthe input
Figure 4.16: Pseudocode of dimension independent incremental convex hull algo-rithm
Convexhull and n -dimensional AHD are provided.
4.3.1.1 Convex hull Computation
The convex hull of a set of points is the smallest region that contains the points.
An incremental algorithm called IncConvexhull [Barber et al., 1996, De Berg
et al., 2008] is used for the convex hull computation. For some 2D points, the In-
cConvexhull algorithm starts with the triangle goes through the first three points,
which is their convex hull. Other points are inserted one by one into the construc-
tion and after each insertion, the convex hull is modified: if the inserted point
is inside the convex hull, no change is needed in the configuration of the current
convex hull; If it is outside, however, its opposite edges are dropped from the con-
vex hull and new edges are added, which are constructed by adding the inserted
point to the border points of the dropped edges. Extension of IncConvexhull
algorithm to n-dimensional points is possible by using the concept of simplexes
as shown by Karimipour et al. [2008]. Figure 4.16 shows the pseudocode of the
n-dimensional convex hull algorithm.
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION59
Input: A set of n−simplexes R = S0, ...,Sm
1 if sptR is empty
2 then return T
3 where
4 T = Node (chR , [])
5 else return T
6 where
7 T = Node (chR , apply build on all elements of the
list sptR)
8 sptR = split (symD(R, chR))
9 chR = Convexhull (vs)
10 vs = union of vertices of input n-simplexes R
11 symD(x,y) = (x-y)+(y-x)
Output: A tree of signed convex regions that construct R
Figure 4.17: n-Dimensional AHD
4.3.1.2 n-Dimensional AHD
Using the concept of simplexes, the AHD algorithm is implemented for polytopes
of any dimension. The dimension independent convex hull algorithm of Figure
4.16 and the split function of Table 6.2 are used for a dimension independent
AHD, the pseudocode of which is shown in Figure 4.17.
4.4 AHD: Level of Detail
AHD supports level of detail because of the hierarchical tree data structure. Both
the build and eval functions of AHD use the level of detail for building the CHT of
an arbitrary region and for evaluating the CHT for forming the arbitrary region
up to the specified level respectively. AHD therefore provides a level of detail
approximation of an arbitrary region by convex decomposition of the region into
its convex hull components. The component convex hulls are always smaller and
overlap the parent convex hull.
Figure 4.19 demonstrates the level of detail approximation of the given 2D
nonconvex polygon using AHD (Figure 4.18).
The support for level of detail facilitates in designing progressive algorithms
for performing Boolean operations. The level of detail together with lazy evalua-
tion (as supported by some programming languages e.g. Haskell) leads to efficient
geometric algorithms for performing operations like point in polygon test, inter-
section, union and so forth. The LoD in AHD facilitates reduced number of
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION60
(a) An example nonconvexpolygon
(b) Convex hull components ofexample polygon
(c) CHT of input polygon
Figure 4.18: AHD example for a nonconvex polygon
computations. Special cases are tested for Boolean operations like intersection,
inclusion, etc., which are best treated by operations with understanding of the
geometric data structure used. For example, test for the “intersection detection”
for two polygons represented by AHD stops at the first level if the convex hulls of
the given polygons do not overlap further exploration is unnecessary. Similarly,
the test for “point in polygon” stops at the first level, if the ’test point’ is not
inside the convex hull of the given polygon.
The progression in such algorithms is controlled either by specifying level ex-
plicitly or specifying an error bound for the result. The LoD in AHD is based
on a geometric criterion, i.e., the difference between the convex hull and the true
region. Applications may need to indicate the required quality of the approxima-
tion, for example, in case of AHD different criteria may apply;
� percent error in area (volume),
� absolute error in area (volume),
� maximum absolute error as the distance of boundary points from the ap-
proximate boundary.
The advantage of AHD compared to a region quadtree is (a) it represents exactly
polytopes of n-dimension bounded by n-flats, and (b) it is invariant to shift or
rotation of a coordinate system.
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION61
(a) LoD0 region approximation and the corresponding CHT node is highlighted
(b) LoD1 region approximation and the corresponding CHT nodes arehighlighted
(c) LoD2 region approximation and the corresponding CHT nodes arehighlighted
Figure 4.19: The concept of LoD in AHD for example polygon in Fig. 4.18. Theincreasing level of detail from LoD0 to LoD2 increases the quality of approxima-tion
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION62
4.5 Solution Summary
4.5.1 Robustness
The solution avoids the problem of numerical nonrobustness [Franklin, 1984,
Mehlhorn and Naher, 1999, Li et al., 2005] by using exact algebra. This has
been achieved through the use of big integers allowing arbitrary precision arith-
metic. Thus, topological inconsistencies are avoided, as the underlying algorithms
have been implemented keeping in mind the arbitrary precision arithmetic for nu-
merical computations involved.
4.5.2 Closed Boolean operations
The solution supports Boolean operations. The computation of Boolean opera-
tions (e.g. intersection between two regions) is done by recursively traversing the
trees and applying the intersection operation on the convex hulls. The operations
are discussed in detail in next chapter.
4.5.3 Dimension independence
To the best of my knowledge, the algorithms for decomposition problem are not
dimension independent as most work for 2-D and some deal with 3-D objects.
Although for demonstration purposes 2-D polytopes are used, the solution is
extendable to n-dimensional polytopes. This is achievable if the algorithm for
computing convex hull is dimension independent. For this any of the existing
dimension independent convex hull algorithms [Barber et al., 1996, Karimipour
et al., 2008] can be used.
4.5.4 Level of detail
Depending on the applications, the complete decomposition of nonconvex poly-
gons to convex polygons may not be needed [Lien and Amato, 2006]. In such cases
approximations are used by halting the decomposition process using some thresh-
old.The solution is flexible in that the decomposition process can be stopped at
any level depending on application thus allowing level of detail approximation.
Chapter 4 - SOLUTION: ALTERNATE HIERARCHICAL DECOMPOSITION63
4.5.5 Number of components
The solutions to the problem of decomposing a nonconvex polygon into its optimal
(minimum) number of convex components vary depending on whether Steiner
points allowed [Lien and Amato, 2006]. Although the solution was not aiming to
optimize the number of resultant components, however the solution is not as bad
as the mostly used decompositions for near optimal number of components. The
number of components is the number of CHT nodes, which is 13 and partitioning
the example in Figure 4.5 results in 14 convex components.
Chapter 5
AHD BASED GEOMETRIC
OPERATIONS
Boolean operations are fundamental to geographic information processing. In
GIS, spatial operations like convex hulls, buffering, distances and Boolean or set-
theoretic operations (intersection, union difference and symmetric difference) are
extensively used for extracting useful spatial information out of stored spatial
data.
The chapter shows that the AHD approach allows optimal (or near optimal)
computation of Boolean and other geometric operations (point in a polygon and
line intersecting a polygon) using AHD approach. In a GIS a separate data
structure is not possible for every operation but the same structure can be used
for all operations.
5.1 Intersection Operation
The intersection operation is of fundamental importance [Shamos and Hoey, 1976]
as other Boolean operations union, difference and symmetric difference can be re-
duced to the computation based on intersection and complement. Also, it is the
most expensive operation computationally, roughly taking 80 % of the running
time of the program [Greiner and Hormann, 1998]. The intersection operation
has two problems, intersection detection and intersection computation. The in-
tersection detection between two convex objects is a basic geometric operation
[Chazelle and Dobkin, 1987] and a description of the topic is provided by Mount
[1997]. The intersection detection problem is not considered that itself is a chal-
65
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 66
Figure 5.1: Different polygon intersection scenarios
lenging problem and algorithms exist in literature (e.g. [Shamos and Hoey, 1976,
Chazelle and Dobkin, 1987]). For simplicity it is assumed that the objects under
consideration for intersection computation do intersect.
The Problem: Given two arbitrary regions (convex or nonconvex,
with non intersecting edges, and with or without holes), compute
the intersection region that may be;
(a) Empty (Figure 1a)
(b) Convex (Figure 1b)
(c) Nonconvex (Figure 1c)
(d) A set of convex and (or) nonconvex disjoint regions (Figure 1d)
5.1.1 AHD based intersection computation
The intersection computation consists of three operations; 1) Intersection between
two convex objects. 2) Intersection between a convex and a CHT (nonconvex ob-
ject). 3) Intersection between two CHTs (two nonconvex objects). This gives for
(1) the basic and simple operation of intersection computation between two con-
vex hulls. For (2), the traversal with basic operation in (1) and for (3), the CHT
traversal with operation in (2). Only the basic operation of intersection of two
convex hulls is geometric (for which well known algorithms exist e.g. [O’Rourke
et al., 1982]) and the other operations are repeated application of this by travers-
ing tree structures.
Figure 5.2 shows the intersection computation approach. The basic intersec-
tion operation between two convex hulls, convex-convex intersection (CCINT) is
used to compute the convex-nonconvex intersection (CNINT) between a convex
and nonconvex object. The CNINT operation is in turn used for computing a
generic intersection operation (GINT).
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 67
Figure 5.2: Intersection computation approach
(a) Two convex polygons (b) CHTs of A and B (c) Resul-tant CHTof A∩B
Figure 5.3: Convex-convex intersection computation
5.1.1.1 Basic intersection operation: Convex-convex intersection
The basic operation is closed under intersection and the result is always a con-
vex region. The intersection computation between two convex polygons is the
basic operation (for which known solutions exist e.g. [Toussaint, 1985, Chazelle,
1989]). Any of the existing solutions can be used but an algorithm for intersection
computation between two convex polygons using convex hulls is provided.
Figure 5.3a shows two intersecting convex objects. The intersection is per-
formed by computing the intersection of the convex hulls ( (a1 and b1) of the two
convex regions as shown by two leaf nodes of CHT in Figure 5.3b. The result
is a convex hull (a1∩b1) region represented by a leaf node of CHT as shown in
Figure 5.3c.
Pseudocode of convex-convex intersection algorithm: The pseudocode
of the basic intersection operation between two convex polygons is shown in Figure
5.4. The algorithm uses the QuickHull algorithm [O’Rourke, 1998] for convex
hull computation and compIntPoints routine computes the intersection points by
pairing the intersecting delta edges.
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 68
Input: Two convex polytopes poly1 and poly2 in vertex notation
1 Initialize set of intersection region points I as an empty
set
2 ch = convex hull (poly1 + poly2)
3 de = (poly1 edges + poly2 edges) - ch edges
4 ip = compIntPoints de
5 irp = (poly1 -poly2) + (poly2 -poly1) + ip
6 I = convexHull irp
Output: Set of intersection region points I (always a convex hull)
Figure 5.4: Pseudocode of convex-convex intersection
Complexity of convex-convex intersection algorithm: For two convex
objects with number of points m and n, the worst case time complexity of the
algorithm is O((m + n) log(m + n)).
An Example of convex-convex intersection: The algorithmic steps of
basic intersection are shown in the example in Figure 5.5.
5.1.1.2 Convex-nonconvex intersection
The process of intersection computation between a convex and a nonconvex object
consists the repeated application of the basic operation, CCINT, between the
convex hull of the convex object and the convex hull components of the nonconvex
object that are traversed recursively in CHT. The resulting convex hull of each
basic operation is then placed in the resultant CHT at the position same as
the position of component convex hull of the nonconvex polygon involved in the
basic operation. The further recursive evaluation of children trees of a component
convex hull is stopped if the basic operation involving that component convex hull
is null. The support for LoD together with lazy evaluation avoids unnecessary
computations reducing the overall number of operations. For example, Figure
5.6a shows a convex polygon A and a nonconvex polygon B and Figure 5.6b shows
the decomposed inputs with their components numbered. Figure 5.6c shows the
data structure representation of the input polygons A and B. The process starts
by the computation of basic the operation between convex hull of convex polygon
A (that is a1 ) and the convex hull of the nonconvex polygon B (that is the root
b1 ). Since b1, the component convex hull of the nonconvex polygon is at root the
result of basic operation between a1 and b1 will form the root of the resultant
CHT as shown in Figure 5.6e. The process is then recursively applied to the
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 69
(a) Input convex polygons (b) Convex hull of (poly1 + poly2)
(c) Delta edges (d) Intersecting edges
(e) Poly1-poly2 (f) Poly2-poly1
(g) Computed intersection points (h) Intersection region
Figure 5.5: Convex-convex intersection example
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 70
(a) Input polygons (b) Decomposed inputs
(c) CHTs of input poly-gons
(d) Intersection of convexhull of A, a1 and the CHTof B
(e) CHT of intersection
(f) Intersection CHT compo-nent hulls
(g) Intersection region (A∩B)
Figure 5.6: Convex-nonconvex intersection computation
children trees of b1. Since the intersection between convex hull a1 and b3 is null,
the child tree of b3 containing convex hull b4 will not be evaluated further (Figure
5.6e). The resultant intersection is a nonconvex region that is represented by two
convex hull components as shown in Figure 5.6f. The evaluation (merging) of the
resultant CHT, containing two component convex hulls results in the resultant
intersection region as shown in Figure 5.6g.
Since the intersection of a convex and a nonconvex polygon may be a set of
disjoint convex or (and) nonconvex regions (e.g. Figure 5.12f ), the resultant CHT
(if it does not represent a single convex region) is further evaluated for simplifi-
cation, that results in multiple CHTs each representing the disjoint intersection
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 71
Input: A convex polytope poly1 and a nonconvex polytope poly2 (CHT): poly1= ch and poly2 = Node x xts where x is the convex hull of poly2 and xts is theset of children CHTs
1 Initialize I as an empty tree
2 if xts is empty then
3 t = Node (intersection of ch and x) [ ]
4 I = set of trees containing only t
5 else
6 ct = map (CNINT ch ) xts // set of trees by
recursively traversing xts with (CNINT ch)
7 it = Node (CCINT of ch and x) ct
8 eir = process it // set of edges of intersection
region obtained by processing it
9 ir = splitRegions eir // set of disjoint regions
obtained by separating closed regions
I = build ir //set of trees obtained by
building each region in ir
Output: Set of disjoint intersection regions I represented in CHT
Figure 5.7: Pseudocode of convex-nonconvex intersection algorithm
region.
Pseudocode of convex-nonconvex intersection: The pseudocode of the
CNINT is shown in Figure 5.7.
Complexity of convex-nonconvex intersection algorithm: For a con-
vex object with number of points m and a nonconvex object with number of
points n and notches p, the worst case time complexity of the algorithm is
O(p((m + n) log(m + n))).
An Example of convex-nonconvex intersection: The example in Figure
5.8 shows the algorithmic steps for convex-nonconvex intersection computation.
5.1.1.3 Generic intersection operation
The intersection operation for a convex and a nonconvex object, CNINT, is then
used to compute the intersection between two nonconvex objects represented
by CHTs. Suppose two intersecting nonconvex polygons A and B as shown in
Figure 5.9a. Figure 5.9b shows the decomposed convex hull components of both
nonconvex polygons and their CHTs are shown in Figure 5.9c. The intersection
process starts by taking the convex-nonconvex intersection between convex hull
of nonconvex polygon B (b1 ) and the CHT of the nonconvex polygon A (Figure
5.9d).
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 72
(a) Input polygons (b) Decomposed inputs
(c) a1 and b1 (d) a1∩ b1
(e) a1 and b2 (f) a1 ∩ b2
(g) a1 ∩ b3 = NULL (h) A ∩ B
Figure 5.8: Convex-nonconvex intersection example
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 73
(a) Two nonconvex polygons Aand B
(b) Decomposed components ofA and B. A∩B is highlighted
(c) CHTs of Aand B
(d) Intersectionof A and b1 (con-vex hull of B)
(e) ResultantCHT of (d)
(f) a2∩b1is null
(g) Intersect the re-sult of (f) with chil-dren of node b1 recur-sively
(h) Resultant CHTof (g)
(i) Result of h as b1∩b2 = b2is null
(j) Add or graft resultof i to the its parentintersection CHT
Figure 5.9: Generic intersection operation
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 74
Since the intersection between b1 and a2 convex components is null (Figure
5.9e ), the result is a single node CHT (Figure 5.9g) having the convex hull
(a1∩b1) which is then intersected recursively with all the children trees of convex
hull b1 of the nonconvex polygon B. Since, b1 has only one child (which is also
a single node representing a convex hull) the intersection of (a1∩ b1) and child
tree (Figure 5.9h) results in a convex hull (a1∩b2) represented in a single node
CHT as shown in Figure 5.9i and Figure 5.9j. The resultant CHT containing
(a1∩b2) is then added (grafted) back to the parent intersection CHT containing
single node (a1∩ b1). Thus, the result is a nonconvex region shown as shaded
in Figure 5.9b and the resultant CHT shown in Figure 5.9j. If the result of
convex-nonconvex intersection is a set of disjoint regions (a set of CHTs, each
representing a region) then for each CHT of the disjoint region, the same process
is repeated independently of the other CHTs representing the disjoint intersection
regions.
Pseudocode of generic intersection operation: The algorithm for generic
intersection operation incorporates the previous two intersection operations. It
is generic in the sense that it computes the intersection of two polygons (convex
or nonconvex and with or without holes). The pseudocode of the GINT is shown
in Figure 5.10.
Complexity of generic intersection algorithm: For two arbitrary regions
with number of points m and n and notches q and p respectively, the worst case
time complexity of the algorithm is O(qp((m + n) log(m + n))).
An example of generic intersection operation: An example of generic
intersection operation between two nonconvex polygons is shown in Figure 5.11.
5.1.2 Intersection special cases
The intersection of two polygons may be a point (Figure 5.12a), a line (Figure
5.12b), polygon or a combination consisting of the three (Figure 5.12f). In cases
where the intersection region is not a polygon or a set of polygons, special treat-
ment is needed. A variety of special cases for the basic intersection have been
shown by O’Rourke [1998] and few are shown in Figure 5.12.
Since the approach is based on basic intersection operation, there is no need
to tackle the special cases for convex-nonconvex or nonconvex-nonconvex inter-
section operations. If the special cases are tackled for basic operation, other
operations based on it will also tackle the special cases. For special cases like
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 75
Input: Two polygonal regions (convex or nonconvex) in CHT notation: poly1 =Node x xts and poly2 = Node y yts where x and y are the convex hulls and xtsand yts are the set of children trees of polygon 1 and polygon2 respectively
1 Initialize I as an empty tree
2 if xts is empty //poly1 is convex
3 then return CNINT x poly2
4 else if yts is empty //poly2 is convex
5 then return CNINT y poly1
6 else
7 oi = CNINT y poly1
8 if oi is empty
9 then return I
10 else for every tree tx in oi
11 for every tree ty in yts
12 ct = compute GINT tx ty
13 nt = addtrees ct in tree tx
14 I = addtree nt into I
15 return I
Output: Set of disjoint intersection regions I in CHT
Figure 5.10: Pseudocode of generic intersection algorithm
complete or partial overlap, the approach do not need any special treatment.
However, for cases when the intersection is a point or a line, special treatment is
needed which is achieved by making changes not in the main intersection algo-
rithms but in the QuickHull module and the eval function of the AHD. Table 5.1
lists special cases and whether these cases are specially treated in the algorithm
or not.
5.1.3 Intersection solution characteristics
The proposed solution to perform Boolean intersection operations on general
polygons has following properties;
1. The approach is robust as it uses homogeneous big integer coordinates for
representing points. Robustness is an important issue in geometric compu-
tations [Mount, 1997]. Most of the algorithms, as discussed in section 2.3
on page 24 , do not cater the robustness issue. Only few address the issue
[Mount, 1997, Smith and Dodgson, 2007].
2. The approach is based on the basic intersection operation between two
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 76
(a) Input nonconvexpolygons
(b) Decomposed inputs (c) CNINT: CHTof A and convexhull of B
(d) Result of (c)
(e) Result in (c) issimplified and twodisjoint regions result
(f) r1 ∩ first child ofB
(g) Result of (f) (h) Remove (g) fromr1
(i) Result in (h) ∩ sec-ond child tree of B
(j) r2 ∩ first child of B (k) Result of (j) (l) Remove result in(k) from r2
(m) Result in (l) ∩second child tree of B
(n) Result of (m) (o) Remove (n)from (l)
(p) Disjoint intersec-tion regions
Figure 5.11: Nonconvex-nonconvex intersection example
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 77
(a) (b) (c)
(d) (e) (f)
Figure 5.12: Some special cases for intersection computation
Table 5.1: Special case treatment
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 78
(a) A given region R (b) Complement of given region R
Figure 5.13: Complement R of a given region R
convex polygons. Thus it exploits the benefits associated with convexity.
3. The tree data structure of AHD supports LoD and allows reduced number
of computations. The children trees of a node are evaluated only if the
intersection is not null. The algorithm proceeds recursively in a branch
of CHT only if the intersection is not null. Rule is “an object A can not
intersect with the component hulls of object B, if convex hull of object A is
not intersecting with the convex hull of object B”.
4. In some applications, the result of Boolean operation may be needed to
be decomposed. For example, Vatti [1992] mentions a case when decom-
posed output is useful. So the approach is useful for applications where the
decomposed output is needed in the form of convex components.
5. The special cases need no special treatment like overlapping edges etc. Some
require few modifications like when the result has dangling edges or points
etc.
6. The approach is dimension independent.
5.2 Complement Computation
Is it possible to compute the complement of a region represented in a CHT data
structure using AHD approach? Since AHD can deal regions with holes, the
complement of a region is just an open region with the given region as a hole as
shown in Figure 5.13.
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 79
At data structure level, the complement of a given region R represented by
AHD is simple. The CHT of R is grafted as to a special node called universal
node, that represents the convex hull of the universal region (which is an open
region) which makes region R a hole of the universal open region. For the sake of
simplicity and understanding assume that the universal region is finite and large
enough to contain the given region or regions.
The complement of a polygon represented by AHD is computed in constant
time as it includes the addition or removal of the universal node from the CHT
representing the polygon.
5.3 Union Computation
Intersection of two arbitrary shape regions represented using AHD is a simple,
efficient and provable algorithm because the intersection of two convex regions
is always convex. Unfortunately, the union of two convex regions is not always
convex and thus no straightforward algorithm for union has been found.
Duality is an approach that allows to replace a complicated algorithm into a
simpler dual one using the equation;
dual(a◦b) = (dual a) ⊗ (dual b)
where
⊗ is the dual operation to ◦ and
dual.dual=id
Here a convincing approach for union computation is presented that exploits
duality. The approach shows that eventually the same algorithm (for intersec-
tion) can be used to compute union with the same low complexity. From school
mathematics (De Morgan’s laws) the intersection and union operations are dual
of each other with respect to complement:
A∩B = A∪B
A∪B = (A∩B) (5.1)
Since the union is computed using the intersection algorithm, it has the same
complexity as the intersection algorithm.
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 80
Figure 5.14: A4B
5.4 Difference and Symmetric Difference Computa-
tion
In mathematics, the symmetric difference (logical XOR) of two sets is the set
of elements which are in either of the sets and not in their intersection. The
symmetric difference of two regions A and B, denoted by A4B, is the region
containing both of the regions A and B (AUB) excluding the intersection region
(A∩B) as shown in Figure 5.14.
The symmetric difference of two regions is defined as;
A4B = (A∪B)\ (A∩B) (5.2)
A4B = (A\B)∪ (B\A) (5.3)
So symmetric difference of two regions results in two disjoint regions, one
representing A\B and other B\A.
Like union, the closure property also do not hold for difference operation.
Like for union operation, the duality is used for computing the difference using
the intersection operation for computing A\B as shown in Figure 5.15.
We know that:
A\B = A∩B (5.4)
B\A = B∩A (5.5)
Using equations 5.1 and 5.4, the equation 5.2 for A4B can be rewritten solely
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 81
(a) Two intersecting regionsA and B
(b) Complement of region B
(c) Compute A∩B (d) A\B = A∩B
Figure 5.15: Difference computation (A\B example)
in intersection and complement terms as;
A4B = ((A∩B))∩ (A∩B) (5.6)
Using equations 5.1, 5.4 and 5.5 the equation 5.3 for A4B can be rewritten
solely in intersection and complement terms as;
A4B = (A∩B)∩ (B∩A) (5.7)
Either of the equations 5.6 or 5.7 can be used for the implementation of
symmetric difference computation.
5.5 Line-Region Intersection Computation
The problem of determining the intersection of a line segment and a region (as
shown in Figure 5.16a) is important for computations like removal of hidden faces,
collision detection, Boolean operations etc [Segura and Feito, 1998].
For the computation of line-region intersection, another data structure called
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 82
segment tree is introduced. The segment tree is an arbitrary tree, like CHT, every
node of which represents a list of line segments. Figure 5.17a shows the CHT of
region R shown in Figure 5.16a.
The algorithm for the line-region computation uses the algorithm for com-
putation of a line with a convex hull. The algorithm starts by computing the
intersection of L and the convex hull of region R which is the root of the CHT
represented as node 0. If the result is empty meaning line is not intersecting the
given region, then the algorithm stops. The result of L∩0 which is a line segment
L0 as shown in Figure 5.16c is stored as the root of the segment tree (Figure
5.17b). The process is then repeated for every child tree of CHT.
LoD enables that sub-trees of a node in CHT are not evaluated if the inter-
section of line L and the convex hull represented by that hull is empty e.g. nodes
L3 and L5 are empty and will not be further evaluated even if these would have
children trees. Figure 5.17 shows the intersection computation steps, the result-
ing component line segments are shown in Figure 5.17d and Figure 5.17c shows
the corresponding final segment tree.
For the output generation the segment tree is traversed in preorder. The
evaluation involves the subtraction of a line segment from a given line. The line
segments in the children nodes of a given node are sequentially subtracted (from
left to right) from the line segments represented by that node. This is done
recursively on the children trees of a node and is shown in Figure 5.17e. Root
contains segment L0 from which evaluated left most child tree is to be subtracted
first. Since it is a single node containing segment L1, simply subtract it from
which results in segment L01 updating node 0. Now second sub-tree is evaluated
which means L4 is subtracted from L3 resulting in two segments L34 which are
then subtracted from L01 at node 0 resulting in L34 and the root node 0 is updated
with the result. The single node then contains set of segments which are the result
of the intersection computation.
5.6 Point in Polygon Test
The hierarchical data structure of AHD offering LoD facilitates the spatial search-
ing problems. This is shown with point in polygon (PIP) test which is a special
case of point location problems and is most basic of spatial operations ([De Smith
et al., 2007], page 86).
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 83
(a) A Line L intersecting a regionR
(b) Convex Decomposition of re-gion R
(c) L0 = L∩0 (d) L1 = L∩1 (e) L3 = L∩3 (f) L4 = L∩4
Figure 5.16: Line-region intersection
(a) CHT of region Rand a single segmentnode containing theline L
(b) Line segment tree (c) Line segmenttree
(d) Component linesegments
(e) Segment tree manipulation
Figure 5.17: Line-region intersection steps
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 84
(a) A nonconvex region and a point forPIP test
(b) Decomposed input polygon for PIPtest
(c) CHT of input polygon for PIP test
Figure 5.18: AHD based point in a polygon test
5.6.1 The PIP problem
Given a region R represented by an arbitrary polygon and a point p as shown
in Figure 5.18, the problem is to determine whether p lies inside or outside of
region R.
5.6.2 AHD based PIP computation
Given a region R represented by AHD and a point p (Figure 5.18) a simple and
efficient algorithm for PIP test is provided. The convex decomposition of the given
region R is shown in Figure 5.18b while Figure 5.18c shows the corresponding
CHT data structure. The algorithm for PIP test utilizes the point in convex
polygon (PICP) algorithm for which efficient solutions already exist.
Because of the LOD structure of AHD, determining that the given point is not
inside a given region is immediate, this needs just testing the point p’s inclusion
in the convex hull of the region which is at the root of the CHT. If the point p is
not in the convex hull 0, the tree is not evaluated further and algorithm returns
with False. If the point is in the convex hull of region, then evaluate the children
trees with PICP from left to right recursively. Pruning is done by filtering those
Chapter 5 - AHD BASED GEOMETRIC OPERATIONS 85
sub-trees whose roots do not contain the point p i.e. test PICP for point p and
the root of the sub-tree is FALSE. Any child tree whose root does not contain
p is pruned and not evaluated further. The pruning greatly reduces the search
space thus adding to search efficiency. The result True or False depends on the
level of node where the search stopped. If odd level, search returns False else it
returns true.
For the example in Figure 5.18, the test PICP for p and convex hull 0 returns
True. Therefore, children of hull 0 are evaluated from left to right. Pruning
eliminates all children trees except two children trees whose roots are hull 4 and
hull 7. These trees are then recursively evaluated until the point is found in hull
5. Since the node containing hull 5 is at even level (level 2), True is returned by
child tree whose root is hull 4.
5.7 Operations Summary
Boolean operations are of fundamental importance for GIS processing. The chap-
ter shows that AHD allows optimal (or near optimal) computation of Boolean and
other geometric operations (point in a polygon and line intersecting a polygon)
using AHD approach.
The intersection operation is of fundamental importance as other Boolean
operations union, difference and symmetric difference can be reduced to the com-
putation based on intersection and complement. The complement of a polygon
is computed in constant time by adding or removing a universal node as the root
of the CHT representing the polygon. The union, difference and symmetric dif-
ference are then computed using the principle of duality in terms of intersection
and complement operations.
The algorithms for the computation of these Boolean operation, point in poly-
gon test and line-region intersection are provided. The algorithms are efficient
because the LoD structure of AHD offers pruning and children nodes are processed
only if they qualify some test. LoD together with lazy evaluation (provided by
some programming languages e.g. Haskell) make the geometric algorithms by
reducing the number of computations.
Chapter 6
AHD AND OPERATIONS:
IMPLEMENTATION IN
HASKELL
In previous chapters, I provided details of AHD and associated operations. In
this chapter I will provide implementation of AHD and related issues. All the
algorithms mentioned in previous chapter have been implemented using Haskell
[Jones and Jones, 2003]. The selection of Haskell for implementation provides ease
of implementation giving simple and compact code. Haskell is selected because
it supports (see appendix... for Haskell overview);
� Lazy evaluation: Lazy evaluation means an expression, a function or an
argument passed to a function, or a data structure is evaluated only when its
value is actually used. Lazy evaluation increases performance by avoiding
unnecessary calculations. In case of AHD, lazy evaluation is useful because
the CHT is not built or evaluated unnecessarily until its result is actually
used. The CHT is built or evaluated only to the LoD needed as opposed to
eagerly building or evaluating the CHT entirely.
LoD data structures and lazy evaluation work well together: levels of detail
not necessary to determine the result are not explored [Samet, 2006]. This
effect is automatic in languages with a lazy execution strategy, e.g., Haskell
[Hudak et al., 1992]; it is then possible to chain multiple operations on LoD
structures together and the whole expression is evaluated only to the degree
necessary. In ordinary languages, the expressions must be folded into the
87
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL88
ultimate evaluation code (gradual for each level) to achieve the same effect,
namely avoiding building intermediate results to a degree not necessary.
� High order functions - Functors: Functors are higher order functions
that take another function (s) as an argument. Functors are a principled
way to extend a given algebra in one domain to another domain e.g. Karim-
ipour et al. [2008] has shown lifting functors to lift a domain (e.g., static 2D)
to another domain (e.g., moving 2D). In case of AHD, lifting functors are
needed to achieve dimension independence by extending the AHD algebra
for 2D to n-dimensions.
� Big Integers: Haskell supports arbitrary precision arithmetic. The built
in Integer (big integer) data type implements arbitrary precision arithmetic.
In case of AHD, the point coordinates are represented as homogeneous big
integers to ensure robustness in geometric computations.
� Compact code: Haskell code is concise and compact. Schrage et al. [2005]
have shown that a system developed in Java when reproduced with Haskell
resulted in reduction of lines of code from 220, 000 to 10,000, a reduction
ratio of 22 : 1. Less code leads to less errors, ease of management and
maintenance.
6.1 Input Representation
In this section the representation of input domain region is discussed. This is
important because the AHD algorithm design depends on the representation of
the region.
A region is a collection of polygons which in turn are defined by set of rings.
The outer ring (having counter clockwise orientation) makes the boundary of
the polygon while inner rings make the holes of the polygon and have opposite
orientation than outer boundary. A ring is defined by set of points. Figure shows
the hierarchy of geometric data types needed for the representation of an input
region.
Figure 6.2 shows a region represented with two polygons, polygon1 and poly-
gon2. Polygon1 (Figure 6.2a) has two rings, the outer ring with 12 points and
counter clockwise orientation represents the boundary of the polygon while an
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL89
Figure 6.1: Hierarchy of geometric data types for region representation
(a) Polygon1 (b) Polygon2
Figure 6.2: A Region R represented by two polygons
inner ring with 6 points and clockwise orientation represents a hole. Polygon2
(Figure 6.2b) has only one ring with 8 points and counter clockwise orientation.
The input data types in Haskell are shown below.
type Point = (Int , Int)
data Ring = Ring {_lpoints :: [Point]}
deriving (Show , Eq, Ord)
unRing (Ring ring) = ring
data Poly = Poly {_lrings :: [Ring]}
deriving (Show , Eq, Ord)
unPoly (Poly poly) = poly
data Reg = Reg {_reg :: [Poly]}
deriving (Show , Eq, Ord)
unReg (Reg reg) = reg
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL90
6.2 Implementation of AHD
AHD is based on the convex decomposition of nonconvex polytopes into their
convex components, convex polytopes. The implementation details of the AHD
algorithms for 2D case are provided and shown how AHD is implemented dimen-
sion independently using simplicial decompositions.
6.2.1 2D implementation of AHD
The implementation of AHD needs the definition of geometric data types that
AHD must support. I defined the geometric data types for AHD based represen-
tation model using the algebraic data types provided by Haskell.
6.2.1.1 Geometric data types
For robustness in geometric computations using AHD the point coordinates are
defined by homogeneous big integer (Haskell Integer type) coordinates. It is
named as BPoint to differentiate it from Haskell built in Point type). AHD uses
convex hull as the geometric building block and is defined as the list of BPoint.
The definition of aforementioned data types using Haskell’s algebraic data
types is shown below;
data BPoint = Bpt {_w , _x , _y :: Integer}
deriving (Show , Eq, Ord)
unPoint (Bpt w x y) = (w, x, y)
data CHull = CHull {_chull :: [BPoint ]}
deriving (Show , Eq, Ord)
unChull (CHull ch) = ch
6.2.1.2 Geometric data structure
AHD is based on the convex decomposition of a given nonconvex region into
its convex components. The component convex hulls are stored in a tree data
structure, which is called convex hull tree. The code below shows the definition
of convex hull tree data structure in Haskell.
data Tree p = Node {_this :: p, _childtrees :: [Tree p],
_level :: Level}
| LeafNode {_this :: p, _level :: Level}
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL91
| UniNode { _this :: p, _childtrees :: [Tree p],
_level :: Level}
| NullTree {}
deriving (Eq, Show , Read)
where
p is the geometric primitive stored by nodes of the Tree which is a convex hull
in case of AHD
Node is the constructor that creates a node and every node has three type
variables
this of a node is the variable of type p and keeps the primitive that
node holds
childtrees of a node is the variable of type list of Tree p and keeps the
children trees of that node
level of a node is the variable of type Level
The Level type is useful for level of detail and is defined as
newtype Level = Level Int
deriving (Eq, Ord , Show , Read , Num)
unLevel (Level i) = i
Some of the useful operations on Level are shown below. The function isEven
is useful for PIP test, where a point contained by a node of CHT at even level
means the point in the region represented by the CHT otherwise it is not included.
The inc and dec functions are used for updating the levels of a CHT and are
needed for grafting and complement operations.
isEven (Level l) = even l
inc l = l+1
dec l = l-1
Figure 6.3 shows relationship hierarchy between input representation and
AHD.
6.2.1.3 AHD Algebra
Haskell’s type classes are used to define AHD algebra which will provide the
prototypes of operations or functions of AHD. The two functions of AHD are build
and eval. The build function builds the tree structure by decomposing the given
nonconvex polygon into its convex components and populates the data structure.
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL92
Figure 6.3: Relationship between input representation and AHD
On the other hand, the eval function retrieves the polygon by evaluating the
given tree structure.
class AHD inp rep prim where
build :: Level -> Level -> inp -> [rep prim]
eval :: Level -> Level -> [rep prim]-> inp
The AHD type class defines two functions build and eval for three type vari-
ables inp, rep and prim. The type variable inp represents the type of input which
may be a ring, a polygon or a region. The “rep prim” type variables define the
type of data structure and type of values in that data structure respectively. In
AHD, type variable rep represents the Tree data structure and prim type variable
represents the type of the value in the node of the tree which is the convex hull,
CHull.
The two Level arguments in the build and eval functions are needed for im-
plementing level of detail. The first argument of type Level specifies the level up
to which the tree is to be built or evaluated and the second argument of type
Level is useful for performing recursion and also specifies the starting level num-
ber (which is 0 by default meaning root of any tree has level 0) for build and eval
functions.
The implementation of AHD class relies on functions of the Input and InpPrim
classes which are shown below.
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL93
Table 6.1: Functions of Input and Prim classes
Function Description Example
merge Given two overlapping polygons (e.g. p1 andp2), the function merges the two and returnsthe resulting new polygon (e.g. p)
split Given a polygon (e.g. p), the function splitsthe polygon and returns its extracted deltasor nonconvexties (e.g. d1 and d2 ).
makePrim Given a polygon (e.g. p), the function makesthe its primitive which is the convex hull ofthe object (e.g. ch)
mergPrim Given two convex hulls (one parent and onechild e.g. ch0 and ch1 respectively), thefunction merges the two hulls and returnsthe resulting polygon ring r as shown in theexample figure on right.
class Input inp where
merge :: inp -> inp -> inp
split :: inp -> [inp]
class InpPrim inp prim where
makePrim :: inp -> prim
mergePrim :: prim -> prim -> inp
Table 6.1 shows a brief description of the functions in both classes for the case
of an input polygon (inp) and convex hull (prim) with some examples.
6.2.1.4 Instances for AHD algebra
The algebra defined above provides a framework for generic implementation of
build and eval functions of AHD. This is achieved by instantiating the defined
classes for different input geometric types. For example, the code snippet below
shows the instantiation of the build and eval functions for a Ring, Poly and Reg
data types.
instance AHD Ring Tree CHull where
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL94
build mlev clev ring = if (drings == []) || (mlev ==
clev) then [( LeafNode p clev)]
else [Node p (concat (map (build mlev nl)
drings)) clev]
where
p = makePrim ring
drings = split ring
nl = inc clev
eval mlev [( LeafNode (CHull p) l)] = Ring p
eval mlev [(Node (CHull p) lts l)] = if (mlev == l)
then Ring p
else (foldl mergeFunc (Ring p) llr)
where
llr = (map (eval mlev) (toLL
lts))
instance AHD Poly Tree CHull where
build mlev clev (Poly []) = []
build mlev clev (Poly lr) = [foldl mergTrees (head
ltrees) (tail ltrees)]
where
ltrees = concat (map (build mlev clev
) lr)
eval mlev [( LeafNode (CHull p) l)] = Poly [Ring p]
eval mlev [(Node (CHull p) lts l)] = if (mlev == l)
then Poly [Ring p]
else (foldl mergeFunc (Poly [Ring p]) lp)
where
lp = (map (eval mlev) (toLL
lts))
instance AHD Reg Tree CHull where
build mlev clev (Reg []) = []
build mlev clev (Reg reg) = ltrees
where
ltrees = concat (map (build mlev clev
) reg)
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL95
eval mlev lts = Reg lp
where
lp = (map (eval mlev) (toLL lts))
The implementations of the functions in the the two classes, Input and InpPrim,
are shown in the code snippet below.
instance Input Ring where
merge (Ring r0) (Ring r1) = Ring (edge2vertex nmre)
where
er0 = mkEdges r0
er1 = mkEdges r1
oes = [(ce0 , ce1)|ce0 <- er0 , ce1 <-
er1 , isEdgeOn ce1 ce0]
re0 = er0\\ (map fst oes)
de1 = er1 \\ (map snd oes)
re1 = map revOrient de1
nmes = concat(map mergEdges oes)
nmre = re1 ++ re0 ++ nmes
split ring = map listedges2Ring dre
where
dre = sepReg fde (CHull (nub( (
unChull ch) ++ dph)))
dph = nub.flatten $ (de \\ fde)
fde = filter (\y -> fst(isEdgeOnReg y
he)==False) de
de = (mkEdges.unRing $ ring) \\ he
he = mkEdges.unChull $ ch
ch = quickHull ring
instance Input Poly where
merge (Poly p0) (Poly p1) = Poly (mergeRingLists p0
p1)
instance InpPrim Ring CHull where
makePrim ring = quickHull ring
mergePrim (CHull r0) (CHull r1) = [Ring (edge2vertex
nmre)]
where
er0 = mkEdges r0
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL96
er1 = mkEdges r1
oes = [(ce0 , ce1)|ce0 <- er0 , ce1 <-
er1 , isEdgeOn ce1 ce0]
re0 = er0\\ (map fst oes)
de1 = er1 \\ (map snd oes)
re1 = map revOrient de1
nmes = concat(map mergEdges oes)
nmre = re1 ++ re0 ++ nmes
instance InpPrim Poly CHull where
makePrim poly = quickHull (head.unPoly $ poly)
6.2.2 Dimension independent implementation of AHD
AHD can be implemented dimension independently by lifting operations [Karim-
ipour et al., 2008]. Table 6.2 lists some of the operations on the simplexes, which
are used for developing the decomposition algorithm.
For dimension independent implementation of AHD, the functions of Input
and InputPrim need dimension independent implementation. The code below
shows the definition of geometric data types and the instances for the classes of
AHD algebra.
-- An n dimnesional point defined as a list of homogenous big
integer points
type ND_BPoint = [BPoint]
-- An n- simplex is defined as a list of n dimensional
vertexes
type Simplex = [ND_Vertex]
-- Canonical representation of a n-simplex
type CnSimplex = (Simplex , Bool)
instance AHD [CnSimplex] Tree [CnSimplex] where
build mlev clev reg = if (dreg == []) || (mlev ==
clev) then [LeafNode chReg clev]
else [Node chReg (concat (map (build mlev (
inc clev)) dreg)) clev]
where
chReg = makePrim reg
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL97
Table 6.2: Operations on simplexes
Operation Description Example
addVertex(v,S) Adding a vertex v to an
n-simplex S , which pro-
duces an (n+1)-simplex
S1 =< v0,v1 > addVertex(v,S1) =< v0,v1 >
cw(p,S) Testing the position of a
point p w.r.t. a simplex
S (clockwise or counter-
clockwise) S1 =< v0,v1 >
cw(p,S1) = f alseccw(p,S1) = true
ptInSimplex(p,{Si}) Testing if a point p is in-
side the region delimited
by a set of simplexes {Si}
R = {S11,S
51}
ptInSimplex(p1,R) = trueptInSimplex(p2,R) = f alse
Border({Si}) Extracting the border of
a set of connected sim-
plexes {Si}
R = {Si}({Si}=all edges)
Border(R) = {S j}({S j}=bold edges)
Split({Si}) Splitting a set of
simplexes{Si} to the
connected subsets
{R1, ...,Rk}
Split({S11, ...,S
91}) =
({S11, ...,S
51},{S6
1, ...,S81},{S9
1})
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL98
dreg = split $ symDiff reg
chReg
eval mlev [( LeafNode s l)] = s
eval mlev [(Node s lts l)] = if (mlev == l) then s
else symDiff s rs
where
rs = concat(map (eval mlev) [
lts])
instance InputPrim [CnSimplex] [CnSimplex] where
makePrim ls = convexHull (nub.concat.map cnSimp2simp
$ ls)
instance Input [CnSimplex] where
split = splitJoints.splitOrdRegs
6.3 Implementation of AHD based Geometric Oper-
ations
6.3.1 AHD operations algebra
The algebra of AHD operations is shown below. The class Prim defines the ccint
(convex-convex intersection) between two convex geometric primitives and the
result is also convex. The class will be instantiated for the AHD primitive convex
hull as the value of type variable prim, the intersection of two convex hulls is
always convex.
The class AHDOperations defines the functions to be implemented for AHD.
class Prim prim where
ccint :: prim -> prim -> prim
class AHDOperations rep prim where
cnint :: prim -> rep prim -> [rep prim]
gint :: rep prim -> rep prim -> [rep prim]
comp :: [rep prim]-> [rep prim]
uni :: rep prim -> rep prim -> [rep prim]
diff :: rep prim -> rep prim -> rep prim
symDiff :: rep prim -> rep prim -> rep prim
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL99
pir :: BPoint -> rep prim -> Bool
6.3.2 Implementation of intersection operation
As discussed earlier, intersection computation is the most important operation
as other operations like union difference and symmetric difference are computed
using the intersection operation. The code for an instance of the Prim class that
implements the ccint function for a convex hull primitive is shown below.
instance Prim CHull where
ccint (CHull p1) (CHull p2) = if (p1\\pinp2 == [])
then CHull p1
else if (p2\\pinp1 == []) then CHull p2
else if (( length.unChull $ ip) <=2)
then CHull []
else ip
where
(CHull chp) = quickHull (Ring (union
p1 p2))
che = mkEdges chp
ine1 = (mkEdges p1) \\ che
ine2 = (mkEdges p2) \\ che
vip = fip $ pie ine1 ine2
pinp2 = pointsInOrOn (CHull p1) (
CHull p2)
pinp1 = pointsInOrOn (CHull p2) (
CHull p1)
ip = quickHull (Ring (pinp1 ++ vip ++
pinp2))
The intersection of a convex and a nonconvex object can be convex or noncon-
vex. Using the ccint function, the intersection of a convex hull and a nonconvex
object is computed by recursively applying ccint over the tree of nonconvex ob-
ject. The generic intersection computation function gint is then implemented
using cnint that computes intersection between two trees which may represent a
convex or nonconvex object as shown in the code below.
instance AHDOperations Tree CHull where
cnint ch (LeafNode x l) = [LeafNode (ccint ch x) l]
cnint ch (Node x xts l) = if (ohi == CHull [] ) then
[]
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL100
else [( Node ohi (concat(map (cnint ohi) xts)
) l)]
where
ohi = (ccint ch x)
gint t (LeafNode y ly) = cnint y t
gint (LeafNode x lx) t = cnint x t
gint (Node x lct lx) t = if (oi == []) then []
else [foldl mergTrees (head oi) lt]
where
oi = cnint x t
lt = concat (map (gint (head
oi)) lct)
intReg r1 r2 = concat [(gint p x)| p<- r1, x<- r2]
...
6.3.3 Implementation of complement operation
Using the principle of duality, the union and symmetric difference operations are
computed in terms of intersection and complement operations (De Morgans Law).
The implementation of complement operation is shown below.
instance AHDOperations Tree CHull where
...
compReg lpt = if (length lpt ==1) && (is_Uni.head $
lpt) then dts
else [UniNode (uniHull) its 0]
where
its = map (incTLevel 0) lpt
dts = map (decTLevel 0) (tail
lpt)
...
The function compReg complements the given region by adding all the polygon
trees of the region as the children of the universal node, UniNode. The universal
node contains the universal hull. For simplification the drawing space could be
assumed as the finite universal space, in which case , the convex hull of the
drawing space becomes the universal hull. According to the principle of duality,
the dual of a dual function is identity. The compReg implements this property,
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL101
and the complement of a complemented region returns the original region.
6.3.4 Implementation of union operation
The union operation is implemented as a dual function using intersection and
complement functions (r1∪ r2 = r1∩ r2 ) , the implementations of which have
already been given.
The function uniReg computes the union of two regions. The complement of
the intersection of complemented regions gives the union region.
instance AHDOperations Tree CHull where
...
uniReg r1 r2 = compReg (intReg cr1 cr2)
where
cr1 = compReg r1
cr2 = compReg r2
...
6.3.5 Implementation of difference and symmetric difference op-
erations
The implementation of symmetric difference needs the implementation of differ-
ence operation. The difference of two regions is implemented using the comple-
ment and intersection operations.
r1\ r2 = r1∩ r2
r2\ r1 = r2∩ r1
The symmetric difference of two regions is then computed by taking the union
of two difference regions r1\ r2 and r2\ r1.
r14 r2 = (r1\ r2)∪ (r2\ r1)
instance AHDOperations Tree CHull where
...
difReg r1 r2 = intReg r1 (compReg r2)
symDReg r1 r2 = uniReg r1minr2 r2minr1
where
r1minr2 = difReg r1 r2
r2minr1 = difReg r2 r1
...
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL102
6.3.6 Implementation of point in polygon
The point in polygon uses the geometric function isPointinHull which determines
whether a given point is in the hull or not. The children of a tree are evaluated
only if the given isPointinHull returns true for the given point and the convex
hull at the root. This is lazily done because the Haskell’s built in lazy evaluation,
thus pruning children trees whose root hulls do not include the point.
instance AHDOperations Tree CHull where
...
pip p (LeafNode ch l) = (isPointinHull p ch) && (
isEven l)
pip p t = (isPointinHull p ch) && (all (== False) lb)
&& (isEven l)
where
ch = this t
l = level t
cts = childtrees t
lb = map (pip p) cts
pir p r = (any (== True) lpip)
where
lpip = map (pip p) r
6.4 Implementation of Generalized LoD Data Struc-
ture
In the previous sections I showed that Tree data structure of AHD supports LoD
and design of progressive algorithms for performing Boolean operations. In this
section I will show that the AHD’s LoD based tree data structure works for other
LoD based data structures and operations. This is shown with two examples,
first in which the data is divided over LoD by delta is the segment tree, and the
other in which the data is divided over the LoD by replacement is the strip tree
[Ballard, 1981].
1. Segment tree - Intersection computation between a region represented in
the CHT of AHD and a line segment is done progressively by intersecting
the line segment with convex hulls stored in the CHT. The intersection
between a line and a convex hull (assuming that both intersect) is always
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL103
Figure 6.4: The function f(x) with three approximations
Figure 6.5: One piece of a linear function approximating the given function
a line segment, which is stored in a tree every node of which contains a list
of line segments, the tree called as segment tree.
2. Strip tree - Strip trees approximate lines with straight line segments describ-
ing for each level the error bounds on the approximation. The approximated
line segments are called a strip organized in a tree called a strip tree. This
can be seen as an approximation of a continuous function in one variable
y = f (x) by a hierarchy of piecewise linear functions (Fig. 6.4). The linear
function in an interval (Fig. 6.4) is gradually approximated with linear
functions, for smaller and smaller intervals Fig. 6.5).
6.4.1 Implementation of segment-region intersection
The implementation of algorithm for computing the intersection between a line
segment and a region is shown below. The function intLinReg computes the
intersection between a given line segment (which is stored in data type LinSeg)
and a region which is achieved by mapping the function intLinPoly over the list
of trees of the polygons of the region.
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL104
data LinSeg = LinSeg {_sp , _ep :: BPoint}
deriving (Show , Eq, Ord)
linSeg (LinSeg p1 p2) = (p1, p2)
intLinReg :: Level -> Level ->LinSeg -> [Tree CHull] -> [Tree
[LinSeg ]]
intLinReg ml cl s rts = concat (map (intLinPoly ml cl s) rts)
intLinPoly :: Level -> Level -> LinSeg -> Tree CHull -> [Tree [
LinSeg ]]
intLinPoly ml cl s (LeafNode ch l) = [LeafNode (intLinHull s
ch) l]
intLinPoly ml cl s (Node ch cts l) = if (int == []) then []
else if (ml == cl) then [LeafNode (intLinHull s ch)
cl]
else [Node int cs cl]
where
int = intLinHull s ch
cs = concat (map (intLinPoly
ml (cl+1) s) cts)
The result is a list segment trees, each node of which stores a line segment.
Each segment is then evaluated, resulting in a segment or list of segments which
is achieved by line segment (primitive) merging. The instance of merge for LinSeg
is is shown below which is used by the eval function of the AHD.
instance Input [LinSeg] where
merge ls1 ls2 = ns
where
os = [(cs1 , cs2)|cs1 <- ls1 , cs2 <-
ls2 , isEdgeOn (s2e cs1) (s2e cs2)
]
rs1 = ls1\\ (map fst os)
rs2 = ls2 \\ (map snd os)
os2e = [(s2e s1, s2e s2)| (s1,s2) <-
os]
ms = concat(map mergEdges os2e)
oe2s = [e2s e| e <- ms]
ns = rs1 ++ rs2 ++ oe2s
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL105
6.4.2 Implementation of strip tree
Two data types are defined for the implementation of strip tree algorithm. The
input is a curve defined by a list of point, so the type Curve is a synonym to the
list of points type. The primitive is a strip which is defined as the data type Strip
having six elements where spx and spy denote the beginning of the strip, epx and
epy denote the end of the strip and wl and wr are the left and right distances of
the strip borders.
type Curve = [BPoint]
data Strip = Strip { spx :: Double , spy :: Double , epx ::
Double , epy :: Double , wl :: Double , wr :: Double}
deriving (Show , Read)
The instantiation of the AHD algebra for strip tree implementation is shown
below.
instance AHD Curve Tree Strip where
build mlev clev c = if (mlev == clev) || (
is_Splittable c == False) then [( LeafNode p clev)
]
else [Node p (concat (map (build
mlev nl) parts)) clev]
where
p = makePrim c
parts = split c
nl = clev+1
eval mlev [( LeafNode
s l)] = getCurve
s
eval mlev [(Node s lts l)] = if (mlev == l) then (
getCurve s)
else (merge (getCurve s) (concat llr))
where
llr = (map (eval mlev) (toLL
lts))
instance Input Curve where
is_Splittable curve = length curve > 2
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL106
split c = parts
where
iline = getline.makePrim $ c
eps = getxtrmpoints iline c
breakpoint = fst (maximumBy (compare
‘on‘ snd) eps)
parts = splitList breakpoint c
merge _ lc = nub $ lc
instance InputPrim Curve Strip where
makePrim (p1: p2: []) = Strip { spx = p1x , spy =
p1y , epx = p2x , epy = p2y , wl = 0.0, wr = 0.0}
where
(1.0, p1x , p1y) = h2c (unPoint p1)
(1.0, p2x , p2y) = h2c (unPoint p2)
makePrim lp = Strip { spx = p1x , spy = p1y , epx = p2x
, epy = p2y , wl = wll , wr = wrr}
where
(1.0, p1x , p1y) = h2c (unPoint p1)
(1.0, p2x , p2y) = h2c (unPoint p2)
(p1, p2) = getinitline lp
wrr = snd.head $ lxt
wll = snd.last $ lxt
lxt = getxtrmpoints (p1, p2) lp
The build function works for an input of type Curve and results in an output
of type [Tree Strip] i.e. a list of strip trees that contains only one strip tree of the
input curve. Similarly the eval function works when the argument is of type [Tree
Strip] and returns the output of the type Curve i.e. the strip tree is evaluated up
to the level specified and returns the corresponding curve. This is shown with an
example curve, expcurve, which is a list of points in the code in the code below.
--example curve of type Curve
expcurve = [(Bpt 1 3 7), (Bpt 1 9 12), (Bpt 1 15 4), (Bpt 1
18 5) ,(Bpt 1 20 7)] :: Curve
-- build strip tree of expcurve up to 3 levels , the root of
strip tree should have level 0
bst = build 3 0 expcurve :: [Tree Strip]
Chapter 6 - AHD AND OPERATIONS: IMPLEMENTATION IN HASKELL107
-- evaluate bst , the strip tree of expcurve up to the level 2
est = eval 2 bst :: Curve
Chapter 7
CONCLUSION AND
FUTURE WORK
7.1 Conclusion
The research showed that a systematic exploration is possible for the represen-
tation of arbitrary objects (convex or nonconvex, with or without holes) by de-
coupling metric (numeric) and geometry considerations. The proposed solution
separates the metric and geometric representation and uses the combination of
big integers homogeneous coordinates for metric and convex polytopes for topo-
logical processing. The literature revealed that combination of big integers and
convex polytopes was unexplored and the research showed that the combination
is useful because it provides opportunity to investigate independently the metric
operations and coordinate representations from representation of topology of the
geometric objects. The advantage of metric calculations of big integers with ho-
mogeneous coordinates is that they are accurate for straight line structures, not
approximates and all topological relations derived are consistent. Intersection
points calculated are exactly on the line, not only “very near”, etc. The solution
provides the geometric model AHD that has the following properties:
� ease of implementation
� support for operations
� level of detail support
� extendable to n-Dimensions
109
Chapter 7 - CONCLUSION AND FUTURE WORK 110
The following questions were answered;
7.1.1 Is the use of big integers practical? How slow are big inte-
gers?
The question was investigated thoroughly in section 4.1. The experiment results
showed that metric computations using homogeneous coordinates with big inte-
gers were conceptually simpler and easy to program. Also, the performance was
acceptable for GIS geometry. The experiment results showed that;
� Solutions using big numbers are realistic for GIS. For applications including
cascaded construction of geometric objects, the growth in the size of the
big numbers representing the coordinates may be problematic. The average
size doubles with each cascaded computation. Typical GIS operations do
not construct long chains of intersections of lines defined as intersections of
other lines.
� Contrary to the widespread belief, programming with big numbers was
straightforward and easier than writing code using conventional Integer or
Float data types. I showed that the performance penalty is relatively small
and affects only the part of GIS operations computing metric relations for
cascaded construction of points geometrically.
As a summary, metric computations with big numbers lead to straightforward
programs that are guaranteed to be consistent. The performance penalty seems
bearable and the memory utilization can be controlled.
7.1.2 Is the combination of convex polytope and big integers
attractive?
I showed that the novel combination of big integers with convex polytopes is prac-
tical for the representation of GIS geometry that is guaranteed to be consistent.
The solution advocated in 1986, namely, floating point arithmetic with simpli-
cial complex [Frank and Kuhn, 1986] can be characterized as “eager”, leading to
too many computations that are not necessarily meaningful. A similar approach
[Thompson, 2007]uses regular polytopes (defined as union of convex polytopes),
but half spaces are defined by finite precision representations. The proposed ap-
proach used convex polytopes in the dual space in which the coordinates of the
point are homogeneous big integers.
Chapter 7 - CONCLUSION AND FUTURE WORK 111
7.1.3 Does the model support Boolean operations?
The model offers support for Boolean operations like intersection, union and sym-
metric difference etc. The convex decomposition and the hierarchical tree data
structure facilitate the computations necessary for performing operations. Inter-
section is the only closed operation (as the intersection of two convex polytopes
is always convex) and most important operation (other operations depend on it).
The solution provided an algorithm for the computation of the intersection of
two general polygons polygons (convex or nonconvex and with or without holes).
The union of two convex polytopes is may be nonconvex. The principle of duality
was to tackle the issue, taking the inverse of given polygons and intersecting them
and then taking the inverse of result provided the union of the two polygons (i.e.
A ∪ B = (A ∩ B)).
The other operations, symmetric difference, point in a polygon and line in-
tersections with a given point were then computed using the intersection and
complement operations.
Once the questions regarding the basic requirements were answered, I explored
additional requirements the model fulfills;
7.1.4 Support for level of detail
Level of detail is an important concept in geometric and graphical description
of objects. The hierarchical tree data structure of the proposed solution, AHD,
provides support for geometric level of detail. I have argued that the use of
convex polytopes as geometric primitives and the hierarchical data structure lead
to algorithms which can be implemented as progressive algorithms. The LoD in
AHD facilitates reduced number of computations.
I concluded that LoD data structures and progressive evaluation are necessary
for GIS and other information systems, where data are stored and requested from
applications with very differing quality requirements. Progressive LoD evaluation
can be achieved generally in programming languages with a lazy evaluation strat-
egy, with the programmer only avoiding depth first recursive calls; this shows that
the ’progressive’ concept can be separated from and dealt with independently of
details of LoD. The analysis showed a further step of generalization, as it sep-
arates generic LoD data structure and operations from operations on primitives
used to construct the approximations.
Chapter 7 - CONCLUSION AND FUTURE WORK 112
Table 7.1: Description of the three test data sets
DatasetNo
Data Type No. ofpoly-gons
No. ofrings
Totalno ofpoints
1 buildingfootprints
5060 5093 49003
2 cadastre 1238 1238 15549
3 cadastre (datasource [Kui Liu
et al., 2007])
972 974 9871
7.1.5 Dimension independence
Most of the current decomposition algorithms are dimension specific needing di-
mension specific implementations. The proposed approach AHD is dimension
independent that is achieveable by lifting the AHD algebra for 2D to nD. For
dimension independent AHD, the geometric primitive convex polytope is imple-
mented as a n- convex polytope (convex polytope of any dimension) and the oper-
ations on the primitive are generalized. The most important geometric operation
is the convex hull computation, for which any existing dimension independent
convex hull algorithm can be used.
In order to show the dimension independence, AHD is implemented by decom-
posing a polytope into a simplicial complex. The approach uses the dimension
independent convex hull, called IncConvexHull algorithm and a lifted split func-
tion that extracts the deltas or notches.
7.1.6 Performance analysis
Finally, to verify the proposed solution, I tested AHD with three real datasets
having large number of vertexes. This proved useful for getting an objective
insight into the strengths and weaknesses of the model. The description of the
three datasets used is given in Table 7.1.
In order to measure the intersection computation performance for AHD, the
test dataset regions r1, r2 and and r3 are intersected against randomly drawn
(clip) regions c1, c2 and c3 respectively as shown in Figure 7.1.
The three intersection regions i1 (r1∩ c1), i2 (r2∩ c2) and i3 (r3∩ c3) are
shown in Figure.
Chapter 7 - CONCLUSION AND FUTURE WORK 113
(a) Dataset1 region r1 and clip region c1
(b) Dataset2 region r2 and clip region c2
(c) Dataset3 region r3 and clip region c3
Figure 7.1: Test datasets and clip regions
Chapter 7 - CONCLUSION AND FUTURE WORK 114
(a) Intersection region i1 (r1∩ c1)
(b) Intersection region i2 (r2∩ c2)
(c) Intersection region i3 (r3∩ c3)
Figure 7.2: Intersection regions
Chapter 7 - CONCLUSION AND FUTURE WORK 115
Table 7.2: Region statistics where r1, r2, and r3 are input regions: c1, c2 and c3are the clip regions intersected with regions r1, r2, and r3 respectively: i1=(r1∩c1), i2=(r2∩ c2) and i3 =(r3∩ c3)
Region No. ofPoints
No. of trees Max.no. ofnodes
Max.no. oflevels
r1 49003 5060 48 4
r2 15549 1238 135 6
r3 9871 972 89 6
c1 8 1 5 2
c2 6 1 2 1
c3 15 2 4 2
i1 3462 1777 76 5
i2 663 249 42 5
i3 1679 442 104 6
Table 7.3: Time for AHD and intersection computation
Data-
set
Buildtime(bt)
bt/point Evaltime(et)
et/point Int.time(it)
it/point
1 2.58s 52.67µs 1.83s 37.26µs 8.37s 2417.48µs2 1.82s 116.95µs 1.11s 71.71µs 0.48s 717.92 µs3 0.93s 94.48µs 0.68s 69.33µs 4.22s 2516.52µs
Table 7.2 summarizes the region statistics for the test dataset, the clip and
the intersection regions.
Table1 7.3 shows the time elapsed in seconds for build and eval functions of
AHD and the computation of intersection regions using a system having a dual
core processor of 2 x 1 GHz and 2GiB RAM.
From the results it is deduced that the performance for build, eval and inter-
section depends on the number of nodes and levels in the region trees or in other
words on the shape of the region polygons. This is shown in Figure 7.3 where the
build and eval times per point for the regions r1, r2 and r3 are plotted together
with the maximum number of nodes in a tree.
The same relationship is shown for time per point for intersection and the
maximum number of nodes in a tree as shown in Figure 7.4.
1All values are rounded to two decimal places for better presentation
Chapter 7 - CONCLUSION AND FUTURE WORK 116
Figure 7.3: Relationship between time per point (in µs) and max. no. of nodesfor AHD operations build and eval
Figure 7.4: Relationship between time per point (in µs) on logarithmic scale andmax. no. of nodes for AHD intersection operation
Chapter 7 - CONCLUSION AND FUTURE WORK 117
7.2 Future Work
7.2.1 Support for non simple polygons
AHD assumes the polygons are simple (without self intersecting edges) which
could be problematic for some real datasets. The possible solution to the prob-
lem is introducing a preprocessing step in which a non simple polygon is splitd
into two or more polygons depending on the shape. This is done by computa-
tion the intersection points of intersecting edges, updating the intersecting edges
and then ultimately segregating the closed regions formed, so as to represent it
with multiple simple polygons. The idea is to convert or decompose a non simple
polygon with multiple simple polygons or use any established method to deal
non simplicity first and then continue with proposed solution. Since AHD sup-
ports multiple polygons, the preprocessing will overcome the problem and would
increase the usability and applicability of AHD for practical GIS.
7.2.2 Implemention of AHD in a database
In some large scale applications it may be more efficient to have a separate geo-
metric database to store large number of geometric objects [Gunther, 1988a]. In
future, the AHD will be implemented in three different database environments
relational (e.g. mySQL), object oriented (e.g. Oracle) and document-oriented
(couchDB2). This will need support of indices that are used to optimize spatial
queries. Popular spatial indices are quadtree [Samet, 1984], R-tree [Guttman,
1984], R+tree [Sellis et al., 1987], cell tree[Gunther, 1988a] and B-tree [Bayer
and McCreight, 1970] etc. The index should be implemented as a secondary
storage data structure that facilitates efficient search queries.
7.2.3 AHD based geometric modeling of buildings
A 3D city model may consist of many different objects needed to be modeled
depending on the application domain. Some of the most important city objects
are buildings, roads, streets, trees, etc. However, buildings are the most important
part of a city model for many applications (Lancelle and Fellner, 2004).
The approach for using AHD for the geometric modeling of buildings involves;
the computation of the convex hull of the building 3D points which makes the
root of the CHT. The concavities are extracted and the process is repeated for
2http://couchdb.apache.org/
Chapter 7 - CONCLUSION AND FUTURE WORK 118
Figure 7.5: AHD based geometric modeling of buildings
each of the extracted nonconvex part, which are ultimately used to approximate
extruding structures like chimneys, or surface structures like doors windows etc
or sometimes just “noise space”. By noise space, I mean the space which is not
part of the actual object space, but is part of the “approximation space”.
Figure 7.5 shows an example for geometric modeling of buildings (in 2D) using
AHD. For dem
newtype Level = Level Int
deriving (Eq, Ord , Show , Read , Num)
unLevel (Level i) = i
onstration purposes the input is a 2D point set representing a building in 2D (Fig-
ure 7.5a). The decomposed convex components and the corresponding CHT are
shown in Figures 7.5b and 7.5c respectively. The convex hulls 1 and 2 represent
noise space of the overall approximation space represented by the CHT.
The application of AHD for geometric modeling of city objects seems attrac-
tive as an additional geometric model, however it needs further investigation.
Some of the issues for future investigation are;
� Integration of ontological information into AHD
� Relating geometric and semantic levels of detail
Bibliography
T. Ai, Z. Li, and Y. Liu. Progressive transmission of vector data based on changes
accumulation model. Developments in Spatial Data Handling, pages 85–96,
2005.
P. Alexandroff and PS Aleksandrov. Elementary concepts of topology. Dover
Publications, 1961.
O. El Badawy and M.S. Kamel. Hierarchical representation of 2-d shapes using
convex polygons: a contour-based approach. Pattern Recognition Letters, 26
(7):865 – 877, 2005. ISSN 0167-8655.
C.L. Bajaj and T.K. Dey. Convex decomposition of polyhedra and robustness.
SIAM J. Comput., 21(2):339–364, 1992.
Dana H. Ballard. Strip trees: A hierarchical representation for curves. ACM
Comm., 24(5 (May)):310–321, 1981.
C. Bradford Barber, David P. Dobkin, and Hannu Huhdanpaa. The quickhull
algorithm for convex hulls. ACM Trans. Math. Softw., 22(4):469–483, 1996.
ISSN 0098-3500.
R. Bayer and E. McCreight. Organization and maintenance of large ordered
indices. In Proceedings of the 1970 ACM SIGFIDET (now SIGMOD) Workshop
on Data Description, Access and Control, SIGFIDET ’70, pages 107–141, New
York, NY, USA, 1970. ACM.
J.L. Bentley and T.A. Ottmann. Algorithms for reporting and counting geometric
intersections. IEEE Transactions on Computers, 28(9), 1979.
Jules Bloomenthal and Jon Rokne. Homogeneous coordinates. The Visual Com-
puter, 11:15–26, 1994. ISSN 0178-2789.
119
Bibliography 120
M. Breunig. On the way to component-based 3D/4D geoinformation systems.
Springer Verlag, 2001.
C. Burnikel, R. Fleischer, K. Mehlhorn, and S. Schirra. Efficient exact geometric
computation made easy. In Proceedings of the fifteenth annual symposium on
Computational geometry, page 350. ACM, 1999.
B. Chazelle. An optimal algorithm for intersecting three-dimensional convex
polyhedra. Foundations of Computer Science, Annual IEEE Symposium on, 0:
586–591, 1989.
B. Chazelle and DP Dobkin. Intersection of convex objects in two and three
dimensions. Journal of the ACM (JACM), 34(1):1–27, 1987.
B. Chazelle and H. Edelsbrunner. An optimal algorithm for intersecting line
segments in the plane. Journal of the ACM (JACM), 39(1):1–54, 1992.
B. Chazelle and L. Palios. Triangulating a nonconvex polytope. Discrete and
computational geometry, 5(1):505–526, 1990.
N.R. Chrisman, J.A. Dougenik, and D. White. Lessons for the design of polygon
overlay processing from the Odyssey Whirlpool algorithm. In Proc. 5th Int.
Symp. on Spatial Data Handling, volume 2, pages 401–410, 1992.
M. De Berg, O. Cheong, M. Van Kreveld, and M. Overmars. Computational
geometry: algorithms and applications. Springer-Verlag New York Inc, 2008.
L. De Floriani, P. Magillo, and E. Puppo. Applications of computational geometry
to geographic information systems. Handbook of computational geometry, pages
333–388, 1999.
M.J. De Smith, M.F. Goodchild, and P. Longley. Geospatial analysis: a compre-
hensive guide to principles, techniques and software tools. Troubador Publish-
ing, second edition, 2007.
Herbert Edelsbrunner and Ernst Peter Mucke. Simulation of simplicity: a tech-
nique to cope with degenerate cases in geometric algorithms. ACM Trans.
Graph., 9(1):66–104, 1990. ISSN 0730-0301.
M. Egenhofer, A. Frank, and J. Jackson. A topological data model for spatial
databases. Design and Implementation of Large Spatial Databases, pages 271–
286, 1990.
Bibliography 121
F. R. Feito and M. Rivero. Geometric modelling based on simplicial chains.
Computers & Graphics, 22(5):611 – 619, 1998. ISSN 0097-8493.
J. Fernandez, L. Canovas, and B. Pelegrin. Algorithms for the decomposition of
a polygon into convex polygons. European Journal of Operational Research,
121(2):330 – 342, 2000. ISSN 0377-2217.
E. Fogel, R. Wein, B. Zukerman, and D. Halperin. 2D Regularized Boolean Set-
Operations. CGAL Editorial Board, CGAL-3.2 user and reference manual edi-
tion, 2006.
S. Fortune. Stable maintenance of point set triangulations in two dimensions.
In Proceedings of the 30th Annual Symposium on Foundations of Computer
Science, pages 494–499. IEEE Computer Society, 1989.
S. Fortune. Polyhedral modelling with exact arithmetic. In Proceedings of the
third ACM symposium on Solid modeling and applications, pages 225–234.
ACM, 1995.
Steven Fortune and Christopher J. Van Wyk. Static analysis yields efficient
exact integer arithmetic for computational geometry. ACM Trans. Graph., 15
(3):223–248, 1996. ISSN 0730-0301.
A. U. Frank and W. Kuhn. Cell graphs: A provable correct method for the
storage of geometry. In D. Marble, editor, Second International Symposium on
Spatial Data Handling, pages 411 – 436, Seattle, Wash., 1986.
Andrew U. Frank. Spatial concepts, geometric data models, and geometric data
structures. Computers & Geosciences, 18(4):409 – 417, 1992. ISSN 0098-3004.
W.R. Franklin. Cartographic errors symptomatic of underlying algebra problems.
In Proc. First Intl. Symposium on Spatial Data Handling, Z
”urich, volume 1, pages 190–208, 20–24 August 1984.
Christine Goldenhuber. Aufdeckung von Numerischen Problemen in geodatischer
Software, volume 11 of Geoinfo Series. Institute for Geoinformation, 1997.
Daniel H. Greene and F. Frances Yao. Finite-resolution computational geometry.
In Foundations of Computer Science, 1986., 27th Annual Symposium on, pages
143 –152, 1986.
Bibliography 122
Gunther Greiner and Kai Hormann. Efficient clipping of arbitrary polygons. ACM
Trans. Graph., 17(2):71–83, 1998. ISSN 0730-0301.
L. Guibas, L. Ramshaw, and J. Stolfi. A kinetic framework for computational
geometry. In Foundations of Computer Science, 1983., 24th Annual Symposium
on, pages 100–111, 1983.
O. Gunther. Efficient Structures for Geometric Data Management, volume 337
of Lecture Notes in Computer Science. Springer Verlag, 1988a.
O. Gunther. A dual approach to detect polyhedral intersections in arbitrary
dimensions. In Oliver Gunther, editor, Efficient Structures for Geometric Data
Management, volume 337 of Lecture Notes in Computer Science, pages 49–64.
Springer Berlin / Heidelberg, 1988b.
O. Gunther and E. Wong. Convex polyhedral chains: A representation for geo-
metric data. Computer-aided Design, Butterworth & Co., 21(3):157–164, 1989.
copy.
Ralf Guting and Markus Schneider. Realms: A foundation for spatial data types
in database systems. In David Abel and Beng Chin Ooi, editors, Advances in
Spatial Databases, volume 692 of Lecture Notes in Computer Science, pages
14–35. Springer Berlin / Heidelberg, 1993.
R.H. Guting and M. Schneider. Realm-based spatial data types: the ROSE
algebra. The VLDB Journal, 4(2):243–286, 1995.
Antonin Guttman. R-trees: A dynamic index structure for spatial searching. In
Proceedings of the 1984 ACM SIGMOD international conference on Manage-
ment of data, SIGMOD ’84, pages 47–57, New York, NY, USA, 1984. ACM.
ISBN 0-89791-128-8.
I. Hanniel and R. Wein. An exact, complete and efficient computation of arrange-
ments of Bezier curves. In Proceedings of the 2007 ACM symposium on Solid
and physical modeling, page 263. ACM, 2007.
John Herring. TIGRIS: A data model for an object-oriented geographic informa-
tion system. Computers and Geosciences, 18:443–452, 1991.
Christoph M. Hoffmann. The problems of accuracy and robustness in geometric
computation. Computer, 22(3):31–40, 1989. ISSN 0018-9162.
Bibliography 123
Christoph M. Hoffmann. Robustness in geometric computations. Journal of
Computing and Information Science in Engineering, 1(2):143–155, 2001.
C.M. Hoffmann, J.E. Hopcroft, and M.S. Karasick. Towards implementing robust
geometric computations. In Proceedings of the fourth annual symposium on
Computational geometry, page 117. ACM, 1988.
C.M. Hoffmann, J.E. Hopcroft, and M.J. Karasick. Robust set operations on
polyhedral solids. Computer Graphics and Applications, IEEE, 9(6):50 –59,
nov. 1989. ISSN 0272-1716.
P. Hudak, S. L. Peyton Jones, P. L. Wadler, Arvind, B. Boutel, J. Fairbarn,
J. Fasel, M. Guzman, K. Hammond, J. Hughes, T. Johnson, R. Kieburtz, R. S.
Nikhil, W. Partain, and J. Peterson. Report on the functional programming
language Haskell, version 1.2. ACM SIGPLAN Notices, 27(5):1–164, 1992.
Annie Hui, Lucas Vaczlavik, and Leila De Floriani. A decomposition-based rep-
resentation for 3D simplicial complexes. In SGP ’06: Proceedings of the fourth
Eurographics symposium on Geometry processing, pages 101–110, Aire-la-Ville,
Switzerland, Switzerland, 2006. Eurographics Association. ISBN 30905673-36-
3.
S.P. Jones and S.L.P. Jones. Haskell 98 language and libraries: the revised report.
Cambridge Univ Pr, 2003.
Michael Karasick, Derek Lieber, and Lee R. Nackman. Efficient delaunay tri-
angulation using rational arithmetic. ACM Trans. Graph., 10(1):71–91, 1991.
ISSN 0730-0301.
Farid Karimipour, Andrew U. Frank, and Mahmoud R. Delavar. An operation-
independent approach to extend 2D spatial operations to 3D and moving ob-
jects. Proceedings of the 16th ACM SIGSPATIAL International Conference on
Advances in Geographic Information Systems (ACM GIS 2008), pages 51–56.
The Association for Computing Machinery, Inc., 2008.
R. Karinthi, K. Srinivas, and G. Almasi. A parallel algorithm for computing
polygon set operations. pages 115 –119. IEEE, apr. 1994.
J.M. Keil. Polygon decomposition. Handbook of Computational Geometry, pages
491–518, 2000.
Bibliography 124
L. Kettner, K. Mehlhorn, S. Pion, S. Schirra, and C. Yap. Classroom examples
of robustness problems in geometric computations. Computational Geometry,
40(1):61–78, 2008.
D.E. Knuth. Axioms and hulls. Springer, 1992.
H.P. Kriegel, H. Horn, and M. Schiwietz. The performance of object decomposi-
tion techniques for spatial query processing. In Advances in Spatial Databases,
pages 257–276. Springer, 1991.
Y. Kui Liu, X. Qiang Wang, S. Zhe Bao, M. Gombosi, and B. Zalik. An algorithm
for polygon clipping, and for determining polygon intersections and unions.
Computers & geosciences, 33(5):589–598, 2007.
U. Lauther. An O (n log n) algorithm for Boolean mask operations. In Proceedings
of the 18th Design Automation Conference, pages 555–562. IEEE Press, 1981.
C. Li, S. Pion, and C.K. Yap. Recent progress in exact geometric computation.
Journal of Logic and Algebraic Programming, 64(1):85 – 111, 2005. ISSN 1567-
8326.
Jyh-Ming Lien and Nancy M. Amato. Approximate convex decomposition. In
SCG ’04: Proceedings of the twentieth annual symposium on Computational
geometry, pages 457–458, New York, NY, USA, 2004. ACM. ISBN 1-58113-
885-7.
Jyh-Ming Lien and Nancy M. Amato. Approximate convex decomposition of
polygons. Computational Geometry, 35(1-2):100 – 123, 2006. ISSN 0925-7721.
Rong Liu, Hao Zhang, and James Busby. Convex hull covering of polygonal
scenes for accurate collision detection in games. In GI ’08: Proceedings of
graphics interface 2008, pages 203–210, Toronto, Ont., Canada, Canada, 2008.
Canadian Information Processing Society. ISBN 978-1-56881-423-0.
A. Margalit and G.D. Knott. An algorithm for computing the union, intersection
or difference of two polygons. Computers & Graphics, 13(2):167–183, 1989.
F. Martınez, A.J. Rueda, and F.R. Feito. A new algorithm for computing Boolean
operations on polygons. Computers & Geosciences, 35(6):1177–1185, 2009.
Bibliography 125
Kurt Mehlhorn and Stefan Naher. LEDA: A platform for combinatorial and
geometric computing. Cambridge University Press, Cambridge, 1999.
V.J. Milenkovic. Practical methods for set operations on polygons using exact
arithmetic. In Proc. 7th Canad. Conf. Comput. Geom, pages 55–60. Citeseer,
1995.
David M. Mount. Geometric intersection. In Handbook of Discrete and Compu-
tational Geometry, chapter 33, pages 615–632. CRC Press LLC, Boca, 1997.
J. Nievergelt and F. P. Preparata. Plane-sweep algorithms for intersecting geo-
metric figures. Commun. ACM, 25(10):739–747, 1982. ISSN 0001-0782.
Joseph O’Rourke. Computational geometry in C (2nd ed.). Cambridge University
Press, New York, NY, USA, 1998. ISBN 0-521-64976-5.
Joseph O’Rourke, Chi-Bin Chien, Thomas Olson, and David Naddor. A new
linear algorithm for intersecting convex polygons. Computer Graphics and
Image Processing, 19(4):384 – 391, 1982. ISSN 0146-664X.
Leonidas Palios. Decomposition Problems in Computational Geometry. PhD
thesis, Department of Computer Sceinces, Princeton University, April 1992.
A. Paoluzzi, F. Bernardini, C. Cattani, and V. Ferrucci. Dimension-independent
modeling with simplicial complexes. ACM Trans. Graph., 12(1):56–102, 1993.
ISSN 0730-0301.
Yu Peng, Jun-Hai Yong, Wei-Ming Dong, Hui Zhang, and Jia-Guang Sun. A new
algorithm for boolean operations on general polygons. Computers & Graphics,
29(1):57 – 70, 2005. ISSN 0097-8493.
Frisco Penninga. 3D Topography: Simplicial Complex-based Solution in a Spatial
DBMS. PhD thesis, The research and development center for Geo-Information
technology of Delft University of Technology, 2008.
Ken Perlin and Luiz Velho. Live paint: Painting with procedural multiscale
textures. 1995.
Eric Persoon and King-Sun Fu. Shape discrimination using fourier descriptors.
Systems, Man and Cybernetics, IEEE Transactions on, 7(3):170 –179, 1977.
ISSN 0018-9472.
Bibliography 126
Ari Rappoport. The n-dimensional extended convex differences tree (ECDT) for
representing polyhedra. In SMA ’91: Proceedings of the first ACM symposium
on Solid modeling foundations and CAD/CAM applications, pages 139–147,
New York, NY, USA, 1991. ACM. ISBN 0-89791-427-9.
P. Rigaux, M.O. Scholl, and A. Voisard. Introduction to spatial databases: with
application to GIS. Morgan Kaufmann Inc., 2002.
M. Rivero and F. R. Feito. Boolean operations on general planar polygons.
Computers & Graphics, 24(6):881 – 896, 2000. ISSN 0097-8493.
A. Sabau. Signedintersection-A new algorithm for finding the intersection of two
simple polygons. INFORMATICA, page 83, 2008.
D. Salesin, J Stolfi, and L. Guibas. Epsilon geometry: building robust algorithms
from imprecise computations. In SCG ’89: Proceedings of the fifth annual
symposium on Computational geometry, pages 208–217, New York, NY, USA,
1989. ACM. ISBN 0-89791-318-3.
H. Samet. The quadtree and related hierarchical data structures. ACM Comput-
ing Surveys, 16(2):187–260, 1984.
Hanan Samet. Foundations of Multidimensional and Metric Data Structures.
Morgan Kaufmann, 2006.
B. Schachter. Decomposition of polygons into convex sets. IEEE Transactions
on Computers, 100(27):1078–1082, 1978.
S. Schirra. Robustness and precision issues in geometric computation. Max Planck
Institut Fur InformatiK-Report-MPI I, 1998.
Markus Schneider. Spatial data types - A survey. In Markus Schneider, editor,
Spatial Data Types for Database Systems, volume 1288 of Lecture Notes in
Computer Science, pages 11–83. Springer Berlin / Heidelberg, 1997a.
Markus Schneider. Introduction. In Markus Schneider, editor, Spatial Data Types
for Database Systems, volume 1288 of Lecture Notes in Computer Science,
pages 1–9. Springer Berlin / Heidelberg, 1997b.
Markus Schneider. Spatial data types: Conceptual foundation for the design and
implementation of spatial database systems and gis. Online Tutorial, 2002.
Bibliography 127
Martijn M. Schrage, Arjan van IJzendoorn, and Linda C. van der Gaag. Haskell
ready to dazzle the real world. In Proceedings of the 2005 ACM SIGPLAN
workshop on Haskell, Haskell ’05, pages 17–26, New York, NY, USA, 2005.
ACM. ISBN 1-59593-071-X.
Rafael J. Segura and Francisco R. Feito. An algorithm for determining intersec-
tion segment-polygon in 3d. Computers & Graphics, 22(5):587 – 592, 1998.
ISSN 0097-8493.
Timos K. Sellis, Nick Roussopoulos, and Christos Faloutsos. The R+-Tree: A
dynamic index for multi-dimensional objects. In Proceedings of the 13th Inter-
national Conference on Very Large Data Bases, VLDB ’87, pages 507–518, San
Francisco, CA, USA, 1987. Morgan Kaufmann Publishers Inc. ISBN 0-934613-
46-X.
Michael Ian Shamos and Dan Hoey. Geometric intersection problems. In SFCS
’76: Proceedings of the 17th Annual Symposium on Foundations of Computer
Science, pages 208–215, Washington, DC, USA, 1976. IEEE Computer Society.
J. R. Shewchuk. Lecture notes on geometric robustness, November 2009.
J.M. Smith and N.A. Dodgson. A topologically robust algorithm for boolean
operations on polyhedral shapes using approximate arithmetic. Computer-
Aided Design, 39(2):149 – 163, 2007. ISSN 0010-4485.
Julian M. Smith. Towards robust inexact geometric computation. PhD thesis,
Computar Laboratory, University of Cambridge, St. Edmunds College, De-
cember 2009.
W.C. Thibault and B.F. Naylor. Set operations on polyhedra using binary space
partitioning trees. In Proceedings of the 14th annual conference on Computer
graphics and interactive techniques, pages 153–162. ACM, 1987.
Rodney James Thompson. Towards a Rigorous Logic for Spatial Data Represen-
tation. PhD thesis, The research and development center for Geo-Information
technology of Delft University of Technology, 2007.
Godfried Toussaint. A simple linear algorithm for intersecting convex polygons.
The Visual Computer, 1:118–123, 1985. ISSN 0178-2789.
Bibliography 128
Peter van Oosterom and Jantien Stoter. 5D data modelling: Full integration
of 2D/3D space, time and scale dimensions. In Sara Fabrikant, Tumasch Re-
ichenbacher, Marc van Kreveld, and Christoph Schlieder, editors, Geographic
Information Science, volume 6292 of Lecture Notes in Computer Science, pages
310–324. Springer Berlin / Heidelberg, 2010.
Bala R. Vatti. A generic solution to polygon clipping. Commun. ACM, 35(7):
56–63, 1992. ISSN 0001-0782.
K. Weiler and P. Atherton. Hidden surface removal using polygon area sorting. In
Proceedings of the 4th annual conference on Computer graphics and interactive
techniques, pages 214–222. ACM, 1977.
Jianning Xu. Morphological decomposition of 2-D binary shapes into modestly
overlapped octagonal and disk components. IEEE Transactions on Image Pro-
cessing, 16(2):337–348, 2007.
Chee K. Yap. Robust geometric computation. pages 653–668, 1997.
Charles T. Zahn and Ralph Z. Roskies. Fourier descriptors for plane closed curves.
Computers, IEEE Transactions on, C-21(3):269 –281, 1972. ISSN 0018-9340.
Appendix A
A Brief Introduction to Haskell
Haskell is a general purpose, purely functional programming language. A func-
tional program is a single expression, which is executed by evaluating the expres-
sion. Functional programs consist entirely of functions. Each function takes a
number of input types and returns a single output type. Constants are functions,
which always return the same value. Even the program itself is a function. Be-
cause a function call can have no other effect than producing a result, functional
programs are said to have no side effects.
Haskell provides higher-order functions, non-strict semantics, static poly-
morphic typing, user-defined algebraic datatypes, pattern-matching, list com-
prehensions, a module system, a monadic I/O system, and a rich set of primi-
tive datatypes, including lists, arrays, arbitrary and fixed precision integers, and
floating-point numbers.
Haskell has an innovative type system which supports a systematic form of
overloading and a module system. Modules provide a way to control namespaces
and to re-use software in large programs. An expression evaluates to a value and
has a static type. Values and types are not mixed in Haskell . However, the
type system allows user-defined datatypes of various sorts, and permits not only
parametric polymorphism (using a traditional Hindley-Milner type structure) but
also adhoc polymorphism, or overloading (using type classes). More information
about Haskell can be found in [Jones and Jones, 2003].
129
Appendix A - A Brief Introduction to Haskell 130
A.1 High Order Functions-Functors
A functional or high order function is a function whose arguments are functions
or whose result is a function. The map function illustrates how a high order
function works. For example, map which is a higher order function (Functor)
defined as [Jones and Jones, 2003];
class Functor f where
fmap :: (a -> b) -> f a -> f b
instance Functor [] where
fmap = map
map is an instance of Functor typeclass which represents types which can be
mapped over. It defines only one function fmap whose arguments are a function
from one type to another type and a functor applied to one type and returns a
functor applied with another type. This becomes clear if we look at the type
signature of map
map : : ( a −> b) −> [ a ] −> [ b ]
which is an instance of Functor class with f being [] (a list type constructor), so
map is a Functor over lists. It takes a function from one type to another type and
a list of one type and returns a list of another type. The example below shows
the working of map;
Prelude> map (+2) [ 1 , 7 , 18 , −29]
[3 ,9 ,20 , −27 ]
High order functions distinguish functional programming from other program-
ming languages. The partial application of functions and the possibility to define
functions of functions permit abstract specifications. As a result generic code can
be written that is modular and highly reusable.
A.2 Haskell Types
Haskell’s static type system defines the formal relationship between types and
values. The static type system ensures that Haskell programs are type safe, All
type errors are detected at compile-time. The type system also ensures that
user-supplied type signatures are correct.
Appendix A - A Brief Introduction to Haskell 131
A.2.1 Generic data types
Haskell provides a way to define type synonyms; i.e. names for commonly used
types. Type synonyms are created using a type declaration. Some examples are
shown below.
type Point = (Int , Int)
type String = [Char]
type Person = (Name ,Address)
In Haskell, the newtype declaration creates a new type from an existing one.
For example, natural numbers can be represented by the type Integer using the
following declaration (refer gentle introduction). For example, AHD level is rep-
resented by the type Int as shown below,
newtype Level = Level Int
A.2.2 Algebraic data types
Algebraic data type is a datatype each of whose values is data from other datatypes
wrapped in one of the constructors of the datatype. Any wrapped datum is an ar-
gument to the constructor. In contrast to other datatypes, the constructor is not
executed and the only way to operate on the data is to unwrap the constructor
using pattern matching. For example, the tree data structure of AHD is defined
as an algebraic data type;
data Tree p = Node {_this :: p, _childtrees :: [Tree p
], _level :: Level}
| LeafNode {_this :: p, _level :: Level}
| UniNode { _this :: p, _childtrees :: [Tree p
], _level :: Level}
| NullTree {}
Here, the Tree is a type type constructor, p is a type variable, Node,
LeafNode, UniNode, and NullTree are data constructors.
A.3 Classes in Haskell
Classes are a specific feature of the Haskell programming language. A class in
Haskell is a collection of types that support certain overloaded operations called
Appendix A - A Brief Introduction to Haskell 132
methods. Haskell provides a number classes that are built-in to the language
shown in Figure A.1.
Eq
/ \
/ \
Ord \ Text
/ \ \ /
/ \ \ /
/ Enum Num
/ \ / \
Ix \ / \
\ Real Fractional
\ / \ / \
\ / \ / \
Integral RealFrac Floating
\ /
\ /
RealFloat
Figure A.1: Hierarchy of Haskell classes
On the top of the hierarchy there is the class Eq. It defines the collection of
types with which the equality of two elements can be tested. The class declaration
of Eq is given below. It consists of a class Name followed by a signature. The
signature is a list of names and their types. The equality operation (==) takes
two types and returns a Boolean value.
class Eq a where
(==) :: a -> a -> Bool
In order to make an arbitrary type a member of the class Eq an implementation
has to be provided. Members of a type class are called instances. Haskell
defines instances for Int, Float, Bool and Char for Eq. Furthermore a default
implementation for the class Eq is given that can be overwritten with a new
implementation.
The equality between two different types may differ. The equality of two
natural numbers may be implemented by simply comparing their values, while
the equality of two strings may be based on comparing the length of the two
strings. Different instances may be defined. The appropriate implementation
will be overloaded for the corresponding type. This way type classes implement
Appendix A - A Brief Introduction to Haskell 133
the ad-hoc polymorphism mentioned earlier.
instance Eq Int where
(==) a b = a == b
instance Eq String where
(==) a b = length a == length b
For example here is the declaration for the type class Points;
data Point = Point Float Float
class Points p where
getX : : p −> Float
getY : : p −> Float
instance Points Point where
getX ( Point x y ) = x
getY ( Point x y ) = y
The first line declares the name of the class, and constraints on what types
can be instances of this class. Then is the list of type declarations of the class
functions that must be implemented in order to be an instance of the class.
Haskell also offers a mechanism of inheritance, similar to object oriented pro-
gramming languages. The class Ord can be derived from the class Eq. Ord defines
the class of ordered types. The class Ord defines the operations to compare types
like <, <=, >, >= . The definition of equality, defined by the == operator is
inherited of the class Eq.
class Eq a => Ord a where
(<),(<=),(>),(>=) :: a -> a -> Bool
...
In the sample code above the => operator indicates that Ord is derived from Eq.
The operator => refers to the context of a class. In the example that means for
any type a that is declared and implemented belonging to Ord there has to be
also a declaration and implementation belonging to Eq.
Appendix B
Prototype Application
The prototype application for AHD is developed in Haskell using leksah1, an open
source Haskell IDE, that uses Glasgow Haskell Compiler (ghc) and Gtk . The
GUI is developed in Glade2 (RAD tool for quick and easy development of user
interfaces) using the gtk2hs3 binding for Haskell.
A screen-shot of the AHD prototype GUI is shown in Figure B.1 and Table
B.1 describes different parts of the application GUI.
Table B.1: Description of AHD prototype GUI
1 Select the region to draw
2 Drawing area where the region polygons are drawn
3 Here the text view of the application status and region statistics are provided
4 Set width of the region edges
5 Set color for the region
6 Checkbox for directed or undirected edges
7 File selection or browser button to load data from files
8 Saving regions for further processing
9 Button provides the AHD build function
10 Specify level of detail for build
11 Button provides the AHD eval function
12 Specify level of detail for eval
13 Buttons to perform operations on the two regions represented by AHD
14 Buttons for zooming operations
1http://leksah.org/2http://glade.gnome.org/3http://www.haskell.org/haskellwiki/Gtk2Hs
135
Figure B.1: Screen-shot of AHD prototype
Appendix C
Haskell Code
module AHD. Algebra where
import AHD.Geom
import AHD. Leve l
class AHD inp inpRep prim where
bu i ld : : Leve l −> Level−> inp −> [ inpRep prim ]
eva l : : Leve l −> [ inpRep prim ] −> inp
class Input inp where
i s S p l i t t a b l e : : inp −> Bool
merge : : inp −> inp −> inp
sp l i t : : inp −> [ inp ]
class InputPrim inp prim where
makePrim : : inp −> prim
mergePrim : : prim −> prim −> [ inp ]
class Prim prim where
c c i n t : : prim −> prim −> prim
mergeComp : : prim −> [ prim ] −> [ prim ]
class AHDOperations rep prim where
cn int : : prim −> rep prim−> [ rep prim ]
137
Appendix C - Haskell Code 138
g in t : : rep prim−> rep prim −> [ rep prim ]
intReg : : [ rep prim ] −> [ rep prim ] −> [ rep prim ]
compReg : : [ rep prim]−> [ rep prim ]
uniReg : : [ rep prim ] −> [ rep prim ] −> [ rep prim ]
di fReg : : [ rep prim ] −> [ rep prim ] −> [ rep prim ]
symDReg : : [ rep prim ] −> [ rep prim ] −> [ rep prim ]
pip : : BPoint −> rep prim −> Bool
p i r : : BPoint −> [ rep prim ] −> Bool
Appendix C - Haskell Code 139
module AHD.Geom where
import Data . List
import Language . ConvsAlg
import AHD. ListAux
−−D e f i n i t i o n o f geometr ic data t y p e s
data BPoint = Bpt { w , x , y : : Integer}deriving (Show, Eq, Ord)
unPoint ( Bpt w x y ) = (w, x , y )
data Ring = Ring { l p o i n t s : : [ BPoint ]}deriving (Show, Eq, Ord)
unRing ( Ring r ing ) = r ing
data Poly = Poly { l r i n g s : : [ Ring ]}deriving (Show, Eq, Ord)
unPoly ( Poly poly ) = poly
data CHull = CHull { c h u l l : : [ BPoint ]}deriving (Show, Eq, Ord)
unChull ( CHull ch ) = ch
data Reg = Reg { r eg : : [ Poly ]}deriving (Show, Eq, Ord)
unReg (Reg reg ) = reg
data LinSeg = LinSeg { sp , ep : : BPoint}deriving (Show, Eq, Ord)
s2e ( LinSeg p1 p2 ) = ( p1 , p2 )
e2s ( p1 , p2 ) = ( LinSeg p1 p2 )
data RegStats = RegStats { npo lys : : Int , n r i n g s : : Int ,
npts : : Int}deriving (Show)
Appendix C - Haskell Code 140
−− For s t r i p t r e e
type Curve = [ BPoint ]
type Triang le = ( BPoint , BPoint , BPoint )
data St r i p = St r i p { spx : : Double , spy : : Double ,
epx : : Double , epy : : Double , wl : : Double , wr : :
Double}deriving (Show, Read)
−− n dimensiona l AHD
−− An n dimnesiona l p o i n t d e f i n e d as a l i s t o f p o i n t s
type ND BPoint = [ BPoint ]
−− Define a v e r t e x as a p o i n t
type Vertex = ND BPoint
−− An n− s imp lex i s d e f i n e d as a l i s t o f n dimensiona l
v e r t e x e s
type Simplex = [ ND BPoint ]
−− Canonical r e p r e s e n t a t i o n o f a n−s imp lex
type CnSimplex = ( Simplex , Bool )
−−Auxi lary geometr ic f u n c t i o n s
−− i s the edge on the g iven reg ion edges
isEdgeOnReg : : ( BPoint , BPoint ) −> [ ( BPoint , BPoint ) ] −>(Bool , ( BPoint , BPoint ) )
isEdgeOnReg e lhe = i f ( f e == [ ] ) then ( False , nu l l edge )
else (True , head f e )
where
f e = [ x | x <− lhe , ( isEdgeOn e x ) && (
isInRange e x ) ]
nu l l edge = ( ( Bpt 0 0 0) , ( Bpt 0 0 0) )
isInRange : : ( BPoint , BPoint ) −> ( BPoint , BPoint ) −> Bool
Appendix C - Haskell Code 141
isInRange ( p3 , p4 ) ( p1 , p2 ) = ( ( ( isbw p3x ( p1x , p2x ) ) && (
isbw p4x ( p1x , p2x ) ) ) && ( ( isbw p3y ( p1y , p2y ) ) && ( isbw
p4y ( p1y , p2y ) ) ) )
where
(p1w , p1x , p1y ) = h2c . unPoint $ p1
(p2w , p2x , p2y ) = h2c . unPoint $ p2
(p3w , p3x , p3y ) = h2c . unPoint $ p3
(p4w , p4x , p4y ) = h2c . unPoint $ p1
−− i s v a l u e in between the g iven range
isbw : : Double −> (Double , Double) −> Bool
isbw v ( v1 , v2 ) = ( ( v >= l r ) && ( v <= ur ) )
where
l r = min v1 v2
ur = max v1 v2
−− Given e1 and e2 , f u n c t i o n t e l l s whether e1 i s on e2
isEdgeOn : : ( BPoint , BPoint ) −> ( BPoint , BPoint ) −> Bool
isEdgeOn ( a , b) e = ( ( isPointOn a e ) && ( isPointOn b e ) )
−− Given an Edge and reg ion edges , f u n c t i o n t e l l s whether
the edge i n s i d e reg ion or not
i sEdgeIn : : [ ( BPoint , BPoint )]−> ( BPoint , BPoint ) −> Bool
i sEdgeIn re ( a , b ) = ( ( isPointInOn a re ) && ( isPointInOn
b re ) )
isPointInOn : : BPoint −> [ ( BPoint , BPoint ) ] −> Bool
i sPointInOn p l e = ( a l l(==True) i n t e s t )
where
i n t e s t = map (\y −> i sPointLeftOn y p) l e
i sPointLeftOn : : ( BPoint , BPoint )−> BPoint −> Bool
i sPointLeftOn (a , b) p = ( ( twiceOfArea p a b) >= 0)
−− Given a p o i n t and an edge , f u n c t i o n t e l l s whether p o i n t
i s ”On” the edge
isPointOn : : BPoint −> ( BPoint , BPoint )−> Bool
isPointOn p ( a , b) = ( ( twiceOfArea p a b) == 0)
−− i s h u l l 1 i n t e r s e c t s wi th h u l l 2
i s H u l l I n t : : CHull −> CHull −> Bool
i s H u l l I n t ( CHull h1 ) ( CHull h2 ) = (any (== True) onTest )
where
Appendix C - Haskell Code 142
onTest = map fst (map (\y −> isEdgeOnReg y
h2e ) h1e )
h1e = mkEdges h1
h2e = mkEdges h2
−−−− Given a l i s t o f d e l t a reg ion edges and a convex h u l l ,
the f u n c t i o n
−−−− s e p a r a t e s the c l o s e d r e g i o n s
sepReg : : [ ( BPoint , BPoint ) ] −> CHull −> [ [ ( BPoint , BPoint
) ] ]
sepReg [ ] ch = [ ]
sepReg dr ( CHull [ ] ) = sepHoles dr −− means s e p a r a t e h o l e s
sepReg dr ch = i f ( s e s == [ ] ) then sepReg dr ( CHull [ ] )
else sep s e s dre ch
where
s e s = f i l t e r ( i sS ta r tE ch ) dr −− l i s t o f
s t a r t i n g edges
dre = dr \\ s e s −− remove s t a r t i n g edges from
l i s t o f d e l t a reg ion edges
−− Given a l i s t o f s t a r t i n g edges , d e l t a reg ion edges and
convex h u l l , f u n c t i o n
−− s e g r e g a t e s the c l o s e d r e g i o n s
sep : : [ ( BPoint , BPoint )]−> [ ( BPoint , BPoint ) ] −> CHull
−> [ [ ( BPoint , BPoint ) ] ]
sep [ ] [ ] = [ ]
sep [ ] dr ch = sepHoles dr
sep ( x : xs ) dr ch = ( x : sex ) : sep xs udr ch where
udr = dr \\ sex −− update dr edges by removing
succeed ing edges
sex = succEdges x dr ch −− succeed ing edges o f x
−− Given an edge (a , b ) and a convex h u l l , the f u n c t i o n
t e l l s whether the g iven edge
−− i s a s t a r t i n g edge f o r a c a n v i t y or not
i s S ta r tE : : CHull −> ( BPoint , BPoint ) −> Bool
i s S ta r tE ch (a , b) = any (== a ) ( unChull ch )
−− Given a s t a r t i n g edge , l i s t o f reg ion edges and convex
h u l l , the f u n c t i o n
Appendix C - Haskell Code 143
−− g i v e s the succd ing edges o f the g iven edge
succEdges : : ( BPoint , BPoint ) −> [ ( BPoint , BPoint ) ] −>CHull −> [ ( BPoint , BPoint ) ]
succEdges e [ ] ch = [ ]
succEdges e l e ch = i f ( ne == [ ] ) then [ ]
else (head ne ) : succEdges (head ne ) n l ch
where
ne = ( nxtEdge e l e ch )
n l = l e \\ne
−− Given an edge , reg ion edges and convex h u l l , f u n c t i o n
r e t u r n s the next edge
nxtEdge : : ( BPoint , BPoint ) −> [ ( BPoint , BPoint )]−> CHull
−> [ ( BPoint , BPoint ) ]
nxtEdge ( a , b) re ch = [ e | e <− re , ( f s t e == b) && (
i sS ta r tE ch e == False ) ]
−− Given l i s t o f h o l e edges , f u n c t i o n s e p a r a t e s the h o l e
edge r e g i o n s
sepHoles : : [ ( BPoint , BPoint )]−> [ [ ( BPoint , BPoint ) ] ]
sepHoles [ ] = [ ]
sepHoles ( ( a , b ) : xs ) = i f ( eee == [ ] ) −− means unwanted
edge
then [ ( a , b ) ] : sepHoles xs
else hr : sepHoles nre
where
se = (a , b)
eee = f i l t e r (\y −> snd y == a ) xs
ee = head eee
r r e = xs \\ [ ee ]
hre = ( ( se : ( getIntEdges se ee r r e ) ) ++ [ ee ] )
i r e = r r e \\ hre
e inhr = f i l t e r (\ e −> ( i sEdgeIn hre e ) == False )
i r e
hr = hre ++ e inhr
nre = i r e \\ e inhr
−− Given a s t a r t i n g edge s , ending edge e , and a l i s t o f
reg ion edges re , f u n c t i o n computes the
Appendix C - Haskell Code 144
−− l i s t o f i n t e r m e d i a t e edges connect ing s and e .
getIntEdges : : ( BPoint , BPoint ) −> ( BPoint , BPoint ) −> [ (
BPoint , BPoint ) ] −> [ ( BPoint , BPoint ) ]
getIntEdges s e [ ] = [ ]
getIntEdges s e re = i f ( ( ne == [ ] ) | | ( nx == e ) ) then [ ]
else nx : getIntEdges nx e n l where
ne = ( f i l t e r (\y −> f s t y == snd s ) re )
nx = head ne
n l = re \\ne
−− Given a reg ion and a h u l l edges the f u n t i o n t e l l
whether the g iven reg ion i s a h o l e or not
i sRegHole : : [ ( BPoint , BPoint ) ] −> [ ( BPoint , BPoint )]−>Bool
i sRegHole he re = (any (==True) onTest )
where
onTest = map fst (map (\y −> isEdgeOnReg y he ) re )
−− Function t h a t computes t w i c e the area o f t r i a n g l e by
g iven t h r e e p o i n t s abc
twiceOfArea : : BPoint −> BPoint −> BPoint −> Integer
twiceOfArea ( Bpt b0 b1 b2 ) ( Bpt c0 c1 c2 ) ( Bpt a0 a1 a2 )
= a0 * ( b1* c2 − b2* c1 ) −b0 *( a1* c2 − a2* c1 ) +c0 * ( a1*
b2 − a2*b1 )
−−Function t h a t computes t w i c e the area o f t r i a n g l e by
g iven t h r e e p o i n t s abc
twiceOfArea2 : : (Num a ) => ( a , a , a ) −> ( a , a , a ) −> ( a ,
a , a ) −> a
twiceOfArea2 ( b0 , b1 , b2 ) ( c0 , c1 , c2 ) ( a0 , a1 , a2 ) = a0 * ( b1
* c2 − b2* c1 ) − b0 *( a1* c2 − a2* c1 ) + c0 * ( a1*b2 − a2*
b1 )
h2c : : (NumConvs Integer ) => ( Integer , Integer , Integer ) −>(Double , Double , Double)
h2c (w, x , y ) = ( toDouble w/ toDouble w, toDouble x/ toDouble
w, toDouble y/ toDouble w)
l i s t e d g e s 2 R i n g : : [ ( BPoint , BPoint ) ] −> Ring
l i s t e d g e s 2 R i n g lp = Ring (nub . f l a t t e n $ lp )
−− Helping f u n c t i o n t h a t r e v e r s e s the o r i e n t a t i o n o f a
Appendix C - Haskell Code 145
g iven edge
revOr ient : : ( BPoint , BPoint )−> ( BPoint , BPoint )
revOr ient ( a , b ) = (b , a )
−−Edge to p o i n t n o t a t i o n
edge2po ints : : ( BPoint , BPoint ) −> [ BPoint ]
edge2po ints ( p1 , p2 ) = [ p1 , p2 ]
−− Given a l i s t o f edges , t h i s f u n c t i o n s p l i t s c l o s e d
r e g i o n s
sepR : : [ ( BPoint , BPoint ) ] −> [ [ ( BPoint , BPoint ) ] ]
sepR [ ] = [ ]
sepR ( ( a , b) : xs ) = i f ( e e l == [ ] ) then [ ( a , b ) ] : sepR xs
else c re : sepR nre where
se = (a , b) −− s t a r t i n g edge
e e l = f i l t e r (\y −> snd y == a ) xs −−p r o b a b l e ending edge / s l i s t
ee = head e e l −− ending edge
r r e = xs \\ [ ee ] −− remaing reg ion edges
when ee i s removed
c re = ( ( se : ( getIntEdges se ee r r e ) ) ++ [
ee ] ) −− c l o s e d reg ion edges
nre = r r e \\ c re −− new reg ion edges f o r
next i t e r a t i o n
−−edge to v e r t e x n o t a t i o n
edge2vertex : : [ ( BPoint , BPoint )]−> [ BPoint ]
edge2vertex [ ] = [ ]
edge2vertex ( ( p1 , p2 ) : [ ] ) = [ p1 ]
edge2vertex ( ( p1 , p2 ) : e s ) = i f se == [ ] then [ ]
else p1 : edge2vertex ( se++ne l )
where
se = sucedge ( p1 , p2 ) es
ne l = es \\ ( se ++[(p1 , p2 ) ] )
−− Given an edge and l i s t o f edges r e t u r n s the next edge
sucedge ( p1 , p2 ) l e = [ e | e <− l e , ( f s t e == p2 ) ]
revRing : : Ring −> Ring
revRing ( Ring r ing ) = Ring ( reverse r i ng )
xOrd : : BPoint −> BPoint −> Ordering
Appendix C - Haskell Code 146
xOrd ( Bpt z1 x1 y1 ) ( Bpt z2 x2 y2 ) = i f ( x1 == x2 ) then EQ
else i f ( x1 > x2 ) then GT
else LT
yOrd : : BPoint −> BPoint −> Ordering
yOrd ( Bpt z1 x1 y1 ) ( Bpt z2 x2 y2 ) = i f ( y1 == y2 ) then EQ
else i f ( y1 > y2 ) then GT
else LT
Appendix C - Haskell Code 147
module AHD. Leve l where
newtype Level = Level Int
deriving (Eq, Ord , Show, Read , Num)
z e roLeve l = Level 0
−− Functions on Leve l
unLevel ( Leve l i ) = i
isEven ( Leve l l ) = even l
i n c l = l+1
dec l = l−1
Appendix C - Haskell Code 148
module AHD. AHDTree where
import AHD. Algebra
import AHD.Geom
import AHD. ListAux
import Data . List
import ConvexHull
import AHD. Leve l
import AHDOps
−− Tree data s t rucure , every node o f which co nt a in s a
p r i m i t i v e o f type p
data Tree p = Node { t h i s : : ! p , c h i l d t r e e s : : ! [ Tree p ] ,
l e v e l : : ! Leve l }| LeafNode { t h i s : : ! p , l e v e l : : ! Leve l }| UniNode { t h i s : : ! p , c h i l d t r e e s : : ! [ Tree p ] ,
l e v e l : : ! Leve l }| NullTree {}
deriving (Eq, Show, Read)
−− data s t r u c t u r e t h a t h o l d s t r e e s t a t i s t i c s
data TreeStats = TreeStats { nNodes : : ! Int , nLeve l s : : !
Int}deriving (Eq, Show, Read)
unTreeStats ( TreeStats nn nl ) = (nn , n l )
class AHDTree t p where
c h i l d t r e e s : : t p −> [ t p ]
l e v e l : : t p −> Level
t h i s : : t p −> p
instance AHDTree Tree CHull where
c h i l d t r e e s = c h i l d t r e e s
l e v e l = l e v e l
t h i s = t h i s
Appendix C - Haskell Code 149
instance AHD Ring Tree CHull where
bu i ld mlev c l e v r ing = i f ( d r ing s == [ ] ) | | ( mlev ==
c l e v ) then [ ( LeafNode p c l e v ) ]
else [ Node p ( concat (map ( bu i ld mlev n l ) d r ing s )
) c l e v ]
where
p = makePrim r ing
dr ing s = sp l i t r i ng
n l = inc c l e v
eva l mlev [ ( LeafNode ( CHull p ) l ) ] = Ring p
eva l mlev [ ( Node ( CHull p ) l t s l ) ] = i f ( mlev == l )
then Ring p
else ( fo ld l merge ( Ring p) l l r )
where
l l r = (map ( eva l mlev ) ( toLL l t s ) )
instance AHD Poly Tree CHull where
bu i ld mlev c l e v ( Poly [ ] ) = [ ]
bu i ld mlev c l e v ( Poly poly ) = [ fo ld l mergTrees (head
l t r e e s ) ( t a i l l t r e e s ) ]
where
l t r e e s = concat (map ( bu i ld mlev c l e v ) poly )
eva l mlev [ ( LeafNode ( CHull p ) l ) ] = Poly [ Ring p ]
eva l mlev [ ( Node ( CHull p ) l t s l ) ] = i f ( mlev == l )
then Poly [ Ring p ]
else ( fo ld l merge ( Poly [ Ring p ] ) lp )
where
lp = (map ( eva l mlev ) ( toLL l t s ) )
instance AHD Reg Tree CHull where
bu i ld mlev c l e v (Reg [ ] ) = [ ]
bu i ld mlev c l e v (Reg reg ) = concat (map ( bu i ld mlev
c l e v ) reg )
eva l mlev l t s = Reg (map ( eva l mlev ) ( toLL l t s ) )
instance Input Ring where
Appendix C - Haskell Code 150
merge ( Ring r0 ) ( Ring r1 ) = Ring ( edge2vertex nmre
)
where
er0 = mkEdges r0
er1 = mkEdges r1
oes = [ ( ce0 , ce1 ) | ce0 <− er0 , ce1<− er1 ,
isEdgeOn ce1 ce0 ] −− edges t h a t o v e r l a p
re0 = er0 \\ (map fst oes ) −− remaining edges
o f 0 , s u b t r a c t o v e r l a p p i n g edge
de1 = er1 \\ (map snd oes ) −− d e l t a edges o f
h u l l 1 , s u b t r a c t o v e r l a p p i n g edge
re1 = map revOr ient de1 −− r e v e r s e o r i e n t a t i o n
o f deta edges o f h u l l 1
nmes = concat (map mergEdges oes )−−merge
o v e r l a p p i n g edges
nmre = re1 ++ re0 ++ nmes −− new merged reg ion
edges
sp l i t r i ng = map l i s t e d g e s 2 R i n g dre −−d e l t a reg ion
edges
where
dre = sepReg fde ( CHull (nub( ( unChull ch ) ++
dph ) ) ) −− add dph to h u l l po in ts , nub and
then s e p e r a t e c l o s e d r e g i o n s
dph = nub . f l a t t e n $ ( de \\ fde ) −− p o i n t s not
in d e l t a reg ion but on h u l l
fde = f i l t e r (\y −> f s t ( isEdgeOnReg y he )==
False ) de −− d e l t a edges not on convex h u l l
de = ( mkEdges . unRing $ r ing ) \\ he −− d e l t a
edges
he = mkEdges . unChull $ ch −−h u l l edges
ch = quickHul l r i ng
instance InputPrim Ring CHull where
makePrim r ing = quickHul l r i ng
mergePrim ( CHull r0 ) ( CHull r1 ) = [ Ring (
edge2vertex nmre ) ]
Appendix C - Haskell Code 151
where
er0 = mkEdges r0
er1 = mkEdges r1
oes = [ ( ce0 , ce1 ) | ce0 <− er0 , ce1<− er1 ,
isEdgeOn ce1 ce0 ] −− edges t h a t o v e r l a p
re0 = er0 \\ (map fst oes ) −− remaining edges
o f 0 , s u b t r a c t o v e r l a p p i n g edge
de1 = er1 \\ (map snd oes ) −− d e l t a edges o f
h u l l 1 , s u b t r a c t o v e r l a p p i n g edge
re1 = map revOr ient de1 −− r e v e r s e o r i e n t a t i o n
o f deta edges o f h u l l 1
nmes = concat (map mergEdges oes )−−merge
o v e r l a p p i n g edges
nmre = re1 ++ re0 ++ nmes −− new merged reg ion
edges
instance Input Poly where
merge ( Poly p0 ) ( Poly p1 ) = Poly ( mergeRingLists
p0 p1 )
instance InputPrim Poly CHull where
makePrim poly = quickHul l (head . unPoly $ poly )
−−Some a u x i l l a r y f u n c t i o n s
mergeRingLists : : [ Ring ] −> [ Ring ] −> [ Ring ]
mergeRingLists [ ] r l = r l
mergeRingLists r l [ ] = r l
mergeRingLists ( r0 : [ ] ) ( r1 : [ ] ) = ( mergeRings r0 r1 )
mergeRingLists ( r0 : [ ] ) ( r1 : l r ) = ( mergeRings r0 r1 ) ++ (
map revRing l r )
mergeRingLists ( r0 : r0 s ) ( r1 : [ ] ) = ( mergeRings r0 r1 ) ++
r0s
mergeRingLists ( r0 : r0 s ) ( r1 : r1 s ) = ( mergeRings r0 r1 )++r0s
++ (map revRing r1s )
−− Merge two r i n g s
mergeRings : : Ring −> Ring −> [ Ring ]
Appendix C - Haskell Code 152
mergeRings ( Ring r0 ) ( Ring r1 ) = i f ( oes == [ ] ) then ( Ring
r0 ) : [ Ring ( reverse r1 ) ]−−means a r1 i s a h o l e
else map Ring (map edge2vertex ( sepR nmre ) )−−[ Ring (
edge2r ing nmre) ]
where
er0 = mkEdges r0
er1 = mkEdges r1
oes = [ ( ce0 , ce1 ) | ce0 <− er0 , ce1<− er1 ,
isEdgeOn ce1 ce0 ] −− edges t h a t o v e r l a p
re0 = er0 \\ (map fst oes ) −− remaining edges
o f 0 , s u b t r a c t o v e r l a p p i n g edge
de1 = er1 \\ (map snd oes ) −− d e l t a edges o f
h u l l 1 , s u b t r a c t o v e r l a p p i n g edge
re1 = map revOr ient de1 −− r e v e r s e o r i e n t a t i o n
o f deta edges o f h u l l 1
nmes = concat (map mergEdges oes )−−mergEdges (
f s t . head $ oes ) ( snd . head $ oes )−−merge
o v e r l a p p i n g edges
nmre = re1 ++ re0 ++ nmes −− new merged reg ion
edges
−− Given two h u l l edges a and b , a i s on b , i n p t i o n merges
a wi th b
mergEdges : : ( ( BPoint , BPoint ) , ( BPoint , BPoint ) )−> [ (
BPoint , BPoint ) ]
mergEdges ( ( a1 , a2 ) , ( b1 , b2 ) ) = i f ( ( a1 , a2 ) == ( b1 , b2 ) )
−− l i n e s o v e r l a p
then [ ] −− re turn empty so t h a t the edge i s removed
only
else i f ( a1 == b1 && a2 /= b2 ) then [ ( b2 , a2 ) ]
else i f ( a2== b2 && a1 /= b1 ) then [ ( a1 , b1 ) ]
else ( a1 , b1 ) : [ ( b2 , a2 ) ]
−−Graft second t r e e on f i r s t t r e e
gra f tTree : : Tree CHull−> Tree CHull−> Tree CHull
g ra f tTree ( LeafNode y ly ) ( LeafNode x lx ) = (Node x [ (
incTLevel ( inc lx ) ( LeafNode y ly ) ) ] l x )
g ra f tTree ( LeafNode y ly ) (Node x xts lx ) = (Node x ( xts
Appendix C - Haskell Code 153
++ [ ( incTLevel ( inc lx ) ( LeafNode y ly ) ) ] ) l x )
g ra f tTree (Node y yts ly ) ( LeafNode x lx ) = (Node x [ (
incTLevel ( inc lx ) (Node y yts ly ) ) ] l x )
g ra f tTree (Node y yts ly ) (Node x xts lx ) = (Node x ( xts
++ [ ( incTLevel ( inc lx ) (Node y yts ly ) ) ] ) l x )
−−−− Given two t r e e s , one paraent and one c h i l d ring ,
i n p t i o n merges the two
mergTrees : : Tree CHull −> Tree CHull −> Tree CHull
mergTrees ( LeafNode p l ) t1 = gra f tTree t1 ( LeafNode p l )
mergTrees (Node p ct l ) t1 = gra f tTree t1 (Node p ct l )
−−Given a Tree , i n p t i o n r e s e t s the l e v e l o f the t r e e
incTLevel : : Level−> Tree CHull−> Tree CHull
incTLevel l ( LeafNode p ) = ( LeafNode p ( inc l ) )
incTLevel l (Node p l t s ) = (Node p (map ( incTLevel ( inc
l ) ) l t s ) ( inc l ) )
−−Given a Tree , i n p t i o n r e s e t s the l e v e l o f the t r e e
decTLevel : : Level−> Tree CHull−> Tree CHull
decTLevel l ( LeafNode p ) = ( LeafNode p l )
decTLevel l (Node p l t s ) = (Node p (map ( decTLevel ( dec
l ) ) l t s ) l )
process2Edges : : Leve l −> [ Tree CHull ] −> [ [ ( BPoint ,
BPoint ) ] ]
process2Edges l l t s = map ( proces sTree l ) l t s
proces sTree : : Leve l −> Tree CHull −> [ ( BPoint , BPoint ) ]
proces sTree l ev ( LeafNode ( CHull ch ) l ) = mkEdges ch
proces sTree l ev t = i f ( l e v == l e v e l t ) then mkEdges ch
else merg che l l e
where
( CHull ch ) = t h i s t
l c t = c h i l d t r e e s t
che = mkEdges ch
l l e = map ( processTree l ev ) l c t
merg : : [ ( BPoint , BPoint ) ] −> [ [ ( BPoint , BPoint ) ] ] −> [ (
Appendix C - Haskell Code 154
BPoint , BPoint ) ]
merg re [ ] = re
merg re ( x : xs ) = merg im xs
where
im = mergReg re x
mergReg : : [ ( BPoint , BPoint ) ] −> [ ( BPoint , BPoint ) ] −>[ ( BPoint , BPoint ) ]
mergReg re [ ] = re
mergReg re ( ( x , y ) : xs ) = i f b then mergReg ure xs
else (y , x ) : mergReg re xs
where
(b , e ) = isEdgeOnReg (x , y ) re
ues = ( updateEdge e (x , y ) ) −− update e
wi th ( x , y )
ure = ( re \\ [ e ] ) ++ ues
−− Given an edge a and another edge b , f u c n t i o n updates a
wi th b
updateEdge : : ( BPoint , BPoint ) −> ( BPoint , BPoint ) −> [ (
BPoint , BPoint ) ]
updateEdge ( a , b) ( x1 , x2 ) = i f ( ( a , b) == ( x1 , x2 ) ) −− l i n e s
o v e r l a p
then [ ] −− re turn empty so t h a t the edge i s removed
only
else i f ( a == x1 && b /= x2 ) then [ ( x2 , b ) ]
else i f (b== x2 && a /= x1 ) then [ ( a , x1 ) ]
else ( a , x1 ) : [ ( x2 , b) ]
Appendix C - Haskell Code 155
module AHDOps where
import AHD. Algebra
import AHD.Geom
import Data . List
import ConvexHull
import AHD. ListAux
import AHD. AHDTree
import AHD. Leve l
instance Prim CHull where
c c i n t ( CHull p1 ) ( CHull p2 ) = i f ( p1\\p2 == [ ] )
then CHull p1
else i f ( p1\\pinp2 == [ ] ) then CHull p1 −−p1
i s ns ide p2
else i f ( p2\\pinp1 == [ ] ) then CHull p2 −−p2 i s
i n s i d e p1
else i f ( ( length . unChull $ ip ) <=2)
then CHull [ ] −− s p e c i a l case i f i n t e r s e c t i o n i s
not a c l o s e d reg ion
else ip
where
( CHull chp ) = quickHul l ( Ring (union p1 p2 ) )
−− convex h u l l o f both po lygons
che = mkEdges chp −− convex h u l l edges
ine1 = (mkEdges p1 ) \\ che −− p o s s i b l e
i n t e r s e c t i n g edges o f po ly 1
ine2 = (mkEdges p2 ) \\ che −− p o s s i b l e
i n t e r s e c t i n g edges o f po ly 1
vip = f i p $ p i e ine1 ine2 −−p a i r a c t u a l l y
i n t e r s e c t i n g egdes and f i n d i n t p o i n t s
pinp2 = pointsInOrOn ( CHull p1 ) ( CHull p2 ) −−p o i n t s o f p1 in p2
pinp1 = pointsInOrOn ( CHull p2 ) ( CHull p1 ) −−p o i n t s o f p1 in p2
ip = quickHul l ( Ring ( pinp1 ++ vip ++ pinp2 ) )
Appendix C - Haskell Code 156
−− t ake q u i c k h u l l f o r v a l i d o r i e t a t i o n
instance AHDOperations Tree CHull where
cn int ch ( LeafNode x l ) = [ LeafNode ( c c i n t ch x ) l ]
cn in t ch (Node x xts l ) = i f ( ohi == CHull [ ] ) then
[ ] −− means no i n t e r s e c t i o n
else [ ( Node ohi ( concat (map ( cn in t ohi ) xts ) ) l ) ]
where
ohi = ( c c i n t ch x ) −− outer
h u l l i n t e r s e c t i o n
g in t t ( LeafNode x lx ) = cn int x t
g in t ( LeafNode x lx ) t = cn int x t
g in t t1 t2 = i f ( o i == [ ] ) then [ ]
else [ fo ld l mergTrees (head o i ) l t ]
where
x = t h i s t1
l c t = c h i l d t r e e s t1
o i = cn int x t2
l t = concat (map ( g in t (head
o i ) ) l c t )
intReg r1 r2 = concat [ ( g in t p x ) | p<− r1 , x<− r2 ]
compReg l p t = i f ( length l p t ==1) && ( i s Un i . head $
l p t ) then dts
else [ UniNode ( uniHul l ) i t s 0 ]
where
i t s = map ( incTLevel 0) l p t
dts = map ( decTLevel 0) ( t a i l l p t )
uniReg r1 r2 = compReg ( intReg cr1 cr2 )
where
cr1 = compReg r1
cr2 = compReg r2
di fReg r1 r2 = intReg r1 (compReg r2 )
symDReg r1 r2 = uniReg r1minr2 r2minr1
where
r1minr2 = di fReg r1 r2
r2minr1 = difReg r2 r1
Appendix C - Haskell Code 157
pip p ( LeafNode ch l ) = ( i s P o i n t i n H u l l p ch ) && (
isEven l )
pip p t = ( i s P o i n t i n H u l l p ch ) && ( a l l (== False ) lb )
&& ( isEven l )
where
ch = t h i s t
l = l e v e l t
c t s = c h i l d t r e e s t
lb = map ( pip p) c t s
p i r p r = (any (== True) l p i p )
where
l p i p = map ( pip p) r
−− Some A u x i l l a r y f u n c t i o n s f o r i n t e r s e c t i o n computation
−− Function to p a i r i n t e r s e c t i n g edges
p i e : : [ ( BPoint , BPoint ) ] −> [ ( BPoint , BPoint ) ] −> [ ( ( BPoint
, BPoint ) , ( BPoint , BPoint ) ) ]
p i e [ ] [ ] = [ ]
p i e e l 1 e l 2 = [ ( x , y ) | x <− e l1 , y <− e l2 , i s In tE x y ]
−− Given two edges , f u n c t i o n t e l l s whether the two edges
are i n t e r s e c t i n g or not
i s In tE : : ( BPoint , BPoint ) −> ( BPoint , BPoint ) −> Bool
i s In tE (a , b) ( c , d ) = t e s t
where −− t e s t i f the two edges i n t e r s e c t or not
t e s t = ( xor ( i sPo intRight ( a , b) c ) (
i sPo intRight ( a , b ) d) ) &&
( xor ( i sPo intRight ( c , d ) a ) (
i sPo intRight ( c , d ) b) )
−− xor f u n c t i o n
xor : : Bool −> Bool −> Bool
xor x y | x == True && y == False = True
| x == False && y == True = True
| otherwise = False
−− Given a l i s t o f i n t e r s e c t i n g edge pairs , f u n c t i o n
computes the i n t e r s e c t i o n p o i n t s
Appendix C - Haskell Code 158
f i p : : [ ( ( BPoint , BPoint ) , ( BPoint , BPoint ) )]−> [ BPoint ]
f i p [ ] = [ ]
f i p l ep = map (uncurry i n tPo in t ) l ep
−−Given two edges , f u n c t i o n computes the i n t e r s e c t i o n
p o i n t
i n tPo in t : : ( BPoint , BPoint ) −> ( BPoint , BPoint ) −> BPoint
in tPo in t ( a , b) ( c , d ) = i f ( isNeg ip ) then rp . dualPoint $
ip
else rp ip where
ip = crossProduct ( crossProduct a b) (
crossProduct c d)
dualPoint ( Bpt a b c ) = Bpt (−1 * a ) (−1 * b) (−1
* c )
isNeg ( Bpt a b c ) = ( a<0 && b<0 && c <0)
gd ( Bpt a b c ) = gcd a (gcd b c ) −− Given a
homogenous point , f i n d gcd o f c o o r d i n a t e s
rp ( Bpt a b c ) = ( Bpt ( a ‘ div ‘ g ) (b ‘ div ‘ g ) ( c ‘ div ‘
g ) )−− Given a point , reduce i t by gcd
where
g = gd ( Bpt a b c )
−−Function to c a l c u l a t e the c r o s s Produtct o f two
Homgenous Points
crossProduct : : BPoint −> BPoint−> BPoint
crossProduct ( Bpt w1 x1 y1 ) ( Bpt w2 x2 y2 ) = ( Bpt w3 x3 y3
)
where
w3 = ( x1 * y2 − x2 * y1 )
x3 = (w2 * y1 − w1 * y2 )
y3 = (w1 * x2 − w2 * x1 )
−−Given two h u l l s A and B, f u n c t i o n g i v e s the p o i n t s o f A
”in or on ” B
pointsInOrOn : : CHull −> CHull −> [ BPoint ]
pointsInOrOn ( CHull p1 ) ( CHull p2 ) = f i l t e r (\x −>i s P o i n t i n H u l l x ( CHull p2 ) ) p1
−− I s g i ven p o i n t p in g iven convex h u l l ch
i s P o i n t i n H u l l : : BPoint −> CHull −> Bool
Appendix C - Haskell Code 159
i s P o i n t i n H u l l p ( CHull ch ) = ( a l l (==True) pLef )
where
he l = mkEdges ch
pLef = map ( isPointLeftOrOn p) he l
−− Given a p o i n t and an edge , f u n c t i o n t e l l s whether p o i n t
i s ” l e f t or On” the edge
i sPointLeftOrOn : : BPoint −> ( BPoint , BPoint )−> Bool
i sPointLeftOrOn p ( a , b) = ( ( twiceOfArea p a b) <=0)
−−A u n i v e r s a l node
p1 = Bpt 1 10 10
p2 = Bpt 1 10 830
p3 = Bpt 1 1670 830
p4 = Bpt 1 1670 10
uniRing = Ring [ p1 , p2 , p3 , p4 ]
uniHul l = quickHul l uniRing
i s Un i : : Tree CHull −> Bool
i s Un i (Node ch ) = ( ( unChull ch \\ unChull un iHul l )
== [ ] )
i s Un i t = False
Appendix C - Haskell Code 160
module AHD.NDAHD where
import AHD. Algebra
import AHD.Geom
import AHD. AHDTree
import Data . List
import CnSimplex . CnSimplex
import AHD. ListAux
import AHD. Leve l
instance AHD [ CnSimplex ] Tree [ CnSimplex ] where
bu i ld mlev c l e v reg = i f ( dreg == [ ] ) | | ( mlev == c l e v
) then [ LeafNode chReg c l e v ]
else [ Node chReg ( concat (map ( bu i ld mlev ( inc
c l e v ) ) dreg ) ) c l e v ]
where
chReg = makePrim reg
dreg = sp l i t $ symDiff reg chReg
eva l mlev [ ( LeafNode s l ) ] = s
eva l mlev [ ( Node s l t s l ) ] = i f ( mlev == l ) then s
else symDiff s ( concat (map ( eva l mlev ) [ l t s ] ) )
instance InputPrim [ CnSimplex ] [ CnSimplex ] where
makePrim l s = convexHull (nub . concat .map cnSimp2simp
$ l s )
instance Input [ CnSimplex ] where
sp l i t = s p l i t J o i n t s . sp l i tOrdRegs
−− Symmetric d i f f e r e n c e o f x and y , i . e . ( x − y ) + ( y − x )
symDiff : : [ CnSimplex ] −> [ CnSimplex ] −> [ CnSimplex ]
symDiff x y = xSuby ++ ySubx
where
xSuby = deleteFirstsBy eqVs x y −− ( x − y )
ySubx = deleteFirstsBy eqVs y x −− ( y − x )
Appendix C - Haskell Code 161
−− S p l i t a l i s t o f n−s i m p l e x e s to a l i s t connected r e g i o n s
s p l i t R e g s : : [ CnSimplex ] −> [ [ CnSimplex ] ]
s p l i t R e g s [ ] = [ ]
s p l i t R e g s ( cs : c s s ) = ( cs : cntSimps ) : s p l i t R e g s css ’
where
cntSimps = connectedSimps [ cs ] c s s
css ’ = c s s \\ cntSimps
−− Extrac t the n−s i m p l e x e s o f ”cs2 ” which are connected to
at l e a s t one o f the
−− n−s i m p l e x e s o f a g iven l i s t ”cs1 ”
connectedSimps : : [ CnSimplex ] −> [ CnSimplex ] −> [ CnSimplex
]
connectedSimps cs1 cs2 = i f null neighSimps then [ ]
else neighSimps ++ ( connectedSimps neighSimps cs2 ’ )
where
neighSimps = f i l t e r ( isNeighSimp cs1 ) cs2
cs2 ’ = cs2 \\ neighSimps
−− Check i f a g i ven n−s imp lex i s a neighbour o f at l e a s t
ne o f the n−s i m p l e x e s
−− o f a g iven l i s t ”c s s ”
isNeighSimp : : [ CnSimplex ] −> CnSimplex −> Bool
isNeighSimp c s s ( vs2 , ) = not ( null [ ( vs1 , b) | ( vs1 , b )
<− css , length ( vs1 \\ vs2 ) == 1 ] )
−− Count the occurence o f each boundary (n−1)−s imp lex o f a
l i s t o f n−s i m p l e x e s ”cs ”
countOccur : : [ CnSimplex ] −> [ [ CnSimplex ] ]
countOccur = ( groupAllBy eqVs ) . concat .map boundary
−− Group a l i s t o f e lements by a g iven e q u a l i t y d e f i n i t i o n
groupAllBy : : (Eq a ) => ( a −> a −> Bool ) −> [ a ] −> [ [ a ] ]
groupAllBy [ ] = [ ]
groupAllBy eq ( x : xs ) = ( x : ys ) : groupAllBy eq zs
where
ys = [ t | t <− xs , eq t x ]
zs = xs \\ ys
−− Detect j o i n t s o f a l i s t o f n−s i m p l e x e s ”cs ” , i . e . the n
−s i m p l e x e s t h a t
Appendix C - Haskell Code 162
−− connect two or more r e g i o n s
j o i n t s : : [ CnSimplex ] −> [ CnSimplex ]
j o i n t s cs = map (head .nub) . ( f i l t e r ((> n) . length ) .
countOccur ) $ cs
where
n = cnSimpDim . head $ cs
−− Separate ord inary and j o i n t s n−s i m p l e x e s in a g iven
l i s t o f n−s i m p l e x e s ”cs ”
s epJo in t s : : [ CnSimplex ] −> ( [ CnSimplex ] , [ CnSimplex ] )
s epJo in t s cs = ( cs1 , cs2 )
where
cs2 = f i l t e r ( isASubSimp jntSimps ) cs
cs1= cs \\ cs2
jntSimps = j o i n t s cs
−− Check i f the v e r t e x e s o f an n−s imp lex ”cs ” i s s u b s e t o f
any n−s imp lex o f a
−− g iven l i s t ”c s s ”
isASubSimp : : [ CnSimplex ] −> CnSimplex −> Bool
isASubSimp c s s cs = not ( null . ( f i l t e r ( isSubSimp cs ) ) $
c s s )
−− Check i f the v e r t e x e s o f an n−s imp lex ( vs2 , b2 ) i s
s u b s e t o f the v e r t e x e s
−− o f another n−s imp lex ( vs1 , b1 )
isSubSimp : : CnSimplex −> CnSimplex −> Bool
isSubSimp ( vs1 , b1 ) ( vs2 , b2 ) = i sSubs e t vs2 vs1
−− S p l i t the ord inary reg ions , which i s the f i r s t e lement
o f the output pair ,
−− and the j o i n t n−s implexes , which i s the second element
o f the output p a i r
sp l i tOrdRegs : : [ CnSimplex ] −> ( [ [ CnSimplex ] ] , [ CnSimplex
] )
sp l i tOrdRegs c s s = ( s p l i t R e g s ordSimps , jntSimps )
where
( ordSimps , jntSimps ) = sepJo in t s c s s
−− S p l i t the ord inary n−s i m p l e x e s to connected r e g i o n s
s p l i t J o i n t s : : ( [ [ CnSimplex ] ] , [ CnSimplex ] ) −> [ [ CnSimplex
Appendix C - Haskell Code 163
] ]
s p l i t J o i n t s ( cs , [ ] ) = cs
s p l i t J o i n t s ( ( cs : c s s ) , jntSimps ) = ( cs ++ nghSimps ) :
s p l i t J o i n t s ( css , jntSimps ’ )
where
nghSimps = neighSimps cs jntSimps
jntSimps ’ = jntSimps \\ nghSimps
−− Extrac t the ne ighbour ing n−s i m p l e x e s o f a g iven one ( cs
) from a l i s t
−− o f n−s i m p l e x e s ( jntSimps )
neighSimps cs jntSimps = f i l t e r ( isNeighSimp cs ) jntSimps
Appendix C - Haskell Code 164
module Str ipTree where
import Prelude hiding ( toInteger )
import AHD.Geom
import AHD. Algebra
import Data . List
import ConvexHull
import Language . ConvsAlg
import Data . Function
import AHD. ListAux
import AHD. AHDTree
instance Input Curve where
i s S p l i t t a b l e curve = length curve > 2
sp l i t c = part s where
i l i n e = g e t l i n e . makePrim $ c
eps = getxtrmpoints i l i n e c
breakpo int = f s t (maximumBy (compare ‘ on ‘ snd ) eps
)
par t s = s p l i t L i s t breakpo int c
merge l c = nub $ l c
instance AHD Curve Tree S t r i p where
bu i ld mlev c l e v c = i f ( mlev == c l e v ) | | (
i s S p l i t t a b l e c == False )
then [ ( LeafNode p c l e v ) ]
else [ Node p ( concat (map ( bu i ld mlev ( inc
c l e v ) ) par t s ) ) c l e v ]
where
p = makePrim c
par t s = sp l i t c
eva l mlev [ ( LeafNode s l ) ] = getCurve s
eva l mlev [ ( Node s l t s l ) ] = i f ( mlev == l ) then (
getCurve s )
else ( merge ( getCurve s ) ( concat l l r ) ) where
l l r = (map ( eva l mlev ) ( toLL l t s ) )
Appendix C - Haskell Code 165
instance InputPrim Curve S t r i p where
makePrim ( p1 : p2 : [ ] ) = S t r i p { spx = p1x , spy = p1y ,
epx = p2x , epy = p2y , wl = 0 . 0 , wr = 0.0}where
( 1 . 0 , p1x , p1y ) = h2c ( unPoint p1 )
( 1 . 0 , p2x , p2y ) = h2c ( unPoint p2 )
makePrim lp = St r i p { spx = p1x , spy = p1y , epx = p2x ,
epy = p2y , wl = wll , wr = wrr}where
( 1 . 0 , p1x , p1y ) = h2c ( unPoint p1 )
( 1 . 0 , p2x , p2y ) = h2c ( unPoint p2 )
( p1 , p2 ) = g e t i n i t l i n e lp
wrr = snd . head $ l x t
w l l = snd . last $ l x t
l x t = getxtrmpoints ( p1 , p2 ) lp
instance Prim St r i p where
mergeComp c cs = cs
−− Some Auxi lary f u n c t i o n s
−− g iven a curve f u n c t i o n r e t u r n s the i n i t i a l l i n e
g e t i n i t l i n e : : Curve −> ( BPoint , BPoint )
g e t i n i t l i n e ( x : [ ] ) = error ” s i n g l e element in the l i s t ”
g e t i n i t l i n e ( x : xs ) = (x , last xs )
−− g iven a l i n e and points , the f u n c t i o n r e t u r n s l e f t m o s t
p o i n t and r i g h t m o s t p o i n t wi th wl and wr
getxtrmpoints : : ( BPoint , BPoint ) −> Curve −> [ ( BPoint ,
Double) ]
getxtrmpoints ( sp , ep ) cp = i f ( po l == [ ] ) then [ ( ( Bpt 0 0
0) , 0) , ( fpr , wr ) ]
else i f ( po l == [ ] ) then [ ( f p l , wl ) , ( ( Bpt 0 0 0) , 0) ]
else [ ( f p l , wl ) , ( fpr , wr ) ] where
por = f i l t e r ( i sPo intRight ( ep , sp ) ) cp
pol = f i l t e r ( i sPo intRight ( sp , ep ) ) cp
f p l = fa r thPo in t ( sp , ep ) po l
Appendix C - Haskell Code 166
fp r = fa r thPo in t ( sp , ep ) por
wl = t r i h e i g h t ( sp , ep , f p l )
wr = t r i h e i g h t ( sp , ep , f p r )
−− g iven a t r i a n g l e , the f u n c t i o n r e t u r n s the h e i g h t o f
the t r i a n g l e
t r i h e i g h t : : Tr iang l e −> Double
t r i h e i g h t ( p1 , p2 , p3 ) = he ight where
area = toDouble (abs $ twiceOfArea p1 p2 p3 ) /2
he ight = ( area * 2) / base
base = sqrt ( ( ( abs ( p1x − p2x ) ) ˆ 2) + ( ( abs ( p1y − p2y )
) ˆ 2) )
(1 , p1x , p1y ) = h2c ( unPoint p1 )
(1 , p2x , p2y ) = h2c ( unPoint p2 )
(1 , p3x , p3y ) = h2c ( unPoint p3 )
−−g iven a s t r i p , g e t the l i n e
g e t l i n e ( S t r i p p1x p1y p2x p2y w1 w2) = ( c2h ( p1x , p1y ) ,
c2h ( p2x , p2y ) )
getCurve ( S t r i p p1x p1y p2x p2y w1 w2) = c2h ( p1x , p1y ) : [
c2h ( p2x , p2y ) ] : : Curve
c2h (x , y ) = ( Bpt 1 ( toInteger x ) ( toInteger y ) ) : : BPoint
Appendix C - Haskell Code 167
module AHD. ListAux where
import Data .Maybe ( fromJust )
import Data . List (elemIndex )
−− f u n c t i o n t h a t f l a t e n s a g iven l i s t o f t u p l e s
f l a t t e n : : [ ( a , a ) ] −> [ a ]
f l a t t e n [ ] = [ ]
f l a t t e n ( x : xs ) = f s t x : snd x : f l a t t e n xs
−− make CCW edges out o f a g iven Hul l o f p o i n t s
mkEdges : : [ a ] −> [ ( a , a ) ]
mkEdges [ ] = [ ]
mkEdges ( x : [ ] ) = [ ( x , x ) ]
mkEdges ( x : y : [ ] ) = [ ( x , y ) , (y , x ) ]
mkEdges ( x : xs ) = pa i r x xs ++ [ ( last xs , x ) ]
pa i r : : a −> [ a ] −> [ ( a , a ) ]
pa i r a [ ] = [ ]
pa i r a ( x : xs ) = ( a , x ) : ( pa i r x xs )
f s t 3 ( a , b , c ) = a
snd3 ( a , b , c ) = b
trd3 (a , b , c ) = c
f s t 4 ( a , b , c , d ) = a
snd4 ( a , b , c , d ) = b
trd4 (a , b , c , d ) = c
f r t 4 ( a , b , c , d ) = d
toLL : : [ a ] −> [ [ a ] ]
toLL [ ] = [ ]
toLL ( x : xs ) = [ x ] : toLL xs
−− Given a l i s t [ a . . z ] , f u n c t i o n s p l i t s the g i ven l i s t a t
g i ven p o i n t x i n t o two p a r t s i . e . [ a . . x ] , [ x . . z ]
s p l i t L i s t : : Eq a => a−> [ a ] −> [ [ a ] ]
Appendix C - Haskell Code 168
s p l i t L i s t p l = d1 : [ p : d2 ]
where
( d1 , d2 ) = splitAt iop l
iop = ( fromJust (elemIndex p l ) )+1
Biography of the Author
Rizwan Bulbul was born in Gilgit, Pakistan on February 10, 1977. He received
his initial education from Public School and College, Gilgit and college education
from Government College Lahore, Pakistan. He did his BSc. (Hons.) in Com-
puter Sciences from International Islamic University Islamabad, Pakistan in 2001.
He then did MS in Computer Sciences from National University of Computer and
Emerging Sciences, FAST NU, Islamabad, Pakistan in 2005.
He started his professional career as a research associate at Center for Agroin-
formatics Research, CAIR, at FAST NU. He then joined Department of Com-
puter Sciences, Karakorum International University, Gilgit, Pakistan as Assistant
Professor. Since November 2007, he has been a researcher and PhD candidate in
the Geoinformation group of Institute for Geoinformation and Cartography at the
Vienna University of Technology. He is interested in computational geometry for
GIS, functional programming (Haskell), data structures and algorithms, dimen-
sion independent geometric modeling and city modeling. He was involved with
COST Action TU801 (Semantic Enrichment of 3D City Models for Sustainable
Urban Development) as a student member. He was investigating the geometric
issues in city models.
Rizwan has authored and co-authored several research papers published in
conference proceedings. He is a candidate for the “Doktor der technischen Wis-
senschaften” degree from the Vienna University of Technology in 2011.
169
top related