2  Planar Triangulations

“I will be glad if I have succeeded in impressing the idea that it is not only pleasant to read at times the works of the old mathematical authors, but this may occasionally be of use for the actual advancement of science.”
– Constantin Carathéodory (speaking at a meeting of the Mathematical Association of America in 1936)

We will discuss the notion of the convex hull of a finite subset of the Euclidean plane, and about ways of “triangulating” such a convex hull.

2.1 Convexity

We discuss very briefly the notion of convexity, which will often come up in this section and later.

Convexity

1. For \bar{x}, \bar{y} \in \mathbb{R}^n we can parametrize the line segment joining them as \{t \bar{x} +(1-t) \bar{y}: t \in [0,1]\}
2. A subset C \subset \mathbb{R}^n is said to be convex iff for each \bar{x}, \bar{y} \in C the line segment [\bar{x}, \bar{y}] \subset C
3. The convex hull of a subset S \subset \mathbb{R}^n is the smallest convex subset of \mathbb{R}^n that contains S

Exercise 2.1  

1. Show that the convex subsets of \mathbb{R} are the intervals.
2. Show that in \mathbb{R}^n the open balls and closed balls are convex
3. Show that the intersection of convex sets is convex.

2.2 The convex hull of a finite subset of \mathbb{R}^2

We will now look at the convex hull of a finite subset S=\{s_1, s_2, \cdots, s_k\} \subset \mathbb{R}^2.

Theorem:

Given S=\{s_1, s_2, \cdots, s_k\} \subset \mathbb{R}^2, the convex hull of S equals \{\alpha_1 s_1 + \cdots + \alpha_k s_k: \alpha_i \ge 0, \sum \alpha_i =1\}

Exercise 2.2  

1. Give an augument explaining why the convex hull of a finite point set in \mathbb{R}^2 is a convex polygon.
2. Observe that a (filled in) nondegenerate triangle is the convex hull of its three vertices. In fact, if \bar{p},\bar{q},\bar{r} \in \mathbb{R}^2 are non-collinear, then T=\textbf{conv}\{\bar{p},\bar{q},\bar{r}\} = \{\alpha_1 \bar{p}+\alpha_2 \bar{q} + \alpha_3 \bar{r}: \alpha_i \ge 0, \sum \alpha_1=1\} (The \alpha_i in the triangle example are called the barycentric coordinates. We will use this idea extensively later).
3. Show that each point in the convex hull of a finite subset S of \mathbb{R}^2 is contained in the convex hull of some triple of points in S (This is Carathéodory’s theorem in dimension 2)

Constantin Carathéodory was a mathematician with contributions to several areas of mathematics, including convex geometry.

Convex Hull: The Gift Wrapping and Graham Scan Algorithms

1. Gift Wrapping (Chand and Kapur, 1970): Given a planar point set is general position. We can choose a distinguished point (say, top left corner - or bottom right corner) and consider it the base point. From the base point, we connect all other points by a line segment, then choose the next point of the wrapping as the point that makes the largest angle with the horizontal line. Considering the new point to be the base, repeat the procedure. The computational complexity of the algorithm is \mathcal{O}(nm), where n is the number of points and m is the number of points on the hull.

2. Graham Scan (Graham, 1972): This is a modification of the Gift Wrapping Algorithm where - once an initial base point is selected - we sort the remaining points according to the angle they form with the horizontal. Afterwords, we process the points according to the magnitude of the angles. We may need to discard some points from the hull based on whether we take “left turn” or “right turn”. The computational complexity of the algorithm is \mathcal{O}(n \log n), where n is the number of points.

1. An Algorithm for Convex Polytopes, D. Chand, S. Kapur, Published in JACM 1970
2. An efficient algorithm for determining the convex hull of a finite planar set, Graham, R. L., Information Processing Letters, 1972.

For algorithms and their complexity, you are encouraged to look at section 2 (Analysis of Algorithms) and section 3 (Asymptotic Order of Growth) of this course by Sushovan Majhi.

2.2.1 Gift Wrapping (*)

viewof width = Inputs.range([100, 800], {step:50, value: 500, label:'Width:'})
viewof height = Inputs.range([100, 800], {step:50, value: 400, label:'Height:'})
viewof num_points = Inputs.range([3, 100], {step:1, value:30, label:'Number of points:'})
viewof delay_gift = Inputs.range([25, 2000], {value:400, step:25, label:'Animation Delay'});
viewof restart_gift = Inputs.button('Restart Animation');
{
  restart_gift;
  var hull = [];
  let current_idx = leftmost_idx, candidate_idx;
  draw_gift(current_idx, -1, -1, hull);
    do {
      hull.push(points[current_idx]);
      candidate_idx = (current_idx + 1) % num_points;
      
      for (let i = 0; i < num_points; i++) { 
        if (i === candidate_idx || i === current_idx) continue;
        yield Promises.delay(delay_gift, draw_gift(current_idx, candidate_idx, i, hull));
        if (clockwiseness(points[current_idx], points[i], points[candidate_idx]) === -1){
          candidate_idx = i;
        }
      }
      current_idx = candidate_idx;
      if(current_idx == leftmost_idx) hull.push(points[current_idx]);
    } while (current_idx !== leftmost_idx); 
  yield draw_gift(-1, -1, -1, hull);
}

function draw_gift(current_idx, candidate_idx, test_idx, hull){
    var ctx = DOM.context2d(width, height);
  
    for(var i = 0; i < num_points; i++){
      ctx.fillStyle = 'Black'
      if(i == current_idx) ctx.fillStyle = 'Red';
      if(i == candidate_idx){
        draw_segment(ctx, points[current_idx][0], points[current_idx][1],
                     points[i][0], points[i][1], 'Green');
        ctx.fillStyle = 'Green';
      }
      if(i == test_idx){
        draw_segment(ctx, points[current_idx][0], points[current_idx][1],
                     points[i][0], points[i][1], 'Blue');
        ctx.fillStyle = 'Blue'
      }
      ctx.beginPath();
      ctx.arc(points[i][0], points[i][1], 5, 0, Math.PI * 2, false);
      ctx.fill();
    }
    
    for(var j = 0; j < hull.length-1; j++){
        draw_segment(ctx, hull[j][0], hull[j][1], hull[j+1][0], hull[j+1][1], 'Red');
    }
    return ctx.canvas;
    
  }
points = {
  var positions = [];
  for(var i = 0; i < num_points; i++){
    var x_pos = Math.random() * (width-20) + 10;
    var y_pos = Math.random() * (height-20) + 10;
    positions.push([x_pos, y_pos]);
    // positions.push([x_pos, y_pos]);
  }
  return positions;
}
sorted_points = [...points].sort((a, b) => get_slope(a, lowest_point) - get_slope(b, lowest_point))
leftmost_idx = {
  var x_vals = points.map(function(v) {return v[0];});
  return x_vals.indexOf(Math.min.apply(null, x_vals));
}
lowest_idx = {
  var y_vals = points.map(function(v) {return v[1];});
  return y_vals.indexOf(Math.max.apply(null, y_vals));
}
lowest_point = points[lowest_idx]
function draw_segment(ctx, x1, y1, x2, y2, color){
        ctx.linewidth = 4;
        ctx.strokeStyle = color;
        ctx.beginPath();
        ctx.moveTo(x1, y1);
        ctx.lineTo(x2, y2);
        ctx.stroke();
    }
function clockwiseness(p0, p1, p2){
  // returns 1 if clockwise, 0 if colinear, -1 if counterclockwise
  var signed_area = (p1[0] - p0[0])*(p2[1]-p0[1]) - (p2[0] - p0[0])*(p1[1]-p0[1])
  if(signed_area == 0) return 0;
  return (signed_area) > 0? 1 : -1
}
function get_slope(p1, p2){
  return p1[1] === p2[1] ? -Infinity : (p2[0] - p1[0]) / (p2[1] - p1[1])
}

2.2.2 Graham Scan (*)

viewof delay_graham = Inputs.range([25, 2000], {value:400, step:25, label:'Animation Delay'})
viewof restart_graham = Inputs.button('Restart Animation')
{
  restart_graham;
  let hull = [sorted_points[0], sorted_points[1], sorted_points[2]];
  for(var i=3; i<=num_points; i++){
    while(clockwiseness(hull[hull.length-2], hull[hull.length-1], sorted_points[i%num_points]) == 1){
      hull.pop();
      yield Promises.delay(delay_graham, draw_graham(hull, i+1));
    }
    hull.push(sorted_points[i%num_points]);
    yield Promises.delay(delay_graham, draw_graham(hull, i+1));
  }
  return hull;
}

function draw_graham(hull, candidate_idx){
  const ctx = DOM.context2d(width, height);

  for(var i = 0; i < num_points; i++){
    ctx.beginPath();
    ctx.fillStyle = 'Black'
    ctx.arc(sorted_points[i][0], sorted_points[i][1], 5, 0, 2 * Math.PI);
    ctx.fill();
  };
  ctx.lineWidth = 4;

  for(var j = 0; j < hull.length - 1; j++){
    draw_segment(ctx, hull[j][0], hull[j][1], hull[j+1][0], hull[j+1][1], 'Red');
  }
  ctx.lineWidth = 2;
  if(candidate_idx != num_points+1) draw_segment(ctx, hull[hull.length-2][0], hull[hull.length-2][1], sorted_points[candidate_idx%num_points][0], sorted_points[candidate_idx%num_points][1], 'blue')
  
  return ctx.canvas;
}

(*) observablehq by Sylvio Jorge Pastene Helt

2.3 Triangulating the convex hull of a finite planar set of points

Once we have found the convex hull of a finite point set, we explore ways of tringulating it.

A triangulation of a finite subset S of the plane is a subdivision of the convex hull of the set S by a maximal set of non-crossing edges. (The edges can intersect only at the points of S, and the maximality means that any new edge will cross some existing edges in the interior)

Triangulating a convex hull: Two algorithms

1. Triangle Spliting Given a finite set of points in general position. We first find the convex hull of the points using one of the algorithms above, and then triangulate the hull (recall that it is a convex polygon) - while ignoring the internal points. Now considering each internal point, add edegs from the point to the three vertices of the triangle containing it (recell that each such internal point is in the interior of some triangle containing it). Do this for each internal point sequentially till all the points are exhausted.

2. Line Sweep Given a finite set of points in general position (placed so that no two points have the same horizontal coordinate). Start with a vertical line to the left of all the points and start moving it rightwards.vEvery time the linea reaches a new point of the the point set, we addd all possible edges towards the left without creating any edge intersections.

Exercise 2.3  

1. Calculate the complexity of the two triangulation algorithms given above.

Measuring the fatness of triangulations, and flipping edges

1. Let S be a finite set of points in \mathbb{R}^2. In a given triangulation of S (with n triangles), we organize the angles of the triangles in decreasing order of magnitude. Given two triangulations T and T^\prime of S we say that T is fatter than T^\prime iff the angle sequence of T is lexicographically larger than T^\prime.
2. Consider edge e of triangulation T of a point set, and consider the quadrilateral formed by the two triangles that have e as a common edge. If Q is convex, let T^\prime is the triangulation after flipping the edge e. We say that e is a good edge if T is fatter than T^\prime otherwise call it a bad edge.
   

Every triangulation of a finite set of points has the same number of triangles.

Delaunay Triangulation

The Delaunay triangulation of a set of points is a triangulation that has only good edges.

Exercise 2.4  

1. Let T and T^\prime be triangulations of a finite set of points. Show that T can be transformed into T^\prime by a finite sequence of edge flips.
2. Let S be a finite set of points in \mathbb{R}^2 such that no 4 points lie on the same circle. Show that a triangulation T is a Delaunay triangulation iff no point of S is in the interior of any circumcenter of a triangle of the triangulation T.

Boris Delaunay was a mathematician with contributions to geometry and mathematical crystallography.

2.3.1 Delaunay Triangulation and its Dual Voronoi Diagram

viewof config = {
  const form = html`<form style="font: 12px var(--sans-serif);">
  <label style="display: block;">
    <input name=s type=range min=1 max=140 value=140 step=0.1 style="width:180px;">
    cell spacing = <output name=os></output>
  </label>
  <label style="display: block;">
    <input name=l type=range min=1 max=200 value=100 step=1 style="width:180px;">
    displacement = <output name=ol></output>
  </label>
  <label style="display: block;">
    <input name=f type=range min=0.001 max=0.005 value=0.002 step=0.0001 style="width:180px;">
    frequency = <output name=of></output>
  </label>
</form>`;
  form.oninput = () => {
    form.value = {
      s: form.s.valueAsNumber,
      l: form.l.valueAsNumber,
      f: form.f.valueAsNumber
    };
    form.os.value = `${form.value.s.toFixed(1)}px`;
    form.ol.value = `${form.value.l.toFixed(0)}px`;
    form.of.value = `${form.value.f.toFixed(3)}`;
  };
  form.oninput();
  return form;
}
viewof run = checkbox({options: ['Simulate'], value: 'Simulate'})
viewof z = slider({min: 100, max: 200, step: 1, value: 150,
                   title: "Tweak the graph",
                   description: 'Drag the slider to animate delaunay/voronoi dual below'})
canvas = {
  const context = DOM.context2d(width, height);
  context.globalCompositeOperation = 'destination-over';
  context.canvas.id = 'canvasDV';
  const coordinates = new Float64Array(points1.length * 2);
  const delaunay = new d3.Delaunay(coordinates);
  const voronoi = delaunay.voronoi([0, 0, width, height]);
  while (true) {
    context.clearRect(0, 0, width, height);
    for (let i = 0, j = -1; i < points1.length; ++i) {
      const [x, y] = points1[i];
      const t = noise(x * config.f, y * config.f, z*0.005) * Math.PI;
      coordinates[++j] = x + config.l * Math.cos(t);
      coordinates[++j] = y + config.l * Math.sin(t);
    }
    context.beginPath();
    voronoi.update().render(context);
    context.strokeStyle = 'red';
    context.stroke();
    delaunay.update().render(context);
    context.strokeStyle = 'purple';
    context.stroke();
    
    yield context.canvas;
    // if (!run) break;
  }
}
points1 = d3.hexbin()
    .radius(config.s)
    .extent([[-config.l, -config.l], [width + config.l, height + config.l]])
  .centers()
noise = {
  const simplex = new (await require("simplex-noise@2"));
  return simplex.noise3D.bind(simplex);
}
d3 = require("d3-hexbin@0.2", "d3-delaunay@5")
//viewof restart = button({value: 'restart'})
import {slider, button, checkbox} from '@jashkenas/inputs'
import { fullscreen } from '@bryangingechen/fullscreen-canvas'

Delaunay Triangulation: Algorithms

Edge Flip: Given a finite set of points in general position, such that no four points are concylic - create a traingulation using one of the algorithms above. Go through the internal (non-hull) edges of the triangulation and if the edge is not good, flip it to make a good edge. Continue doing this till you go through all the edges.

Exercise 2.5  

1. Calculate the complexity of the algorithm given above.
2. Is it possible - in the course of the Delaunay Triangualtion algorithm above - for an edge to be good and later become bad?

2.4 Voronoi Diagrams

Voronoi regions and Voronoi Diagrams

1. Let S be a finite set in \mathbb{R}^2 ( we will call the points “sites” or “seeds”). The Voronoi region \textbf{Vor} (s_i) is defined as the set of points of \mathbb{R}^2 which are at most as close of s_i than all the other points of S. That is,

\textbf{Vor} (s_i) =\{ z \in \mathbb{R}^2: \|x-s_i\| \le \|x-s_j\|, \text{ where } j=1,2\cdots, n \}

2. The Voronoi diagram of S,\textbf{Vor} (S), is the union of the boundaries of all the Voronoi regions of seeds from S.

Georgy Voronyi was a mathematician with contributions to several areas of mathematics, including tessellations.

Exercise 2.6  

1. Show that a Voronoi region is convex.
2. Show that the edges of a Voronoi diagram for a configuration of points (on the plane) can be line segments, rays or infinite lines. Give examples of each for actual configuarations of points.
3. Can a Voronoi diagram be a disconnected graph?

Voronoi vertices and Voronoi edges

Let S be a finite set of sites in \mathbb{R}^2, with Voronoi diagram \textbf{Vor} (S).

1. A point v \in \mathbb{R}^2 is a Voronoi vertex of \textbf{Vor} (S) iff there exists a circle centered at v with at least three sites on its circumference and with no sites in its interior.

2. A connected subset E of a bisector of two sites s_i and s_j is a Voronoi edge of \textbf{Vor} (S) iff for eachpoint x \in E, the circle centered at x through s_i and s_j contains no other sites in its circumference or interior.