Research Context

• (Blue sky:) Geologically correct terrain representation
  • Soon to be DARPA/DSO/Geo* sponsored (via NGA).
  • Replace Fourier series, fractals with scooping and other nonlinear operators that respect how terrain is formed.
• (Applications) to geometry and cartography
  • e.g., multiple observer siting with intervisibility on terrain
• Efficient fundamental operations, like finding nearest neighboring points, finding connected components, etc.
  • A proper foundation is necessary for the blue sky projects' success.

Today: How to process large geometric databases quickly!

What?

1. Larger and larger geometric databases are available, with tens of millions of primitive components.
2. We want to perform geometric operations, such as
   1. interference detection
   2. boolean operations: intersection, union
   3. planar graph overlay
   4. mass property computation of the results of some boolean operation
3. Apps:
   1. Find the volume of an object defined as the union of many overlapping primitives.
   2. Two object interfere iff the volume of their union is positive.

Themes

1. Limited address space => must minimize storage => I/O is faster.
2. For N>>1,000,000 lg N is nontrivial => deprecate binary trees
3. Minimize explicit topology, expecially in 3D.
4. Plan for 3D; many 2D data structures do not easily extend to 3D.
   1. 2D line sweep easily finds intersecting edge segments. 3D plane sweep would need to maintain a planar graph of the active faces. Messy!
   2. 2D Voronoi diagram is useful. However space for 3D Voronoi diagram is easily O(N²).

Confessions

1. I'm not necessarily a deep philosophical thinker because I'm always seeing holes in the generalities. I like Galileo better than Aristotle. Galileo experimented.
2. I try to do small things well, to lay a foundation of computational observations and tools on which to generalize.
3. This is driven by a deep knowledge of classical geometry, where the order is implicit in the axioms. You don't need to store it explicitly to use it.
4. Example of hidden order: the centroid, circumcenter, and orthocenter of a triangle are collinear.
5. Yes, the above points are a little contradictory.

Theme: Minimum Explicit Topology

1. What explicit info does the application need? Less > simpler
2. Object: polygon with multiple nested components and holes.

   

3. Apps:
   1. area
   2. inclusion testing.
4. Complete topology: loops of edges; the tree of component containments.
5. Necessary info: the set of oriented edges.

Point Inclusion Testing on a Set of Edges

1. "Jordan curve" method
2. Extend a semi-infinite ray.
3. Count intersections.
4. Odd <==> Inside
5. Obvious but bad alternative: summing the subtended angles. Implementing it w/o actually doing arctan to compute the angles, and handling special cases wrapping around 2pi is rather tricky. Then it reduces to the Jordan curve method.

Area Computation on a Set of Edges

1. Each edge, with the origin, defines a triangle.
2. Sum their signed areas.

Advantages of the Set of Edges Data Structure

1. It's simple enough to debug.

   SW can be simple enough that there are obviously no errors, or complex enough that there are no obvious errors.

2. It takes less space to store.
3. Parallelizing is easy.
   1. Partition the edges among the processors.
   2. Each processor counts intersections or sums areas independently, to produce one subtotal.
   3. Total the subtotals.

What About a Set of Vertices Data Structure?

1. Too simple.
2. Ambiguous: two distinct polygons may have the same set of edges.

Set of Vertex-Edge Incidences

Here's another minimal data structure.

The only data type is the incidence of an edge and a vertex, and its neighborhood. For each such:

1. V = coord of vertex
2. T = unit tangent vector along the edge
3. N = unit vector normal to T pointing into the polygon.

Polygon: {(V, T, N)} (2 tuples per vertex)

Perimeter = -(V.T).            
Area = 1/2 ( V.T V.N )

Multiple nested components ok.

Demonstration: Mass Properties of the Union of Millions of Cubes

Demonstration: Mass Properties of the Union of Millions of Cubes

1. This illustrates the above ideas.
2. This is a prototype on an easy subcase (congruent axis-aligned cubes).
3. The ideas extend to general polyhedra.
4. This is not a statistical sampling method - the output is exact, apart from losing 15 significant digits during the computation.
5. This is not a subdivision-into-voxel method - the cubes' coordinates can be any representable numbers.

Application: Cutting Tool Path

1. Represent path of a tool as piecewise line.
2. Each piece sweeps a polyhedron.
3. Volume of material removed is (approx) volume of union of those polyhedra.
4. Image is from Surfware Inc's Surfcam website.

Traditional N-Polygon Union

1. Construct pairwise unions of primitives.     Iterate.
2. Time depends on assumed size of intermediates, and assumed elementary intersection time.
3. Let P = size of union of an M-gon and an N-gon. Then P=O(MN).
4. Time for union (using line sweep) T=(P lg P) .
5. Total T=O(N²lg N).
6. Hard to parallelize upper levels of computation tree.

Problems With Traditional Method

1. lg N levels in computation tree cause lg N factor in execution time. If N>20, that's significant.
2. Intermediate swell: gets worse as overlap is worse. Intermediate computations may be much larger than the final result.

3. 

   The explicit volume has a complicated topology:
   1. loops of edges,
   2. shells of faces,
   3. nonmanifold adjacancies.
      
4. This is tricky to get right.
5. The explicit volume is not needed for computing mass properties.
6. The set of the vertices with their neighborhoods suffices.

Volume Determination

1. Box:
   1. V = _i s_i x_i y_i z_i
   2. s_i: +1 or -1
   
   
2. General rectilinear polygons:
   1. 8 types of vertices, based on neighborhood
   2. 4 are type "+", 4 "-"
   3. Area = (s_ix_iy_i)
   4. Perimeter = (+x_i+y_i)

Rectilinear Polyhedron Volume

1. V = _i s_i x_i y_i z_i
2. 2 classes of vertices: s _i = +1 , or -1.
3. Extends to lower dimensional properties, like area and edge length.
4. Extends to general polyhedra.

Vertex Classes

1. Each vertex has 8 adjacent octants.
2. Each octant is either interior or exterior to the polyhedron.
3. Which octants are interior determines the vertex's sign.
4. Debugging this was tricky.

Properties

The output union polyhedron will be represented as a set of vertices with their neighborhoods.

1. no explicit edges; no edge loops
2. no explicit faces; no face shells
3. no component containment info
4. represents general polygons: multiple nested or separate comps
5. any mass property determinable in one pass thru the set
6. parallelizable
7. amenable to external storage implementation, if necessary.

Volume Computation Overview

1. Find all the vertices of the output object.
2. For each vertex, find its location and local geometry.
3. Sum over the vertices, applying the formula.

Finding the Vertices

3 types of output vertex:

1. vertex of input cube
2. intersection of edge of one cube with face of another,
3. intersection of three faces of different cubes.

Find possible output vertices, and filter them

An output vertex must not be contained in any input cube.

Isn't intersecting all triples of faces, then testing each candidate output vertex against every input cube too slow? No, if we do it right.

Use a 3D Grid - Summary

1. Overlay a uniform 3D grid on the universe.
2. For each input primitive (cube, face, edge), find which cells it overlaps.
3. With each cell, store the set of overlapping primitives.

Grid - Properties:

1. Simple, sparse, uses little memory if well programmed.
2. Parallelizable.
3. Robust against moderate data nonuniformities.
4. Bad worst-case performance: could be defeated be extremely nonuniform data.
5. Ditto any hierarchical method like an octree.

Advantage of the Grid

1. Intersecting primitives must occupy the same cell.
2. The grid filters the set of possible intersections.

Covered Cell Concept

1. Only visible intersections contribute to the output.
2. That's often a small fraction - inefficient.
3. Solution: add the cubes themselves to the grid.

Adding the Cubes Themselves to the Grid

1. For each cubes, find which cells it completely covers.
2. If a cell is completely covered by a cube, then nothing in that cube can contribute to the output.
3. Find the covered cells first.
4. Do not insert any objects into the covered cells.
5. Intersect pairs and triples of objects in the noncovered cells.

If the cell size is somewhat smaller than the edge size then almost all the hidden intersections are never found. Good.

Expected time = (size(input) + size(useful intersections)).

Face-Face Intersection Details

1. Iterate over the grid cells.
2. In each cell, test all triples of faces, where each face must be from a different cube, for intersection.
3. Three faces intersect if their planes intersect, and the intersection is inside each face (2D point containment).
4. If the 3 faces intersect, look up s_i in a table and add a term to the accumulating volume.
5. The implementation is easier for cubes.

Point Containment Testing

P is a possible vertex of the output union polyhedron.

Q: Is point P contained in any input cube?

A:

1. Find which cell, C, contains P.
2. If C is completely covered by some cube then P is inside the covering cube.
3. Otherwise, test P against all the cubes that overlap C.

The expected number of such cubes is constant, under broad conditions.

Simulation of Simplicity

1. We must resolve degeneracies.
2. For any equality, break the tie by comparing the cubes' ID numbers.
3. However, now lower-dimensional mass properties give different results than for a regularized set operation.
4. E.g., if two cubes abut on a face, then either both faces or neither face are included.
5. That is topologically valid.

Execution Time

1. N: number of cubes
2. L: edge length, in a 1x1x1 universe.
3. Expected number of 3-face intersections = (N³ L⁶)

Effect of Grid

1. Good choice of cell size: half the cube size.
2. G: number of grid cells on a side = 2/L.
3. Number of face triples: N³
4. Prob. of a 3-face test succeeding = N^-2L⁶.
5. Depending on asymptotic behavior of L(N), this tends to 0.
6. With the grid, prob. of 3 tested faces actually intersecting = c, independent of N and L(N) .
7. Big improvement!

Effect of Covered Cells

1. Expected number of 3-face intersections = (N³L⁶)
2. However, for uniform i.i.d. input, expected number of visible intersections = (N).
3. With covered cells, probability that a computed intersection that is visible = c, independent of N and L(N) .
4. Big improvement!
5. Time to test if a point is inside any cube is also constant.
6. Total time reduces to (N) .

Parallel Implementation

1. Pentium-class processors allow only 3GiB memory per user process.
2. Machine may have 12GiB real memory, so real memory > virtual memory!
3. Solution: decompose problem into pieces and process in parallel.
4. Cannot just slice up the input: this increases the area and edge length.
5. Observation: inserting objects into cells is quicker than testing the cells for intersections.
6. Solution:
   1. Insert the {cubes, faces, edges} into the cells.
   2. Dustribute the cells among the processors, and test for intersections in parallel.
   3. Each process reads several GiB of data, but writes only 3 words: its part of the volume, area, length.

Implementation

1. Very compact data structures.
2. Random (Tausworth generator) uniform i.i.d. cubes.
3. 1000 executable lines of C++.
4. dual 2.4GHz Xeon with 4GiB memory, Linux.

Small Datasets are Fast

1. 10,000 cubes; 0.03 edge length
2. 60x60x60 grid
3. Not using parallelism
4. Execution time: 2.5 CPU secs
5. Memory: 15 MB.

Large Datasets are Feasible

1. 20,000,000 cubes
2. 0.009 edge length (quite small, few intersections).
3. 500x500x500 grid
4. Using 4 parallel processes.
5. Execution time: 2000 to 2200 CPU secs per process. I.e., the work was evenly divided.
6. Memory per process: 2.6 GB

Why not larger datasets? Address space limitations. Creating more parallel processes thrashes the memory. We're creating separate sequential processes, but still working on the interfacing.

Implementation Validation

1.  
   1. The terms summed for the volume are large and mostly cancel.
   2. Errors would be unlikely to total to a number in [0,1].
2.  
   1. The expected volume is 1-(1-L³)^N.
   2. Compare to computed volume.
   3. Assume that coincidental equivalence is unlikely.
3.  
   1. Construct specific, maybe degenerate, examples with known volume.

Extension to General Polygons

The only data type is the incidence of an edge and a vertex, and its neighborhood. For each such:

1. V = coord of vertex
2. T = unit tangent vector along the edge
3. N = unit vector normal to T pointing into the polygon.

Polygon: {(V, T, N)} (2 tuples per vertex)

Perimeter = -(V.T).            
Area = 1/2 ( V.T V.N )

Multiple nested components are ok.

Extension to 3D

The only data type is the incidence of a face, an edge and a vertex, and its neighborhood. For each such:

1. V = coord of vertex
2. T = unit tangent vector along the edge
3. N = unit vector normal to T in the plane of the face, pointing into the face.
4. B = unit vector normal to T and N, pointing into the polyhedron.

Polyhedron: {(V, T, N, B)}

Perimeter = -1/2 (V.T)          
Area = 1/2 ( V.T V.N )
Volume = -1/6 ( V.T V.N V.B )

Intersection Volumes from Overlaying 3-D Triangulations

Existing 2-D implementation: overlay two planar graphs; find areas of all intersection regions.

E.g. overlay counties of the coterminous US against hydrography polygons. Map stats:

Map #0: 55068 vertices, 46116 edges, 2985 polygons
Map #1: 76215 vertices, 69835 edges, 2075 polygons

Total time (excluding I/O): 1.58 CPU seconds.

Future: Extend to 3-D.

Application: Interpolate some property from one triangulation to the other.

Extension: More Complicated Boolean Operations

Given a red object defined as the union of many overlapping cubes, and a similar blue one.

1. Do they intersect?
2. Simplified to 2-D: N red edge segments, and N blue ones.
3. There may be O(N²) red-red and blue-blue intersections.
4. You don't know them.
5. There are O(1) red-blue intersections.
6. Find them.
7. No subquadratic worst-case algorithm (SFAIK).

Moral: sometimes heuristics are necessary.

Other Possible Extensions

1. Online algorithm, with incremental input insertions, would be easy.
2. Online deletion of cubes should be only a little harder, if covered cells not used.
3. This could be fast enough to be process hundreds of cubes at video update rates (1/30 sec).

Summary

Guiding principles:

1. Use minimal possible topology, and compact data structures.
2. Short circuit the evaluation of volume(union(cubes)).
3. Design for expected, not worst, case input.
4. External data structures unnecessary, tho possible.

Allows very large datasets to be processed quickly in 3-D with algorithms where sweep lines might have problems.

A paper on this subject is at http://wrfranklin.org/research/union3/union3.pdf.

The html version of this talk is on my website http://wrfranklin.org/Research/geo_ops_millions-jul2004/ . It is designed to be displayed in full screen mode with Opera.

Thanks!

Portions of this material are based upon work supported by the National Science Foundation under Grants ENG-7908139, ECS-8021504, ECS-8351942, CCF-9102553, and CCF-0306502, and by IBM, Sun Microsystems, and Schlumberger--Doll Research.