Lecture 3D Geometry Processing
Delaunay Mesh Generation
Prof. Dr. David Bommes
Computer Graphics Group

Last week: Poisson Surface Reconstruction

  • Summary

Last week: Implementation

  • Given the input point set:
    • Setup adaptive octree
    • Compute vector field
    • Compute indicator function
    • Extract iso-surface

Last week: Implementation

  • Given the input point set:
    • Setup adaptive octree
    • Compute vector field
    • Compute indicator function
    • Extract iso-surface

Last week: Implementation

  • Given the input point set:
    • Setup adaptive octree
    • Compute vector field
      • Define function basis
      • Splat the samples
    • Compute indicator function
    • Extract iso-surface

Last week: Implementation

  • Given the input point set:
    • Setup adaptive octree
    • Compute vector field
      • Define function basis
      • Splat the samples
    • Compute indicator function
    • Extract iso-surface

Last week: Implementation

  • Given the input point set:
    • Setup adaptive octree
    • Compute vector field
      • Define function basis
      • Splat the samples
    • Compute indicator function
    • Extract iso-surface

Last week: Implementation

  • Given the input point set:
    • Setup adaptive octree
    • Compute vector field
      • Define function basis
      • Splat the samples
    • Compute indicator function
    • Extract iso-surface

Last week: Implementation

  • Given the input point set:
    • Setup adaptive octree
    • Compute vector field
    • Compute indicator function
      • Compute divergence
      • Solve Poisson equation
    • Extract iso-surface

Last week: Implementation

  • Given the input point set:
    • Setup adaptive octree
    • Compute vector field
    • Compute indicator function
      • Compute divergence
      • Solve Poisson equation
    • Extract iso-surface

Last week: Implementation

  • Given the input point set:
    • Setup adaptive octree
    • Compute vector field
    • Compute indicator function
    • Extract iso-surface

Outline

  1. Voronoi Diagrams
  2. Delaunay Triangulation
  3. Mesh Generation

Voronoi Diagrams

Definition

  • Given a set of 2D sample points \(\{\mathbf{p}_1, \ldots , \mathbf{p}_n\}\)
  • Partition the plane by assigning each 2D point \(\mathbf{x}\) to its nearest sample.
  • All points assigned to \(\mathbf{p}_i\) form its Voronoi cell \[\mathcal{V}(\mathbf{p}_i) = \left\{ \mathbf{x} \in \mathbb{R}^2 : ||\mathbf{x}-\mathbf{p}_i|| \leq ||\mathbf{x}-\mathbf{p}_j|| \quad \forall j \neq i \right\}\]
  • Edges and vertices of these cells form the Voronoi diagram (VD).

Georgy Voronoi
1868 - 1908

Voronoi Diagram

Voronoi Diagram

  • Each point pair \(\mathbf{p}_i\), \(\mathbf{p}_j\) defines a perpendicular bisector
  • The bisector of \(\mathbf{p}_i\) and \(\mathbf{p}_j\) defines two half-planes
  • Let \(H(\mathbf{p}_i, \mathbf{p}_j)\) be the half-plane containing \(\mathbf{p}_i\)
  • Voronoi cells are intersections of half-planes \[\mathcal{V}(\mathbf{p}_i)=\cap_{i\neq j} H(\mathbf{p}_i,\mathbf{p}_j)\]

Voronoi Diagram

Voronoi Diagram

Complexity of Voronoi Diagram

  • Exactly one Voronoi cell per point, but cells are defined by intersecting \(n\) half-planes \[\mathcal{V}(\mathbf{p}_i)=\cap_{i\neq j} H(\mathbf{p}_i,\mathbf{p}_j)\]
  • Cells could have \(O(n)\) edges
    • VD could have \(O(n^2)\) edges
    • Or is it just \(O(n)\) edges?

Voronoi Diagram

Complexity of Voronoi Diagram

  • Dual graph of VD is a triangle mesh (why?)
    • Dual: vertex \(\rightarrow\) face, edge \(\rightarrow\) edge, face \(\rightarrow\) vertex
    • Euler formula for triangle meshes: \(E \approx 3V = O(n)\)
    • Dual mesh has \(O(n)\) edges
    • Primal mesh has \(O(n)\) edges

Properties of Voronoi Diagram

  • Voronoi cells are convex
  • Voronoi vertices have valence 3
    (if no 4 points are co-circular)
  • Voronoi vertices are circumcenters of its three defining sample points
  • These circumcircles do not contain other sample points

Voronoi Algorithms

  • Incremental construction: \(O(n^2)\)
    • Insert point by point and update VD
    • New point falls into circumcircles of some points
    • Only those regions have to be updated
  • Fortune’s sweep-line algorithm: \(O(n \log n)\)
    • Sort samples along x-direction
    • Sweep plane from left to right
    • Construct VD behind plane
    • Optimal complexity

Voronoi Visualization

  • Draw a cone in z-direction at each sample \(\mathbf{p}_i\)
  • z-values measure distance from samples, i.e. \(z_i(\mathbf{x}) = dist(\mathbf{p}_i,\mathbf{x})\)
  • View cones from below (parallel projection)
  • Cone intersections project to Voronoi edges!

3D cones
side view
bottom view
perspective
bottom view
parallel

Outline

  1. Voronoi Diagrams
  2. Delaunay Triangulation
  3. Mesh Generation

Delaunay Triangulation

Delaunay Triangulation

  • The dual graph of the Voronoi diagram is a planar straight line graph, the
    Delaunay triangulation (DT)

Boris Delaunay
1890 - 1980

  • Circumcircles of DT triangles are empty
    • DT triangles are duals of VD vertices
    • Criterion can be used for DT construction

Delaunay Triangulation

Edge Flipping

  • Check whether an edge is Delaunay by testing circumcircles of its incident triangles

  • Flip edge to make it Delaunay

In-Circle Test

  • How to efficiently compute \(InCircle(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})\)?
    - Lift onto paraboloid: \((x,y) \rightarrow (x, y, x^2 + y^2)\)
    - Lifted points \(A'\), \(B'\), \(C'\) define a plane cutting through
    the paraboloid
    - 3D intersection curve projects to 2D circumcircle

In-Circle Test

  • How to efficiently compute \(InCircle(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})\)?
    - \(\mathbf{D}\) is in/out circumcircle \(\Leftrightarrow\) \(\mathbf{D}'\) is below/above plane
    - \(\mathbf{D}'\) below/above \(\Leftrightarrow\) \(volume(\mathbf{A}',\mathbf{B}',\mathbf{C}',\mathbf{D}') < 0\) or $ > 0$

\[InCircle(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}) \Leftrightarrow -\det \left[ \begin{array}{cccc} A_x & A_y & A_x^2 + A_y^2 & 1 \\ B_x & B_y & B_x^2 + B_y^2 & 1 \\ C_x & C_y & C_x^2 + C_y^2 & 1 \\ D_x & D_y & D_x^2 + D_y^2 & 1 \end{array} \right] < 0 \]

In-Circle Test

  • How to efficiently compute \(InCircle(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})\)?

\[InCircle(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}) \Leftrightarrow - \det \left[ \begin{array}{cccc} A_x & A_y & A_x^2 + A_y^2 & 1 \\ B_x & B_y & B_x^2 + B_y^2 & 1 \\ C_x & C_y & C_x^2 + C_y^2 & 1 \\ D_x & D_y & D_x^2 + D_y^2 & 1 \end{array} \right] < 0 \]

  • Swapping rows in determinants
    • \(InCircle(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}) == InCircle(\mathbf{C}, \mathbf{D}, \mathbf{A}, \mathbf{B})\)
    • \(InCircle(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}) == -InCircle(\mathbf{B}, \mathbf{C}, \mathbf{D}, \mathbf{A})\)

Maximum Minimum Angle

  • Delaunay criterion maximizes the minimum angle
    • Can be seen from Thales’ theorem

  • Delaunay triangulation avoids small angles
    • Leads to numerically preferable meshes
    • Most important advantage of DT

Incremental Algorithm

  • Given a set of 2D sample points \(\{ \mathbf{p}_1, \mathbf{p}_2, \mathbf{p}_3 \}\)
    • Insert points one by one
    • Flip edges to establish Delaunay property
  • Avoid special cases
    • Add a big triangle \((\mathbf{q}_1, \mathbf{q}_2, \mathbf{q}_3)\) containing all \(\mathbf{p}_i\)
    • Those points can later be removed
    • Points always inserted into existing triangles

Incremental Algorithm

Incremental Algorithm

  • For each point \(\mathbf{p}_i\)
    • Find containing triangle
    • Insert point into triangle (1-to-3 split)
    • Flip edges to re-establish Delaunay property
      (see [deBerg] on which edges to check)
  • Incremental algorithm
    • Complexity is \(O(n^2)\), but it’s simple to implement
    • \(O(n \log n)\) with better triangle search
  • Fortune’s sweep-line algorithm
    • Optimal \(O(n \log n)\) complexity, but complicated

Outline

  • Voronoi Diagrams
  • Delaunay Triangulation
  • Mesh Generation

Mesh Generation

Constrained Delaunay Triangulation

  • Enforce certain edges in triangulation
    • Either prevent flipping (\(\rightarrow\) bad triangles)
    • Or subdivide edges sufficiently (\(\rightarrow\) many triangles)

2D Meshing

Input points & edges

Constrained Delaunay triangulation

Remove exterior

Images from J. Shewchuck

2D Meshing

  • 2D Delaunay triangulation
    • Maximizes minimum angle
    • Optimal triangulation for given set of vertices
    • Why can there still be bad triangles?

Delaunay Refinement

  • Delaunay triangulation might contain bad triangles, depending on vertex distribution
  • Refine triangulation
    • Insert new vertices, eliminate bad triangles
  • Measure triangle quality by
    • circumradius / shortest-edge
    • smallest inner angles

Delaunay Refinement

  • Insert new vertices to eliminate bad triangles
  • Eliminate “bad” triangles by inserting their circumcenter into the triangulation
  • Bad triangle will fail the empty-circumcircle test and will therefore be removed

Delaunay Refinement

  • Insert new vertices to eliminate bad triangles
  • Eliminate “bad” triangles by inserting their circumcenter into the triangulation
  • Bad triangle will fail the empty-circumcircle test and will therefore be removed

2D Meshing

Images from J. Shewchuck

Delaunay Refinement

Images from J. Shewchuck

Centroidal Voronoi Diagrams

  • How can we get a more regular triangulation?
    • Delaunay triangulation of Centroidal Voronoi Diagram (CVD)
  • Definition:
    • All points are centroids of their Voronoi cells

Centroidal Voronoi Diagrams

Centroidal Voronoi Diagrams

Centroidal Voronoi Diagrams

Literature

  • Botsch et al., Polygon Mesh Processing, AK Peters, 2010.
    • Chapter 6.4

  • de Berg et al., Computational Geometry: Algorithms and Applications, Springer Verlag, 2008.

  • O’Rourke, Computational Geometry in C, Cambridge University Press, 1998.