{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "polygons.ipynb\n", "\n", "Discussion: This Jupyter notebook investigates polygons.\n", "\n", "Licensing: This code is distributed under the GNU LGPL license.\n", " \n", "Modified: 26 October 2016\n", "\n", "Author: John Burkardt, Lukas Bystricky" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Import necessary libraries and set plot option\n", "%matplotlib inline\n", "%config InlineBackend.figure_format = 'png'\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": [ "# Polygons\n", "\n", "In this module we will look at how our understanding of triangles\n", "allows us to answer questions about polygons.\n", "\n", "We will assume that we have a description of a polygon as a\n", "sum of triangles, or that some special software has been used\n", "which starts with the polygonal vertex coordinates and determines\n", "a triangulation.\n", "\n", "Then we look at how we can turn many questions about polygons\n", "into a question about one of the triangles composing the polygon.\n", "\n", "While some of these procedures will be obvious, the task of\n", "uniform random sampling of a polygon will require a little extra work." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEACAYAAACatzzfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE3xJREFUeJzt3W+sXHWdx/HPt9YSkAWirGBo+LuuWlkVVAqUwLhEFyUB\nn5RI3Bh4QILZZfljDBseQPcJcXkCuxpjRLZZCG6pBgSjWTDKaAoFuxQQaYvERXbLQiOBQBpSWuS7\nD84M93Q6986ZmXPO7/c7v/cruemUe+7cn9d7P/n0MzOtubsAAPFaFvoAAIClEdQAEDmCGgAiR1AD\nQOQIagCIHEENAJFbXuUiM/uDpNckvS1pn7uf1uShAAALKgW1ioDuufurTR4GAHCgqtOHTXEtAKBG\nVcPXJf3MzLaY2WVNHggAsL+q08cad3/RzP5cRWBvd/dNTR4MAFCoFNTu/uLg1z+a2T2STpO0X1Cb\nGX9pCABMyd1t0jUTpw8zO8TMDh3cfo+kz0n67SKfkDd33XDDDcHPEOTttdfkZ54pV7GV3bBqVfgz\nRfKW7fcEX4sl36qq0qiPknTPoDEvl3Snuz9Q+TMgD6+/Ln3+89LDDy/8t3POCXceoEMmBrW7Pyfp\nEy2cBakaF9Jr10rvf3+4MwEdUvXBREyh1+uFPkJ7yiG9YoW0d2/x36+/Xr2XXw57tohk9T0xAV+L\n6dk0O8mSd2Tmdd0XElEO6WOPlU4+WfrpT4s2vXFj6NMB0TMzeYUHEwlqzGY0pO+8Uzr33KJRP/VU\nEdoAllQ1qHm1IaY3GtIPPijddVcR0mvXEtJAzWjUmM64kD7oIOnEE2nTwJSqNmoeTER140L6xBOl\nK66gTQMNolGjmsVC+oUXaNPAjNioUZ/FQlqSvvEN2jTQMBo1lrZUSNOmgbnQqDG/pUJaok0DLaFR\nY7xJIU2bBuZGo8bsJoW0RJsGWkSjxv6qhDRtGqgFjRrTqxLSEm0aaBmNGoWqIU2bBmpDo0Z1VUNa\nok0DAdCoczdNSNOmgVrRqDHZNCEt0aaBQGjUuZo2pGnTQO1o1FjctCEt0aaBgGjUuZklpGnTQCNo\n1DjQLCEt0aaBwGjUuZg1pGnTQGNo1Fgwa0hLtGkgAjTqrpsnpGnTQKNo1JgvpCXaNBAJGnVXzRvS\ntGmgcTTqnM0b0hJtGogIjbpr6ghp2jTQChp1juoIaYk2DUSGRt0VdYU0bRpoDY06J3WFtESbBiJE\no05dnSFNmwZaRaPOQZ0hLdGmgUjRqFNVd0jTpoHW0ai7rO6QlmjTQMQqN2ozWybpvyTtdPcLxryf\nRt2GJkKaNg0E0USjvlLSttmPhLk1EdISbRqIXKWgNrOVkr4g6XvNHgeLaiqkX3hB+u53i9vXXz//\n/QGoXdVGfbOkr0ti2wihqZCWaNNAApZPusDMzpe0y92fMLOepEX3lHXr1r1zu9frqdfrzX/C3DUZ\n0rRpoFX9fl/9fn/qj5v4YKKZ3SjpbyW9JelgSX8m6W53/8rIdTyYWLcmQ1qSrrhC+ta3ija9cWN9\n9wugkqoPJk71PGozO0fS13jWRwuaDmme6QEEx/OoU9Z0SEts00BCeGVibNoIado0EAUadYraCGmJ\nNg0khkYdi7ZCmjYNRINGnZK2QlqiTQMJolGH1mZI06aBqNCoU9BmSEu0aSBRNOpQ2g5p2jQQHRp1\nzNoOaYk2DSSMRt22ECFNmwaiRKOOUYiQlmjTQOJo1G0JFdK0aSBaNOqYhAppiTYNdACNumkhQ5o2\nDUSNRh2DkCEt0aaBjqBRNyV0SNOmgejRqEMKHdISbRroEBp13WIIado0kAQadQgxhLREmwY6hkZd\nl1hCmjYNJING3aZYQlqiTQMdRKOeV0whTZsGkkKjbkNMIS3RpoGOolHPKraQpk0DyaFRNym2kJZo\n00CH0ainFWNI06aBJNGomxBjSEu0aaDjaNRVxRrStGkgWTTqOsUa0hJtGsgAjXqSmEOaNg0kjUZd\nh5hDWqJNA5mgUS8m9pCmTQPJo1HPI/aQlmjTQEZo1KNSCGnaNNAJNOpZpBDSEm0ayAyNeiiVkKZN\nA51Bo55GKiEt0aaBDE1s1GZ2kKRfSVoxeLvX3a8bc12ajTqlkKZNA51StVEvn3SBu79pZp9x9zfM\n7F2SHjKzNe7+UC0nDSmlkJZo00CmptqozewQSX1Jl7j7tpH3pdWoUwtp2jTQObVu1Ga2zMwel/SS\npP5oSCcntZCWaNNA17z1VuVLp23Uh0l6QNK17v7Lkfel0ahTDGnaNJCuffuk3/9eevrp4m3btuLX\nZ56R7dtXz0Zd5u6vm9lPJH1K0i9H379u3bp3bvd6PfV6vWnuvnkphrREmwZSsEQga98+ScVu3J/h\nrqs86+NISfvc/TUzO1jS/ZL+yd1/PnJd3I061ZCmTQNxqRDIBzjuOOmjHy3e3v1u6cYbJUkm1dao\nPyDp383MVGzad4yGdPRSDWmJNg2EMm8gr1pV/PqRj0iHHlq8f88e6fTTi9sXXSRt3FjpKN1/ZWLK\nIU2bBprXRCAv5qtflb7zHemkk6StW2WHH17/Rp2clENaok0DdWozkMfZsKEI6RUrpB/8QDrssMof\n2t1GnXpI06aB2YQO5HGefVY69VRp927p298umrVqfGViklIPaYk2DUwSYyCPs2dP8XO8e3exS19+\n+dR30b1G3YWQpk0DC1IJ5MWM7NLlySPPRt2FkJZo08hT6oE8zhy7dFl3GnVXQpo2ja7rYiCPs8gu\nXZZXo+5KSEu0aXRHLoE8Tg27dFn6jbpLIU2bRopyDuTFLLFLl+XRqLsU0hJtGnEjkKupaZcuS7dR\ndy2kadOIBYE8uwq7dFm3G3XXQlqiTaN9BHK9at6ly9Jr1F0Mado0mkQgt6PiLl3WzUbdxZCWaNOo\nB4EcTgO7dFk6jbqrIU2bxrQI5LhMuUuXdatRdzWkJdo0Fkcgx6/BXbos/kbd5ZCmTUMikFM2wy5d\n1o1G3eWQlmjTuSGQu6XhXbos3kbd9ZCmTXcXgdx9c+zSZWk36q6HtESb7gICOU8t7dJl8TXqHEKa\nNp0WAhllc+7SZWk26hxCWqJNx4pAxiQt7tJl8TTqXEKaNh0egYxZ1LRLl6XVqHMJaYk23SYCGXUJ\nsEuXhW/UOYU0bboZBDKaVuMuXZZGo84ppCXa9LwIZIQQaJcuC9eocwtp2nR1BDJi0cAuXRZ3o84t\npCXa9DgEMmIWeJcua79R5xjSubdpAhkpamiXLouzUecY0lI+bZpARldEsEuXtdeocw3pLrZpAhld\n1vAuXRZXo841pKW02zSBjNxEtEuXNd+ocw7pVNo0gQwUWtily+Jo1DmHtBRfmyaQgcVFtkuXNdeo\ncw/pkG2aQAam0+IuXRa2Uece0lI7bZpABuYX6S5dVn+jJqTrb9MEMtCclnfpstoatZmtlHS7pKMk\nvS3pVnf/17EXE9KFWds0gQy0K+JdumxiozazoyUd7e5PmNmhkh6TdKG77xi5zv3MMwnpKm2aQAbC\nC7RLl9XWqN39JUkvDW7vNrPtko6RtOOAi3MPaWn/Nv2hD0k7dhDIQGwS2KXLptqozex4SX1JJ7v7\n7pH3uece0uU2fdRR0iuvEMhAjALu0mW1P+tjMHv8UNKVoyH9jvvvzzekpeJ//969xe1du4pfCWQg\nLons0mWVgtrMlqsI6Tvc/d7Frlt3ySXSeedJknq9nnq9Xg1HTMjatUWLPvJIAhmI0bPPSpddVty+\n5RbplFNa/fT9fl/9fn/qj6s0fZjZ7ZJedvdrlrimuKd77pG++MWpDwIAjdqzRzr9dOnJJ4tdesMG\nySauDo2qOn1UedbHGkm/kvSUJB+8Xefu/zlyXXFPRxxRbD4nnDDr2QGgfpHs0mW1BfUUn9D9gguk\n++6TPv1padOmYgMCgNA2bJAuvrjIpEceaX3yWEzVoF5W62ddv7548GzLFunaa2u9awCYSeBdug71\nv4T80Uels86S3nqLvRpAWBHu0mVhGrUkrV4t3XRTcfvSS6Xnnqv9UwBAJVdfXYT0SSdJt94aVUhP\no5m/5tS9aNLs1QBCiXSXLgvXqIvPzl4NIJwO7NJlzf5TXOzVANoW+S5dFrZRD7FXA2hbR3bpsub/\ncVv2agBtSWCXLoujURcnYa8G0LyO7dJlzTfqIfZqAE1JaJcui6dRD7FXA2hKB3fpsvYatcReDaB+\nie3SZfE1aom9GkC9OrxLl7XbqIfYqwHMK9FduizORj3EXg1gXh3fpcvCNGqJvRrA7BLepcvibtQS\nezWA2WSyS5eFa9RD7NUAqurALl0Wf6MeYq8GUFVGu3RZ+EYtsVcDmKwju3RZOo1aYq8GsLQMd+my\nOBr1EHs1gFEd26XL0mrUQ+zVAEZlukuXxdWoJfZqAAs6uEuXpdmoJfZqAIXMd+my+Br1EHs1kK8O\n79Jl6TbqIfZqIF/s0vuJt1FL7NVAjjq+S5el36gl9mogN+zSY8XdqIfYq4Huy2SXLutGox5irwa6\nj116UWk0aom9GuiyjHbpsm41aom9GugqdumJ0mnUQ+zVQHdkuEuXda9RD7FXA93BLl1Jeo1aYq8G\nuiDTXbqstkZtZreZ2S4z+009R6sBezWQNnbpqVSZPtZL+pumDzK1975Xuusuafny4v/oH/0o9IkA\nVLFnj7R2rbR7d7FLX3556BNFb2JQu/smSa+2cJbpsVcD6WGXnlqljdrMjpP0Y3f/2BLXtLdRl7FX\nA3F7803p8celzZuLn8+77856ly6rulEvr/OTrlu37p3bvV5PvV6vzrsfb7hXn3rqwl59883Nf14A\n4+3cWYTy5s1FGD/2mLR378L7zaRvfjPLkO73++r3+1N/XPqNeojnVwPtK7fl4dvOnQdet2qVdMYZ\nxdvZZ0sf/GD7Z41Q1UZdNaiPVxHUf7XENWGDWiqa9DXXSEccIW3dKp1wQtjzAF0zqS1L0uGHF48f\nDYN59eriZxIHqC2ozez7knqS3idpl6Qb3H39mOvCBzV7NVCfWdryGWdIH/6wtCy919KFUGujrvgJ\nwwe1JL3ySrFXP/+8dNVV7NVAVbTl1uUb1BJ7NTAJbTkKeQe1xF4NlNGWo0RQs1cjV7TlZBDUEns1\n8kBbThZBPcRejS6hLXcKQV3GXo1U0ZY7jaAuY69GCmjL2SGoR7FXIza05ewR1OOwVyMU2jLGIKgX\nw16NNpTb8ubNxfcabRkjCOrFsFejbrRlzIigXgp7NeZBW0ZNCOpJ2KtRBW0ZDSKoq2CvxijaMlpE\nUFfBXp032jICI6irYq/OB20ZkSGop8Fe3T20ZSSAoJ4We3XaaMtIEEE9LfbqdNCW0REE9SzYq+NE\nW0ZHEdSzYq8Oi7aMjBDU82Cvbg9tGRkjqOfBXt0M2jKwH4J6XuzV86MtA0siqOvAXl3dm28WQfzI\nI7RloCKCui7s1ePRloG5EdR1Ya9eaMvDfzKKtgzUgqCuU257NW0ZaAVBXbeu7tW0ZSAYgroJXdir\nactANAjqJqS2V9OWgagR1E2Jea+mLQNJIaibFMNeTVsGkkdQN63tvZq2DHQOQd20Jvdq2jKQhVqD\n2szOk3SLpGWSbnP3fx5zTV5BLdW3V9OWgSzVFtRmtkzS7ySdK+n/JG2R9CV33zFyXX5BLY3dq/v9\nvnq93vjrM2rLS34dMsPXYgFfiwVVg3p5hfs6TdKz7v784I43SLpQ0o4lPyoXq1dLN91U7NWXXip9\n/OP7fyNm3Jb5gVzA12IBX4vpVQnqYyT9b+n3O1WEN4auukrq94u9+vzzpWOOkS66qNNtGUB7qgQ1\nJjGT1q8vnvmxfXvxNtTRtgygPVU26tMlrXP38wa//0dJPvqAopllOFADwHzqejDxXZKeUfFg4ouS\nfi3pYnffvuQHAgBqMXH6cPc/mdnfS3pAC0/PI6QBoCW1veAFANCMuZ9qYGbnmdkOM/udmV1bx6FS\nZGa3mdkuM/tN6LOEZmYrzewXZva0mT1lZv8Q+kyhmNlBZvaomT0++HrcGPpMoZnZMjPbamb3hT5L\nSGb2BzN7cvC98eslr52nUVd9MUwOzOwsSbsl3e7uHwt9npDM7GhJR7v7E2Z2qKTHJF2Y4/eFJJnZ\nIe7+xuDxnockfc3dHwp9rlDM7GpJn5R0mLtfEPo8oZjZf0v6pLu/OunaeRv1Oy+Gcfd9koYvhsmO\nu2+SNPELngN3f8ndnxjc3i1pu4rn42fJ3d8Y3DxIxc9ctt8nZrZS0hckfS/0WSJgqpjB8wb1uBfD\nZPsDiQOZ2fGSPiHp0bAnCWfwR/3HJb0kqe/u20KfKaCbJX1dEg+OFV+Dn5nZFjO7bKkLeTkcGjOY\nPX4o6cpBs86Su7/t7qdIWinpbDM7J/SZQjCz8yXtGvxpywZvOVvj7qeq+BPG3w3m07HmDeoXJB1b\n+v3KwX9D5sxsuYqQvsPd7w19nhi4++uSfiLpU6HPEsgaSRcMttn/kPQZM7s98JmCcfcXB7/+UdI9\nWuKv5pg3qLdI+gszO87MVkj6kqScH8mlJSz4N0nb3P1fQh8kJDM70swOH9w+WNJnJT0R9lRhuPt1\n7n6su5+oIit+4e5fCX2uEMzskMGfOGVm75H0OUm/Xez6uYLa3f8kafhimKclbcj1xTBm9n1JD0v6\nSzP7HzO7NPSZQjGzNZK+LOmvB0892jr4O81z9AFJDw426kck3efuPw98JoR3lKRNpe+LH7v7A4td\nzAteACByPJgIAJEjqAEgcgQ1AESOoAaAyBHUABA5ghoAIkdQA0DkCGoAiNz/A2HulK78ZXiHAAAA\nAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Representations of a polygon.\n", "#\n", "# There are several ways to represent a polygon.\n", "#\n", "# We will start by assuming that the polygon is simple enough\n", "# that the user can determine a decomposition into triangles.\n", "#\n", "# Here is a diagram of such a polygon, formed by connecting\n", "# A-B-C-D-E-F-A\n", "#\n", "# 5 * E * * * *\n", "# 4 * * * * * *\n", "# 3 * * * * * C\n", "# 2 F * D * * *\n", "# 1 * * * * B *\n", "# 0 * A * * * *\n", "#\n", "# 0 1 2 3 4 5\n", "#\n", "# Notice that, just like we did with triangles, we list the vertices of the\n", "# polygon in counterclockwise order. This is important so that our area\n", "# calculations make sense.\n", "#\n", "p1 = np.array ( [ [1.0,0.0], [4.0,1.0], [5.0,3.0], [1.0,2.0], [2.0,5.0], [0.0,2.0] ])\n", "#\n", "# Plot the outline of the polygon. \n", "# You can almost do this with a single plot command that lists the X and Y components\n", "# of the polygon. However, you need one extra plot statement to take care of the final\n", "# line segment that connects F to A.\n", "#\n", "plt.plot ( p1[0:6,0], p1[0:6,1], linewidth = 2.0, color = 'r' )\n", "plt.plot ( [ p1[5,0], p1[0,0] ], [ p1[5,1], p1[0,1] ], linewidth = 2.0, color = 'r' )\n", "#\n", "# Perhaps we created the vertex coordinates separately.\n", "# Then the polygon can be constructed from them.\n", "#\n", "a = np.array ( [1.0,0.0] ) \n", "b = np.array ( [4.0,1.0] ) \n", "c = np.array ( [5.0,3.0] )\n", "d = np.array ( [1.0,2.0] ) \n", "e = np.array ( [2.0,5.0] )\n", "f = np.array ( [0.0,2.0] )\n", "p2 = np.array ( [ a, b, c, d, e, f ] )\n", "#\n", "# A print statement would show you that P1 and P2 are the same information.\n", "#" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Perimeter and angles #\n", "\n", "There is nothing exciting about computing the perimeter of \n", "a polygon. You simply have more pieces to sum up than you\n", "did for a triangle.\n", "\n", "Similarly, you can use the law of cosines to compute any angle.\n", "To compute the angle alpha associated with vertex A, which\n", "is formed by the vertices F-A-B, you would pretend that you were\n", "looking at a triangle formed by vertices F, A and B, and compute\n", "the lengths of sides F-A, A-B and the \"imaginary side\" F-B. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The Centroid of a polygon #\n", "\n", "We might hope that the calculation of the centroid of a polygon\n", "is as easy as for a triangle, so that the coordinates of the \n", "centroid of our example triangle would be (A+B+C+D+E+F)/6.\n", "This is not correct.\n", "\n", "The convenient thing to keep in mind is that physically, a\n", "triangle behaves the same as though its mass was all located\n", "at the centroid. Now the mass of a triangle is simply its area.\n", "So if triangle T1 has a centroid at (-5,0) and an area of 9 \n", "and triangle T2 has a centroid at (+5,0) and an area of 3, then\n", "we can think of just point masses at those two locations, in\\\n", "which case the balance point would be at 9*(-5,0)+3*(5,0)/(9+3)=-2.5\n", "\n", "In general, if we compute the centroid C and area A of each triangle,\n", "then the centroid of the polygon is\n", " C(polygon) = sum ( A(i) * C(i) ) / sum ( A(i) )\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Areas: 5.000000 8.000000 -2.000000 7.000000\n", "Polygon_area = 18.000000\n", "Polygon_centroid = (2.203704.1.759259)\n" ] } ], "source": [ "# Example computation of centroid\n", "#\n", "# For our 6-vertex example polygon, compute the centroid.\n", "#\n", "# We need to decompose the polygon into triangles. \n", "# One way to do this is to use the following list of vertex triples:\n", "# ABC, ACD, ADE, AEF\n", "# Convince yourself that this divides the polygon into 4 disjoint triangles.\n", "#\n", "# Write a short function that computes triangle areas based on three vertices.\n", "#\n", "# Notice something surprising, (but correct) about one of the areas.\n", "#\n", "# Despite this surprising thing, continue on as though nothing happened,\n", "# and compute the centroid. If we're not sure it's correct because of\n", "# the surprising thing, we will come back and check it later.\n", "#\n", "def triangle_area ( a, b, c ):\n", " area = a[0] * ( b[1] - c[1] )\\\n", " + b[0] * ( c[1] - a[1] )\\\n", " + c[0] * ( a[1] - b[1] )\n", " return area\n", "\n", "abc_area = triangle_area ( a, b, c )\n", "acd_area = triangle_area ( a, c, d )\n", "ade_area = triangle_area ( a, d, e )\n", "aef_area = triangle_area ( a, e, f )\n", "\n", "print ( 'Triangle areas: %f %f %f %f' % ( abc_area, acd_area, ade_area, aef_area ) )\n", "\n", "polygon_area = abc_area + acd_area + ade_area + aef_area\n", "\n", "print ( 'Polygon_area = %f' % ( polygon_area ) )\n", "\n", "polygon_centroid = ( abc_area * ( a + b + c ) / 3.0 \\\n", " + acd_area * ( a + c + d ) / 3.0 \\\n", " + ade_area * ( a + d + e ) / 3.0 \\\n", " + aef_area * ( a + e + f ) / 3.0) / polygon_area\n", "\n", "print ( 'Polygon_centroid = (%f.%f)' % ( polygon_centroid[0], polygon_centroid[1]) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Using only positive area triangles #\n", "\n", "We broke up our example polygon into triangles by simply taking triples\n", "of consecutive vertices. Unfortunately, this led to the disturbing\n", "result that triangle CDE had negative area. As it turns out, we can\n", "still use the sum of the triangle areas to get the polygon area, but\n", "for other calculations, it is actually important that all the triangles\n", "be described with positive area. That means we need a better \n", "triangulation of the polygon.\n", "\n", "Such a triangulation is always possible to find. This is because,\n", "for any \"reasonable\" polygon (we have to forbid polygons whose\n", "lines cross themselves or do other degenerate things) it is always\n", "possible to find at least two \"ears\", that is, three consecutive\n", "vertices ABC such that the line A-C is inside the polygon. We can\n", "imagine slicing ABC from the polygon, which adds a triangle to\n", "our triangulation, and reduces the number of vertices in the polygon\n", "by 1. The smaller polygon will again show at least two \"ears\", and\n", "we can keep slicing one off until the polygon is reduced to a triangle,\n", "which is the last one we need for our collection.\n", "\n", "While the algorithm sounds simple, implementing it can be tricky. Rather\n", "than have you implement it, we will simply have you try it out on our\n", "example polygon. After we have gotten a better triangulation with no\n", "negative areas, we can proceed to the next topic.\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1 5 0]\n", " [3 1 2]\n", " [5 3 4]\n", " [1 3 5]]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEACAYAAACatzzfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFNJJREFUeJzt3X2MZXV9x/H3d11AF1S0ECBFRCVGQXmSoATs3qr1AYwa\nMZGmRtKIf7V20xJTQ0hYIjGKVVobA6HQRkxRI8EKGBXRXnVRHmTBKiuCuz5UC7u4ATc8iMj++se9\nyxyGOzPnzj3nnt855/1KNu4ul9nTyfr2s9+Z2UZKCUlSvtY0/QCSpOUZaknKnKGWpMwZaknKnKGW\npMwZaknK3NoyL4qIXwC/A3YDj6eUTqzzoSRJC0qFmlGgBymlB+p8GEnS05U9fcQUr5UkVahsfBPw\njYi4NSLeX+cDSZKequzp4+SU0r0RcSCjYP8kpbSpzgeTJI2UCnVK6d7xf94fEV8CTgSeEuqI8C8N\nkaQppZRipdesePqIiHURsd/4+/sCbwR+vMQv6LeUOO+88xp/hqa+bd6ceP7zE2vXJjZs6O/7wd8T\nvi/KfCurzKI+CPjSeDGvBf4zpXR96V9BvXH77fCWt8DBB8Ppp8P++zf9RFI3rLioU0o/Tykdm1I6\nLqX0ypTSR+fxYGqXPZH+8IfhvvvgnHOafiKpO/yUuxoMBoOmH2Gu9kT605+GX/5ytKYPP7x/74fl\n+L5Y4PtiejHNnWTZNxSRqnpbao9ipAcDeOlL4bbbRqGWtLyIIFXxwURpKcVIn346XHTRwpqWVB0X\ntVZlcaR37nRNS9NyUas2iyMNrmmpTi5qTWVSpF3T0uq4qFW5SZEG17RUNxe1Slkq0q5pafVc1KrM\nUpEG17Q0Dy5qLWu5SLumpdm4qDWz5SINrmlpXlzUmmilSLumpdm5qLVqK0UaXNPSPLmo9RRlIu2a\nlqrhotbUykQaXNPSvLmoBZSPtGtaqo6LWqWVjTS4pqUmuKh7bppIu6alarmotaJpIg2uaakpLuqe\nmjbSrmmpei5qLWnaSINrWmqSi7pnVhNp17RUDxe1nmY1kQbXtNQ0F3VPrDbSrmmpPi5qPWm1kQbX\ntJQDF3XHzRJp17RULxe1Zoo0uKalXLioO2rWSLumpfq5qHts1kiDa1rKiYu6Y6qItGtamg8XdQ9V\nEWlwTUu5cVF3RFWRdk1L8+Oi7pGqIg2uaSlHLuqWqzLSrmlpvlzUPVBlpME1LeXKRd1SVUfaNS3N\nX+WLOiLWRMTmiLhmtkfTrKqONLimpZytneK1G4AtwHNqehaVUEekd+6Eiy8erWlJ+Sm1qCPiUOBU\n4LJ6H0fLqSPS4JqWcld2UV8EfBB4bo3PomXUFWnXtJS/FUMdEacB21NKd0TEAFjy8L1x48Ynvz8Y\nDBgMBrM/oWqLNLimpXkaDocMh8Op/70VP+sjIj4CvAf4I/As4NnA1Sml9y56nZ/1UYM6I+1nekjN\nKvtZH1N9el5ErAfOTim9bcI/M9QVqzPSAOeeCzt2wKWXVv+2Ja2sbKin+awPzVHdkfY2LbWHX/CS\nobojDa5pKQe1nD5W+AUNdQXmEWlv01Ie/Ls+WmgekQY/00NqGxd1JuYVade0lA8XdYvMK9Lgmpba\nyEXdsHlG2jUt5cVF3QLzjDS4pqW2clE3ZN6Rdk1L+XFRZ2zekQbXtNRmLuo5ayLSrmkpTy7qDDUR\naXBNS23nop6TpiLtmpby5aLOSFORBte01AUu6po1GWnXtJQ3F3UGmow0uKalrnBR16TpSLumpfy5\nqBvUdKTBNS11iYu6YjlE2jUttYOLugE5RBpc01LXuKgrkkukXdNSe7io5yiXSINrWuoiF/WMcoq0\na1pqFxf1HOQUaXBNS13lol6l3CLtmpbax0Vdo9wiDa5pqctc1FPKMdKuaamdXNQ1yDHS4JqWus5F\nXVKukXZNS+3loq5QrpEG17TUBy7qFeQcade01G4u6grkHGlwTUt94aJeQu6Rdk1L7eeinkHukQbX\ntNQnLupF2hBp17TUDS7qVWhDpME1LfWNi3qsLZF2TUvdUXZRry3xhvYBvgPsPf725ZTSObM/Yj7a\nEmlwTUt9VGpRR8S6lNIjEfEM4Ebg7JTSjYte08pF3aZIu6albqn0Rp1SemT83X3G/84DMzxbNtoU\naXBNS3214ukDICLWALcBLwEuSSltqfWp5qBtkd65Ey6+eLSmJfVLqVCnlHYDx0XEc4DrI2J9Sunb\n9T5afdoWaXBNS13wxz/Ctm1w552wZYq5WyrUe6SUdkXEV4ATgKeFeuPGjU9+fzAYMBgMpnnzc9HG\nSLumpXYpBnlPlO+8E+66a8i6dUMOPBC2bi3/9lb8YGJEHAA8nlL6XUQ8C/g6cH5K6ZuLXpf9BxPb\nGGmAc8+FHTvg0kubfhJJRUsF+Z574JBD4Mgj4aijRt+OPBJe9jLYd1/46Efh2mvhe98r98HEMqF+\nJfAZIBh9IPGzKaV/mvC6rEPd1kj7mR5S81Yb5Em++11417vgBz+Aww6rKNRl5RzqtkYaXNPSPFUZ\n5Enuvx+OPx4uuQROO638p+d1PtRtjrRrWqpH3UGeZPfuUZyPPho+9rHRzxlq2h1pcE1Ls2oiyEvZ\nc5ceDmGvvUY/1/tQtz3SrmmpvJyCPEnxLv2CFyz8fK9D3fZIg2tamiT3IE+y+C5d1NtQdyHSrmn1\nXRuDPMmku3RRL0PdhUiDa1r90ZUgL2XSXbqod6HuSqRd0+qirgd5kqXu0kW9CnVXIg2uabVbH4M8\nyXJ36aLehLpLkXZNqy0M8tJWuksX9SLUXYo0uKaVH4M8vZXu0kWdD3XXIu2aVpMMcjXK3KWLOh3q\nrkUaXNOaD4Ncn7J36aLOhrqLkXZNq2oGeb6muUsXdTLUXYw0uKa1egY5D9PcpYs6F+quRto1rTIM\ncr6mvUsXdSrUXY00uKb1VAa5XVZzly7qTKi7HGnXdH8Z5PZb7V26qBOh7nKkwTXdBwa5u1Z7ly5q\nfai7HmnXdLcY5H6Z5S5d1OpQdz3S4Jpuq2mD/PKXw7p1TT+1qjTrXbqotaHuQ6Rd0/kzyJqkirt0\nUStD3YdIg2s6JwZZ06jiLl3UulD3JdKu6WYYZM2qqrt0UatC3ZdIg2u6bgZZdajyLl3UmlD3KdKu\n6eoYZM1L1XfpolaEuk+RBtf0ahhkNa3qu3RR9qHuW6Rd08szyMpRHXfpoqxD3bdIg2t6D4Ostqjr\nLl2Ubaj7GOk+rmmDrDar8y5dlGWo+xhp6PaaNsjqojrv0kXZhbqvke7KmjbI6ou679JFWYW6r5GG\n9q1pg6w+m8dduiibUPc50jmvaYMsPdW87tJFWYS6z5GGPNa0QZbKmddduqjxUPc90vNe0wZZWr15\n3qWLGg113yMN9a1pgyxVa9536aLKQh0RhwJXAAcBu4F/Syl9asLrUkrJSFPNmjbIUv2auEsXVRnq\ng4GDU0p3RMR+wG3A21NKdy16Xdq8OfU+0jDdmjbIUnOauEsX1Xb6iIj/Av41pfTNRT+fDjoo9T7S\nS61pgyzlpam7dFEtoY6Iw4Eh8IqU0kOL/lm66qrU60gDfOhD8POfwxlnGGQpV03epYsqD/X47DEE\nPpxS+vKEf56guv8v5F1w4IFw5pmjaBtkKQ9N36WLyoZ6bck3tha4CvjspEjvccwxG9m5Ex58EB59\ndMATTwzYay844AB40Yvg2GPhta+F170O9tuv/P8xbbJ9O1xzDXz1q/DDH8InPgGf/ORoTR93HLz1\nrfDud8P++zf9pFI/XXgh7NoFF1ww/197OBwyHA6n/vdKLeqIuAL4bUrpH5Z5TXrHOxJXXw0x/t+H\nBx+EG26Ab3979Cl727bBb38Ljz8Oe+89WpwvecnojyDr18Mb3tC9gO/ePfpAxRe/CDfeCD/7GTz6\nKOy77+iWvX79KNyvfvXC+01SPXK4SxdV+VkfJwPfAX7E6LaRgHNSSl9b9Lp0wgmJ97wHNmxY/hft\ne8B37IDPfW5hdW/fPoq0q1uqTy536aJGvuBl69bEa14D110HJ544/dvoa8Bd3VK9crpLFzX2lYlX\nXw1nnw2bN8PznlfJm+5lwF3dUnWa/nzppTT6JeQbNsCvfsVT7tV16FPAXd3S6uR2ly5qNNSPPQan\nnEKpe3Ud+hJwV7e0vBzv0kWN/+1527Yx0726Dl0PuKtbWpDrXbqo8VADtdyr69DlgLu61Ve53qWL\nsgg1zO9eXYcuBry4ujdtgq1bXd3qnpzv0kXZhLrpe3UduhZwV7e6JPe7dFE2oYY879V1mCbggwG8\n/vV5BtzVrbZqw126KKtQQ3vu1XXoQsBd3WqDNtyli7ILNbT7Xl2HNgfc1a3ctOUuXZRlqLt4r65D\nWwPu6lZT2nSXLsoy1NCfe3Ud2hZwV7fmoW136aJsQw39vlfXoU0Bd3Wram27SxdlHWrwXj0PbQi4\nq1uzaONduij7UHuvbk7uAXd1q4y23qWLsg81eK/OTa4Bd3VrsTbfpYtaEWrwXt0GOQbc1d1vbb5L\nF7Um1OC9uq1yCniZ1X3GGaPVrXZr+126qFWh9l7dLbkEfLnVffzxoz86u7rbpQt36aJWhRq8V/dB\n0wF3dbdbV+7SRa0LNXiv7qsmA+7qbo+u3KWLWhlq8F6tBU0E3NWdpy7dpYtaG2rv1VrJvAPu6m5W\n1+7SRa0NNXiv1urMK+Cu7vnp4l26qNWhBu/Vqs48Ar5jB1x5JXzta67uKnXxLl3U+lCD92rVq86A\nu7pn19W7dFEnQu29Wk2oK+Cu7vK6fJcu6kSowXu18lF1wF3dk3X9Ll3UmVCD92rlrcqAu7q7f5cu\n6lSowXu12qdMwI84YvSXSC0V8L6t7j7cpYs6F2rv1eqKWQPe1dXdl7t0UedCDd6r1W2rDXgXVnef\n7tJFnQw1eK9W/6wm4I880q7V3ae7dFFnQw3eqyWYLuDr148CeN11+a3uvt2lizodau/V0tLKBvyI\nI+APf4Df/Aa2bGlmdffxLl3U6VCD92ppWmUCvt9+o3g//DDs2gW//319q7uvd+miykIdEZcDbwW2\np5SOXuZ1cw01eK+WqrBcwNeuHa3sJ54YhRVgzZpqVndf79JFVYb6FOAh4IrcQg3eq6W6LBfwxdas\ngaOOGn0gs8zq7vNduqjS00dEvBC4NsdQe6+W5mtxwO++e3RrnuQVr4CzzoIzz1xY3X2/Sxf1JtTg\nvVrKwZ6AD4dw/fVwzz2TX/fMZ8IHPgAXXjjXx8tSI6E+77zznvzxYDBgMBiUe9oKeK+W8rJrF9xy\nyyjel18+Opvs8dhjow9g9s1wOGQ4HD754/PPP78/i3oP79VSM1IanUC+//2Fb1u3jj6H+6STFr4d\nckjTT5qXqhf14YxC/cplXtN4qL1XS/OxZy3vifLNN8Ozn/3UKB9zTD9X8zSq/KyPK4EB8CfAduC8\nlNJ/THhd46EG79VS1VzL9en8F7wsx3u1tHqu5fnpdajBe7VUhmu5Wb0Ptfdq6elcy3npfajBe7X6\nzbWcP0M95r1afeFabh9DXeC9Wl3jWu4GQ13gvVpt51ruJkO9iPdqtYVruT8M9QTeq5Uj13J/Geol\neK9Wk1zLKjLUS/BerXlyLWs5hnoZ3qtVB9eypmWoV+C9WrNyLWtWhroE79Uqy7WsOhjqErxXaymu\nZc2DoS7Je7Vcy2qKoZ6C9+p+cS0rF4Z6St6ru8m1rJwZ6il5r+6GPWv5pptGUb7pJtey8mWoV8F7\ndbu4ltV2hnqVvFfny7WsrjHUM/Be3TzXsvrAUM/Ae/X8uZbVR4Z6Rt6r6+NalkYMdQW8V1fDtSxN\nZqgr4r16Oq5lqTxDXRHv1ctzLUurZ6gr5L16xLUsVctQV6yP92rXslQvQ12DLt+rXcvS/BnqGnTp\nXu1alppnqGvSxnu1a1nKk6GuUe73atey1A6Guma53Ktdy1J7GeqaNXWvdi1L3WGo56Due7VrWeq2\nSkMdEW8G/hlYA1yeUvrYhNf0LtRQ7b3atSz1S2Whjog1wN3A64H/A24Fzkgp3bXodb0MNTz9Xj0c\nDhkMBsv+O31Yy2XeD33h+2KB74sFZUO9psTbOhG4J6X0y5TS48DngbfP+oBdcuGF8Otfw6c+Nfrx\ncDh82mt27YIbboALLoDTToMDDoA3vQm+/nU46ii4/HJ44AHYtAk+/nF45zvbHWmY/H7oK98XC3xf\nTG9tidf8KfC/hR//mlG8NbbPPvCFL4zu1SedNFrLP/3p0mv5fe+Dyy5rf4glzUeZUKuEF78YLrkE\nTj0VHn4YPvOZhfPFWWd5W5a0emVu1K8BNqaU3jz+8YeAtPgDihHRzwO1JM2gqg8mPgP4KaMPJt4L\n3AL8ZUrpJ1U8pCRpeSuePlJKT0TE3wLXs/DpeUZakuaksi94kSTVo8yn5y0rIt4cEXdFxN0R8Y9V\nPFQbRcTlEbE9Iv6n6WdpWkQcGhHfiog7I+JHEfF3TT9TUyJin4i4OSJuH78/PtL0MzUtItZExOaI\nuKbpZ2lSRPwiIn44/r1xy7KvnWVRl/1imD6IiFOAh4ArUkpHN/08TYqIg4GDU0p3RMR+wG3A2/v4\n+wIgItallB4Zf7znRuDslNKNTT9XUyLi74FXAc9JKb2t6edpSkRsA16VUnpgpdfOuqj9YpixlNIm\nYMV3eB+klO5LKd0x/v5DwE8YfT5+L6WUHhl/dx9G/53r7e+TiDgUOBW4rOlnyUBQssGzhnrSF8P0\n9r+QerqIOBw4Fri52SdpzviP+rcD9wHDlNKWpp+pQRcBHwT84NjoffCNiLg1It6/3AtnvlFLSxmf\nPa4CNoyXdS+llHanlI4DDgX+LCLWN/1MTYiI04Dt4z9txfhbn52cUjqe0Z8w/mZ8Pp1o1lD/Bjis\n8ONDxz+nnouItYwi/dmU0pebfp4cpJR2AV8BTmj6WRpyMvC28W32c8CfR8QVDT9TY1JK947/837g\nSyzzV3PMGupbgSMi4oURsTdwBtDnj+S6Ehb8O7AlpfQvTT9IkyLigIh47vj7zwL+Arij2adqRkrp\nnJTSYSmlFzNqxbdSSu9t+rmaEBHrxn/iJCL2Bd4I/Hip188U6pTSE8CeL4a5E/h8X78YJiKuBL4H\nvDQifhURf930MzUlIk4G/gp43fhTjzaP/07zPjoE+O/xjfom4JqU0jcbfiY17yBgU+H3xbUppeuX\nerFf8CJJmfODiZKUOUMtSZkz1JKUOUMtSZkz1JKUOUMtSZkz1JKUOUMtSZn7fwMyXhtP790WAAAA\nAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Triangulate a polygon\n", "#\n", "# The python function \"polygon_triangulate.py\" should be available on\n", "# the webpage http://people.sc.fsu.edu/~jburkardt/classes/urop_2016/urop_2016.html\n", "#\n", "# Download that file, and use an import statement to access the main function,\n", "# which has the form:\n", "# t = polygon_triangulate ( n, x, y )\n", "# where \n", "# n is the number of polygon vertices, \n", "# x and y are the coordinates of the vertices,\n", "# t is the [n-2,3] array of vertex indices forming triangles.\n", "#\n", "# The array T describes a way of decomposing the polygon into triangles\n", "# all of which have positive area. \n", "#\n", "# Now draw the polygon by drawing the individual triangles.\n", "# To plot a triangle, you need to make lists of the x and y coordinates\n", "# with an extra copy of the first vertex tacked on to the end of the list.\n", "#\n", "# Looking at the plot, you should realize that all the triangles\n", "# have positive area.\n", "#\n", "from polygon_triangulate import polygon_triangulate\n", "\n", "#\n", "# polygon_triangulate wants a vector of x's and a vector of y's.\n", "#\n", "x = [ a[0], b[0], c[0], d[0], e[0], f[0] ]\n", "y = [ a[1], b[1], c[1], d[1], e[1], f[1] ]\n", "t = polygon_triangulate ( 6, x, y )\n", "\n", "print ( t )\n", "\n", "for i in range ( 0, 4 ):\n", " x4 = np.array ( [ x[t[i,0]], x[t[i,1]], x[t[i,2]], x[t[i,0]] ] )\n", " y4 = np.array ( [ y[t[i,0]], y[t[i,1]], y[t[i,2]], y[t[i,0]] ] )\n", " plt.plot ( x4, y4, 'b-' )\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Use the positive area triangle triangulation\n", "#\n", "# Using the triangulation computed by polygon_triangulate,\n", "# recompute the triangle areas, and the polygon arean.\n", "#\n", "# Then recompute the polygon centroid.\n", "#\n", "# We should get the same total area, and centroid as we did\n", "# when one triangle had negative area. It's nice when occasionally\n", "# a weird thing like a negative area triangle pops up, but doesn't\n", "# hurt your calculation.\n", "#\n", "# However, for the next steps, point inclusion and random sampling, \n", "# we definitely need to use only positive area triangles!\n", "#\n", "from polygon_triangulate import polygon_triangulate\n", "#\n", "# polygon_triangulate wants a vector of x's and a vector of y's.\n", "#\n", "x = [ a[0], b[0], c[0], d[0], e[0], f[0] ]\n", "y = [ a[1], b[1], c[1], d[1], e[1], f[1] ]\n", "t = polygon_triangulate ( 6, x, y )\n", " \n", "polygon_area = 0.0\n", "polygon_centroid = np.zeros ( 2 )\n", "for i in range ( 0, 4 ):\n", " v0 = np.array ( [ x[t[i,0]], y[t[i,0]] ] )\n", " v1 = np.array ( [ x[t[i,1]], y[t[i,1]] ] )\n", " v2 = np.array ( [ x[t[i,2]], y[t[i,2]] ] )\n", " area = triangle_area ( v0, v1, v2 )\n", " polygon_centroid = polygon_centroid + area * ( v0 + v1 + v2 ) / 3.0\n", " print ( '%d %d %d %f' % ( t[i,0], t[i,1], t[i,2], area ) )\n", " polygon_area = polygon_area + area\n", "\n", "print ( 'Polygon area = %f' % ( polygon_area ) )\n", "polygon_centroid = polygon_centroid / polygon_area\n", "print ( 'Polygon centroid = %f %f' % ( polygon_centroid[0], polygon_centroid[1] ))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Centroid is inside triangle 1\n" ] } ], "source": [ "# Is a point inside a polygon?\n", "#\n", "# Assuming we have decomposed our polygon into a collection of triangles\n", "# and none of the triangles have negative area, then this problem is easy to solve.\n", "#\n", "# If polygon P consists of N triangles T1, T2, ..., TN, then the point X is\n", "# in P if and only if it is inside of at least one triangle Ti.\n", "#\n", "# We already know how to determine if a point is inside a triangle (refer to\n", "# the triangle notebook!)\n", "#\n", "# We know that a triangle always contains its centroid. This is not necessarily\n", "# true for a polygon, (can you think of an example?) However, let's ask whether\n", "# the centroid of our example polygon is contained in the polygon.\n", "#\n", "from polygon_triangulate import polygon_triangulate\n", "#\n", "# polygon_triangulate wants a vector of x's and a vector of y's.\n", "#\n", "x = [ a[0], b[0], c[0], d[0], e[0], f[0] ]\n", "y = [ a[1], b[1], c[1], d[1], e[1], f[1] ]\n", "t = polygon_triangulate ( 6, x, y )\n", " \n", "polygon_area = 0.0\n", "polygon_centroid = np.zeros ( 2 )\n", "for i in range ( 0, 4 ):\n", " v0 = np.array ( [ x[t[i,0]], y[t[i,0]] ] )\n", " v1 = np.array ( [ x[t[i,1]], y[t[i,1]] ] )\n", " v2 = np.array ( [ x[t[i,2]], y[t[i,2]] ] )\n", " area = triangle_area ( v0, v1, v2 )\n", " polygon_centroid = polygon_centroid + area * ( v0 + v1 + v2 ) / 3.0\n", " polygon_area = polygon_area + area\n", "\n", "polygon_centroid = polygon_centroid / polygon_area\n", "\n", "inside = False\n", "for i in range ( 0, 4 ):\n", " v0 = np.array ( [ x[t[i,0]], y[t[i,0]] ] )\n", " v1 = np.array ( [ x[t[i,1]], y[t[i,1]] ] )\n", " v2 = np.array ( [ x[t[i,2]], y[t[i,2]] ] )\n", " area_c12 = triangle_area ( polygon_centroid, v1, v2 )\n", " area_0c2 = triangle_area ( v0, polygon_centroid, v2 )\n", " area_01c = triangle_area ( v0, v1, polygon_centroid )\n", " if ( 0 <= area_c12 and 0.0 <= area_0c2 and 0.0 <= area_01c ):\n", " inside = True\n", " print ( 'Centroid is inside triangle %d' % ( i ) )\n", " break\n", "if ( not inside ):\n", " print ( 'Centroid is not inside the polygon.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Uniform sampling of a line divided into unequal lengths #\n", "\n", "Before we can consider uniform sampling of polygons, it will be helpful\n", "to consider how to handle the following problem:\n", " \n", "Suppose a line is divided up into N intervals I1, I2, ..., IN of \n", "lengths L1, L2, ..., LN. Then consider the partial sums of the\n", "lengths, S0 = 0, S1 = S0 + L1, S2 = S1 + L2, ..., SN = SN-1+LN.\n", "\n", "Normalize by dividing by SN: P0 = 0, P1 = S1/SN, P2 = S2/SN, ..., PN=SN/SN=1.\n", " \n", "One way to sample the line is to choose an interval and then choose a\n", "point in the interval. If the intervals are of different lengths,\n", "then we need to choose longer intervals more often, and in proportion\n", "to their lengths. To do that, choose a random number r1 between 0 and 1, and then \n", " choose I1 if S0 <= r1 < S1\n", " I2 if S1 <= r1 < S2\n", " ...\n", " IN if SN-1 <= r1 <= SN \n", "Once we have decided on the interval uniformly, we then have to pick a \n", "point inside the interval uniformly. But that just means that if we\n", "have decided to sample interval I7, we need a random value r2\n", "which we then have to linearly transform between S6 and S7.\n", " s = S6 + (S7-S6)*r2 = S6 + L7 * r2." ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[5 1 4]\n", "[ 0. 5. 6. 10.]\n", "[ 0. 0.5 0.6 1. ]\n", "[ 503. 115. 382.]\n" ] } ], "source": [ "# Uniform sampling of a line divided into unequal lengths\n", "#\n", "# Consider the line segment from 10 to 20, divided into the \n", "# subintervals [10,15], [15,16], [16,20].\n", "#\n", "# If we choose 100 points uniformly, we expect about 50 to be in\n", "# the first segment, about 10 in the second, and about 40 in the third.\n", "#\n", "# Put together a code to compute 1000 sample points from this line segment,\n", "# using the two step process in which you first choose which subinterval\n", "# to sample, and then choose a point in that subinterval.\n", "#\n", "# Normally, we'd be interested in the final X values, but for this example,\n", "# we mainly want to show that we can sample uniformly from unevenly sized\n", "# intervals if we know how big they are.\n", "#\n", "l = np.array ( [ 5, 1, 4 ] )\n", "print ( l )\n", "\n", "s = np.zeros ( 4 )\n", "for i in range ( 1, 4 ):\n", " s[i] = s[i-1] + l[i-1]\n", "print ( s )\n", "\n", "p = np.zeros ( 4 )\n", "p = s / s[3]\n", "print ( p )\n", "\n", "count = np.zeros ( 3 )\n", "n = 1000\n", "for j in range ( 0, n ):\n", " r1 = np.random.random ( )\n", " for i in range ( 1, 4 ):\n", " if ( r1 <= p[i] ):\n", " interval = i - 1\n", " break\n", " r2 = np.random.random ( )\n", " x = s[interval] + l[interval] * r2\n", " count[interval] = count[interval] + 1\n", "#\n", "# Print count to verify roughly 500 / 100 / 400 split in interval choice.\n", "#\n", "print ( count )\n", " " ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Uniform sampling of a polygon #\n", "\n", "Suppose our region is a polygon. As discussed in the polygon notebook,\n", "any polygon can be represented as a collection of triangles.\n", "\n", "For sampling, we must insist that all the triangles have positive\n", "(well, nonnegative) area, so to be safe we should triangulate using\n", "a function like polygon_triangulate().\n", "\n", "Now if we wish to sample uniformly from a polygon P, we could \n", "use a two step process:\n", "* decide which triangle Ti to sample;\n", "* decide which point xj in Ti to sample.\n", "\n", "We know how to choose x in T uniformly, because T is a triangle\n", "and we did that in the triangle notebook.\n", "\n", "So how do we do step 1 so that the process is uniform?\n", "All we have to do is choose triangles with a frequency proportional\n", "to their area, since that \"counts\" how many points are in that triangle.\n", "But that is the same as the problem of sampling unequal line segments uniformly.\n", "\n", "In other words, we prepare for step 1 by computing the area Ai \n", "of each triangle Ti, creating an increasing sequence of numbers P \n", "which are spaced proportionally to the triangle areas, select a random\n", "number r1, and use that to select triangle Ti." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.33333333 1. ]\n", "\n", " IT Hits Area CentroidX CentroidY\n", "\n", " 0 342 0.513000 0.359181 0.329075\n", " 1 658 0.987000 1.011538 0.672720\n" ] } ], "source": [ "# Uniform sampling of a polygon #\n", "#\n", "# We'll use a simple quadrilateral made up of two triangles, \n", "# one having twice the area of the other.\n", "#\n", "# 1 (3)---*---(2)\n", "# | | \\ T1 / \n", "# Y | \\ /\n", "# | |T0 \\ /\n", "# 0 (0)--(1)\n", "# \n", "# 0----1-----2---X\n", "#\n", "# Prepare:\n", "#\n", "# Define a list XY of 4 point coordinates.\n", "# Define the polygon P as two triangles, that is, 2 rows of 3 indices into XY.\n", "# Compute the areas of T0, T1 and P.\n", "# Compute the array of S values.\n", "#\n", "# Sample:\n", "#\n", "# Compute 1000 random samples of the polygon.\n", "# Estimate the area and centroid of each triangle.\n", "# (We can't compute the energy, because we didn't \n", "# pick generator points in each triangle)\n", "#\n", "xy = np.array ( [ [0.0,0.0], [1.0,0.0], [2.0,1.0], [0.0,1.0] ] )\n", "\n", "p = np.array ( [ [ 0, 1, 3 ], [ 1, 2, 3 ] ] )\n", "\n", "a0 = 0.5 * ( xy[p[0,0],0] * ( xy[p[0,1],1] - xy[p[0,2],1])\\\n", " + xy[p[0,1],0] * ( xy[p[0,2],1] - xy[p[0,0],1])\\\n", " + xy[p[0,2],0] * ( xy[p[0,0],1] - xy[p[0,1],1]) )\n", " \n", "a1 = 0.5 * ( xy[p[1,0],0] * ( xy[p[1,1],1] - xy[p[1,2],1])\\\n", " + xy[p[1,1],0] * ( xy[p[1,2],1] - xy[p[1,0],1])\\\n", " + xy[p[1,2],0] * ( xy[p[1,0],1] - xy[p[1,1],1]) )\n", "\n", "ap = a0 + a1\n", "\n", "s = np.array ( [ a0 / ( a0 + a1 ), ( a0 + a1 ) / ( a0 + a1 ) ] )\n", "\n", "hit = np.zeros ( 2 )\n", "centroid = np.zeros ( [ 2, 2 ] )\n", "\n", "n = 1000\n", "\n", "for j in range ( 0, n ):\n", "#\n", "# Select triangle IT.\n", "#\n", " r0 = np.random.random ( )\n", " if ( r0 < s[0] ):\n", " it = 0\n", " else:\n", " it = 1\n", " \n", " hit[it] = hit[it] + 1\n", "#\n", "# Select point S in triangle IT.\n", "#\n", " r1 = np.random.random ( )\n", " r2 = np.sqrt ( np.random.random ( ) )\n", "\n", " i = 1 - r2\n", " j = ( 1 - r1 ) * r2\n", " k = r1 * r2\n", " sample = i * xy[p[it,0],:] + j * xy[p[it,1],:] + k * xy[p[it,2],:]\n", " \n", " centroid[it,:] = centroid[it,:] + sample\n", "\n", "print ( '' )\n", "print ( ' IT Hits Area CentroidX CentroidY' )\n", "print ( '' )\n", "for it in range ( 0, 2 ):\n", " centroid[it,:] = centroid[it,:] / float ( hit[it] )\n", " area_est = ap * float ( hit[it] )/ float ( n )\n", " print ( ' %2d %4d %14.6f %14.6f %14.6f' % ( it, hit[it], area_est, centroid[it,0], centroid[it,1] ) )\n", "#\n", "# The correct areas are 0.5 and 1.0\n", "#" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 1000 9.325161\n", " 10000 9.546162\n", " 100000 9.553896\n" ] } ], "source": [ "# Integral estimate over a polygon by sampling.\n", "#\n", "# To approximate the integral of f(x,y) over our quadrilateral polygon using N points,\n", "# we write\n", "# I approximately sum ( f(x,y) ) * area / N\n", "# where we evaluate f(*,*) at N points chosen uniformly at random.\n", "#\n", "# Our function will be f(x,y) = x^2 * exp ( x + y )\n", "#\n", "# Estimate the integral using 1000, 10000, and 100000 samples.\n", "#\n", "xy = np.array ( [ [0.0,0.0], [1.0,0.0], [2.0,1.0], [0.0,1.0] ] )\n", "\n", "p = np.array ( [ [ 0, 1, 3 ], [ 1, 2, 3 ] ] )\n", "\n", "a0 = 0.5 * ( xy[p[0,0],0] * ( xy[p[0,1],1] - xy[p[0,2],1])\\\n", " + xy[p[0,1],0] * ( xy[p[0,2],1] - xy[p[0,0],1])\\\n", " + xy[p[0,2],0] * ( xy[p[0,0],1] - xy[p[0,1],1]) )\n", " \n", "a1 = 0.5 * ( xy[p[1,0],0] * ( xy[p[1,1],1] - xy[p[1,2],1])\\\n", " + xy[p[1,1],0] * ( xy[p[1,2],1] - xy[p[1,0],1])\\\n", " + xy[p[1,2],0] * ( xy[p[1,0],1] - xy[p[1,1],1]) )\n", "\n", "area_polygon = a0 + a1\n", "\n", "s = np.array ( [ a0 / ( a0 + a1 ), ( a0 + a1 ) / ( a0 + a1 ) ] )\n", "\n", "hit = np.zeros ( 2 )\n", "\n", "n = 1000\n", "\n", "for step in range ( 0, 3 ):\n", "\n", " estimate = 0.0\n", "\n", " for j in range ( 0, n ):\n", "#\n", "# Select triangle IT.\n", "#\n", " r0 = np.random.random ( )\n", " if ( r0 < s[0] ):\n", " it = 0\n", " else:\n", " it = 1\n", " \n", " hit[it] = hit[it] + 1\n", "#\n", "# Select point S in triangle IT.\n", "#\n", " r1 = np.random.random ( )\n", " r2 = np.sqrt ( np.random.random ( ) )\n", "\n", " i = 1 - r2\n", " j = ( 1 - r1 ) * r2\n", " k = r1 * r2\n", " sample = i * xy[p[it,0],:] + j * xy[p[it,1],:] + k * xy[p[it,2],:]\n", " \n", " estimate = estimate + sample[0] ** 2 * np.exp ( sample[0] + sample[1] )\n", "\n", " estimate = estimate * area_polygon / float ( n )\n", "\n", " print ( ' %5d %f' % ( n, estimate ) )\n", " n = n * 10" ] }, { "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 }