{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "delaunay.ipynb\n", "\n", "Discussion: This Jupyter notebook investigates triangulations.\n", "\n", "Licensing: This code is distributed under the GNU LGPL license.\n", " \n", "Modified: 29 October 2016\n", "\n", "Author: John Burkardt, Lukas Bystricky" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using matplotlib backend: agg\n" ] } ], "source": [ "# Import necessary libraries and set plot option\n", "%matplotlib\n", "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import math\n", "import scipy.spatial as spatial" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Triangulations #\n", "\n", "In this module we will look at some problems related to the idea of\n", "organizing a set of points into triangles. Such an arrangement is\n", "called a triangulation, and normally, it is understood that it has\n", "the properties that no two triangle edges cross, and that no more\n", "triangles can be can be added.\n", "\n", "We will look at\n", "* how a triangulation can be represented;\n", "* how to write a function to read triangulation files;\n", "* how to plot a triangulation\n", "* how a Delaunay triangulation is the \"best\" of all possible triangulations;\n", "* how the polygonal shape created by the triangulation can be analyzed;\n", "* how the Delaunay triangulation is related to the Voronoi diagram." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Representing a triangulation #\n", "\n", "Ordinarily, before we have a triangulation, we first have a set of N nodes,\n", "whose coordinates are stored in an Nx2 array. The natural way to describe\n", "a triangulation, then, is as an Mx3 array, with each row corresponding to\n", "a particular triangle, described by listing in counterclockwise order the\n", "indices of the nodes that form it.\n", "\n", "Download the files:\n", " http://people.sc.fsu.edu/~jburkardt/classes/urop_2016/ten_tri1.png\n", "This is a picture of a triangulation. There are N=10 nodes and M=10 triangles.\n", "Choose a numbering for the points and for the triangles, and then create\n", "the Nx2 node array and the Mx3 triangle array that correspond to this picture\n", "and print them. \n", "\n", "In the node coordinate array \"xy\", the order of the data is up to you; the order\n", "you choose here will affect the values you put into the triangle array \"t\"." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "xy = np.array ( [ \\\n", " [ 0.0, 4.0 ], \\\n", " [ 1.0, 13.0 ], \\\n", " [ 4.0, 7.0 ], \\\n", " [ 5.0, 2.0 ], \\\n", " [ 6.0, 15.0 ], \\\n", " [ 9.0, 10.0 ], \\\n", " [ 12.0, 3.0 ], \\\n", " [ 13.0, 14.0 ], \\\n", " [ 15.0, 5.0 ], \\\n", " [ 17.0, 11.0 ] ] )\n", "\n", "t = np.array ( [ \\\n", " [ 0, 3, 6 ], \\\n", " [ 0, 6, 2 ], \\\n", " [ 2, 6, 8 ], \\\n", " [ 0, 2, 1 ], \\\n", " [ 2, 9, 5 ], \\\n", " [ 1, 2, 5 ], \\\n", " [ 2, 8, 9 ], \\\n", " [ 1, 5, 9 ], \\\n", " [ 1, 9, 4 ], \\\n", " [ 4, 9, 7 ] ] )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Plotting a triangulation #\n", "\n", "After you get the arrays set up, see if you can make a plot of the triangulation\n", "using the \"triplot()\" function:\n", "\n", " plt.figure ( )\n", " plt.gca().set_aspect('equal')\n", " plt.triplot ( xy[:,0], xy[:,1], t )" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[,\n", " ]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure ( )\n", "plt.gca().set_aspect('equal')\n", "plt.triplot ( xy[:,0], xy[:,1], t )" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Reading a Triangulation from a file #\n", "\n", "For our next experiment, we will try to read node and triangle\n", "information from a file.\n", "\n", "You will need to copy the following files:\n", " http://people.sc.fsu.edu/~jburkardt/classes/urop_2016/ten_nodes.txt\n", " http://people.sc.fsu.edu/~jburkardt/classes/urop_2016/ten_tri2.txt\n", "\n", "The first file contains a 10x2 table of node coordinates.\n", "The second contains a 10x3 table of triangles (triples of node indices)\n", "\n", "Try to write programs\n", "\n", "def read_nodes ( filename ):\n", " ***\n", " return xy\n", "\n", "and\n", "\n", "def read_triangles ( filename ):\n", " ***\n", " return t\n", "\n", "which get the data from the files." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.0, 4.0], [1.0, 13.0], [4.0, 7.0], [5.0, 2.0], [6.0, 15.0], [9.0, 10.0], [12.0, 3.0], [13.0, 14.0], [15.0, 5.0], [17.0, 11.0]]\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Read the node file #\n", "#\n", "def read_nodes ( filename ):\n", " \n", " input = open ( filename, 'r' )\n", "\n", " xy_list= []\n", " for line in input:\n", " line = line.strip ( )\n", " xy_list.append ( [ float ( x ) for x in line.split ( ) ] )\n", "\n", " input.close ( )\n", "\n", " xy = np.asarray ( xy_list )\n", " \n", " return xy\n", "\n", "print ( xy )\n", "#\n", "# Test the function.\n", "#\n", "xy = read_nodes ( 'ten_nodes.txt' )\n", "plt.plot ( xy[:,0], xy[:,1], 'b*' )" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0 3 6]\n", " [0 6 2]\n", " [2 6 8]\n", " [0 2 1]\n", " [2 9 5]\n", " [1 2 5]\n", " [2 8 9]\n", " [1 5 9]\n", " [1 9 4]\n", " [4 9 7]]\n" ] }, { "data": { "text/plain": [ "[,\n", " ]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Read the triangle file #\n", "#\n", "def read_triangles ( filename ):\n", " \n", " input = open ( filename, 'r' )\n", "\n", " t_list= []\n", " for line in input:\n", " line = line.strip ( )\n", " t_list.append ( [ int ( x ) for x in line.split ( ) ] )\n", "\n", " input.close ( )\n", "\n", " t = np.asarray ( t_list )\n", " \n", " return t\n", "\n", "print ( t )\n", "#\n", "# Test the function.\n", "#\n", "xy = read_nodes ( 'ten_nodes.txt' )\n", "t = read_triangles ( 'ten_tri1.txt' )\n", "\n", "plt.figure ( )\n", "plt.gca().set_aspect('equal')\n", "plt.triplot ( xy[:,0], xy[:,1], t )" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# The minimum angle of a triangulation #\n", "\n", "For a given set of points, there are many possible triangulations.\n", "Is there any reason to prefer one triangulation over another?\n", "One feature that is usually undesirable is the existence of\n", "small angles in the triangles. Suppose that for any triangulation\n", "T we let Theta(T) be the smallest angle in any triangle in that\n", "triangulation. Then we could say that triangulation T2 is better\n", "than T1 if Theta(T1) < Theta(T2), that is, the smallest angle in\n", "T2 is bigger than that in T1.\n", "\n", "Write a function which takes the values of x, y, and t, and\n", "computes the value of Theta(T1), of the following form:\n", "\n", "def theta ( xy, t ):\n", " ***\n", " return min_angle\n", "\n", "Test your function on the triangulation T1 defined by the \n", "nodes \"ten_nodes.txt\" and \"ten_tri1.txt\"." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Minimum angle for triangulation T1 is 9.97771 degrees\n" ] } ], "source": [ "# Compute minimum angle of a triangulation\n", "#\n", "# Make output in degrees because no one thinks in radians!\n", "#\n", "\n", "def theta ( xy, t ):\n", " import numpy as np\n", " from numpy.linalg import norm\n", " n = xy.shape[0]\n", " m = t.shape[0]\n", " min_angle = np.pi\n", " for j in range ( 0, m ):\n", " a = t[j,0]\n", " b = t[j,1]\n", " c = t[j,2]\n", " bma = norm ( xy[b,:] - xy[a,:] )\n", " cmb = norm ( xy[c,:] - xy[b,:] )\n", " amc = norm ( xy[a,:] - xy[c,:] )\n", " cosa = ( bma**2 + amc**2 - cmb**2 ) / ( 2 * bma * amc )\n", " cosb = ( cmb**2 + bma**2 - amc**2 ) / ( 2 * cmb * bma )\n", " cosc = ( amc**2 + cmb**2 - bma**2 ) / ( 2 * amc * cmb )\n", " anga = np.arccos ( cosa )\n", " angb = np.arccos ( cosb )\n", " angc = np.arccos ( cosc )\n", " min_angle = min ( min_angle, anga )\n", " min_angle = min ( min_angle, angb )\n", " min_angle = min ( min_angle, angc )\n", " min_angle = min_angle * 180.0 / np.pi\n", " return min_angle\n", "#\n", "# Test on our sample triangulation.\n", "#\n", "xy = read_nodes ( 'ten_nodes.txt' )\n", "t = read_triangles ( 'ten_tri1.txt' )\n", "min_angle = theta ( xy, t )\n", "print ( 'Minimum angle for triangulation T1 is %g degrees' % ( min_angle ))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The Delaunay Triangulation #\n", "\n", "The Delaunay triangulation T of a given set of points maximizes \n", "the value of Theta. For this reason, given a set of points, we \n", "(almost) always choose the Delaunay triangulation.\n", "\n", "It turns out that for the set of points in the file \"ten_nodes.txt\",\n", "the triangles in \"ten_tri1.txt\" are a \"bad\" triangulation, while\n", "the triangles in \"ten_tri2.txt\" are the Delaunay triangulation.\n", "\n", "Use your theta function to evaluate the minimum angle in triangulation T2.\n", "\n", "Also, plot triangulation T2, and look at how it differs from T1.\n", "You should agree that the minimum angle in T2 is bigger than the\n", "minimum angle in T1." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Minimum angle for triangulation T2 is 26.9958 degrees\n" ] }, { "data": { "text/plain": [ "[,\n", " ]" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#\n", "# Test on our sample triangulation.\n", "#\n", "xy = read_nodes ( 'ten_nodes.txt' )\n", "t = read_triangles ( 'ten_tri2.txt' )\n", "min_angle = theta ( xy, t )\n", "print ( 'Minimum angle for triangulation T2 is %g degrees' % ( min_angle ))\n", "\n", "plt.figure ( )\n", "plt.gca().set_aspect('equal')\n", "plt.triplot ( xy[:,0], xy[:,1], t )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Computing Properties of a triangulated region #\n", "\n", "In the polygons notebook, we were able to compute properties\n", "of the polygon, such as area and centroid, by working\n", "with the triangles. The same approach works for a region\n", "created by starting with points and determining their\n", "triangulation. One difference is that a shape created\n", "by triangulating points will always be convex, that is,\n", "there will be no notches or indentations in the boundary,\n", "and the triangulated shape will probably have a number of\n", "points and triangles in the interior; a triangulated\n", "polygon has all its points on the boundary, and it\n", "can have a very nonconvex (indented) boundary. But\n", "the area and centroid are computed the same way.\n", "\n", "For the same reason, sampling is done the same way,\n", "and so integrals can also be approximated over the region.\n", "I won't make you repeat these calculations for the triangulated\n", "region but stop for one moment, go back to the polygon notebook\n", "and convince yourself that everything you did there will still\n", "work for this new but very similar case." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Relationship between Voronoi Diagram and Delaunay Triangulation\n", "\n", "Mathematically, the Voronoi diagram and the Delaunay\n", "triangulation are duals of each other. The information\n", "from one object can be used to reconstruct the other one.\n", "\n", "One way to see this is to start with a set of points and on\n", "a single plot show both the Delaunay triangulation and the\n", "Voronoi diagram. The generator points can be regarded as\n", "the \"capitals\" of states defined by the Voronoi polygons.\n", "The Delaunay triangulation essentially builds roads between\n", "neighboring capital cities. Every road crosses the state\n", "line at a 90 degree angle, and is half in one state and half\n", "in the other.\n", "\n", "As you saw in the 2D geometry notebook, it's possible to compute\n", "the Voronoi diagram and then draw it by the commands\n", " vor = scipy.spatial.Voronoi ( xy )\n", " scipy.spatial.voronoi_plot_2d ( vor )\n", "and the Delaunay triangulation can be compute and drawn by\n", " t = scipy.spatial.Delaunay ( xy )\n", " matplotlib.pyplot.triplot ( xy[:,0], xy[:,1], t )\n", "However, I haven't figure out how to show both figures on a single plot.\n", "\n", "For the xy data we have, compute and plot the Voronoi diagram and\n", "Delaunay triangulation.\n", "\n", "Oh my, now it works!" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[,\n", " ]" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "vor = spatial.Voronoi ( xy )\n", "spatial.voronoi_plot_2d ( vor )\n", "t = spatial.Delaunay ( xy )\n", "plt.triplot ( xy[:,0], xy[:,1], t.simplices.copy() )" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.8" } }, "nbformat": 4, "nbformat_minor": 0 }