2c:null 2d:["$","title",null,{"children":"Understanding point-in-polygon | Tom MacWright | Observable"}] 32:Tbeb,md`# Understanding point-in-polygon In geometry and geography, it’s very often necessary to know whether a given point falls within a polygon. [point-in-polygon](https://github.com/substack/point-in-polygon), by [James Halliday](https://github.com/substack), has become the standard algorithm for figuring that out - within the context of JavaScript. It’s a cool algorithm because it’s incredibly compact, but fairly deep. Here's the full source: \`\`\`js module.exports = function (point, vs) { // ray-casting algorithm based on // http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html var x = point[0], y = point[1]; var inside = false; for (var i = 0, j = vs.length - 1; i < vs.length; j = i++) { var xi = vs[i][0], yi = vs[i][1]; var xj = vs[j][0], yj = vs[j][1]; var intersect = ((yi > y) != (yj > y)) && (x < (xj - xi) * (y - yi) / (yj - yi) + xi); if (intersect) inside = !inside; } return inside; }; \`\`\` The gist of this solution to the [the Point in Polygon](https://en.wikipedia.org/wiki/Point_in_polygon) problem is: 1. Cast a ray (essentially ‘draw a line’) from the left side of the drawing until it hits the point. 2. Count the number of edges it crosses before it gets to the point 3. If that number is odd, then the point is in the polygon. If it's even (including 0) it's not in the polygon. It links to the [PNPoly](https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html) webpage, which then refers to [Simulation of Simplicity](https://arxiv.org/abs/math/9410209), a 1994 paper that explains a very general approach for _coping with degenerate input data_. > “Don't mess with this code unless you're familiar with the idea of Simulation of Simplicity.” PNPoly was authored by [W. Randolf Franklin](https://wrf.ecse.rpi.edu/pmwiki/pmwiki.php), a professor at RPI. If you've dealt with geometry before, you'll immediately understand that degenerate data is the cause of many bad times. I’m going to work through the algorithm, taking notes as I go. Writing this, I don’t understand the algorithm yet, though I’ve used it so many times. So, I’ll start with the structure of the code. There’s - x & y, the point coordinates - an inside boolean, that flips each time the ray hits a line - the loop: what kind of loop is this? it’s an interesting one ### Values of i & j in each iteration Upon inspection, well - it makes sense. The for loop iterates through each edge of the polygon, including the one from the last coordinate to the first. Note that, funny enough, that means that [GeoJSON polygons](https://tools.ietf.org/html/rfc7946#section-3.1.6) and all other geometries that repeat a coordinate at the first and last position are wasting 1 iteration of this loop - because the first iteration, where you're comparing the first and last position (0 and 9 here), will be comparing the same coordinate. Think of the processor time, in total, wasted! It’s tiny, but anyway. `33:T5fb,md`## The second condition \`x < (xj - xi) * (y - yi) / (yj - yi) + xi\` Is tougher to follow. The first level - for me, because despite _knowing_ operator precedence by the book, it doesn’t come naturally - is what order it should be read in. Cluttering it up with lots of parentheses to remove all idea of precendence leaves us with this: \`x < ((((xj - xi) * (y - yi)) / (yj - yi)) + xi)\` So this code is trying to determine: **Is the point in the half-plane to the left of the extended edge?** Buckle up, it's going to get much more difficult from here on out. So, on page 22 of Simulation of Simplicity, this is one of the usage examples. The first paragraph is a review of how point-in-polygon tests work: all good so far. > To perturb the points we use the same ${tex`\epsilon`} expansion as described in Section 3.2, that is, we replace ${tex`v _ { i } = \left( \nu _ { i ,1} ,\nu _ { i ,2} \right)`} by ${tex`v _ { i } ( \epsilon ) = \left( \nu _ { i ,1} ( \epsilon ) ,\nu _ { i ,2} ( \epsilon ) \right)`} where ${tex`\nu _ { i ,j } ( \epsilon ) = \nu _ { i ,j } + \epsilon ( i ,j )`} with ${tex`\in ( i ,j )`} as in 3-b. Anyway, I’ll meditate on that for a bit and do a quick demonstration of what the half-plane to the left of the extended edge means: The dotted lines in the figure below are the _extended edges_ of each edge that has an effect on the point-in-polygon result. Click the numbers to choose which half-plane to show. The half-planes are areas extending from the extended edge to the left. `34:T7dd,{ let relevantEdges = []; for (let i = 0, j = shape.length - 1; i < shape.length; j = i++) { if ((shape[i][1] > point[1]) != (shape[j][1] > point[1])) { relevantEdges.push([shape[i], shape[j]]); } } let extendedEdges = []; for (let edge of relevantEdges) { // ‘rise over run’ let m = (edge[0][1] - edge[1][1]) / (edge[0][0] - edge[1][0]); let b = edge[0][1] - (m * edge[0][0]); extendedEdges.push([[-b / m, 0], [(320 - b) / m, 320]]); } var halfPlane = 0; while (true) { let ui = html`

`; yield ui; halfPlane = await new Promise(resolve => { ui.addEventListener('click', e => { if (e.target.tagName === 'rect') { resolve(e.target.getAttribute('data-value')); } }) }); } }35:T610,{ let relevantEdges = []; for (let i = 0, j = shape.length - 1; i < shape.length; j = i++) { if ((shape[i][1] > point[1]) != (shape[j][1] > point[1])) { relevantEdges.push([shape[i], shape[j]]); } } let extendedEdges = []; for (let edge of relevantEdges) { // ‘rise over run’ let m = (edge[0][1] - edge[1][1]) / (edge[0][0] - edge[1][0]); let b = edge[0][1] - (m * edge[0][0]); extendedEdges.push([[-b / m, 0], [(320 - b) / m, 320]]); } let intercepts = []; for (let edge of relevantEdges) { let xi = edge[0][0], xj = edge[1][0], yi = edge[0][1], yj = edge[1][1]; let y = point[1]; intercepts.push([(xj - xi) * (y - yi) / (yj - yi) + xi, point[1]]); } return html`

`; }36:T477,md`Yep, those are intercepts. But I’m in detective mode and should switch back to analysis. So, connecting the dots between this code, the paper it references, and my relatively light sense of geometry, I re-read the Wikipedia page for this intercept problem, and the [2D vector determinant form is this one](https://en.wikipedia.org/wiki/Linear_equation#2D_vector_determinant_form). The 2D vector determinant form, then, is ${tex.block`\left(x - x_1 \right) \left(y_2 - y_1 \right) = \left(y - y_1 \right) \left(x_2 - x_1 \right)`} _note to self, this is getting hairy and I need to still draw more connections..._ ${tex.block`\frac{(x - x_1)(y_2 - y_1)}{y - y_1} + x_1 = x_2`} Finally! The ‘solve for x’ of the 2D determinant form is clearly the same thing as (xj - xi) * (y - yi) / (yj - yi) + xi Now the mystery is: does the author of PNPoly actually use Simulation of Simplification in a way that changes the implementation of this algorithm? From my naïve understanding, this is roughly a direct translation from textbook math into code. Anyway, I think we should have one final diagram to demonstrate the algorithm! `37:T803,{ let [x, y] = point; var inside = false; var crossings = 0; await visibility(); restart; for (var i = 0, j = shape.length - 1; i < shape.length; j = i++) { var xi = shape[i][0], yi = shape[i][1]; var xj = shape[j][0], yj = shape[j][1]; var intersect = ((yi > y) != (yj > y)) && (x < (xj - xi) * (y - yi) / (yj - yi) + xi); if (intersect) { inside = !inside; crossings++; } let edgeRelevant = ((yi > y) != (yj > y)); let intercept, extendedEdge; if (edgeRelevant) { let y = point[1]; intercept = [(xj - xi) * (y - yi) / (yj - yi) + xi, point[1]]; let m = (yi - yj) / (xi - xj); let b = yi - (m * xi); extendedEdge = [[-b / m, 0], [(320 - b) / m, 320]]; } console.log(`M ${[shape[i], shape[j]].map(s => s.join(' ')).join('L')}`); yield html`

edge: ${edgeRelevant ? 'relevant' : 'irrelevant'}
crossings: ${crossings}
current state: ${inside ? 'inside (odd count)' : 'outside (even count)'}

` await Promises.delay(800); } }38:T4ca,{/*`I re-read the Wikipedia article about [Linear equations](https://en.wikipedia.org/wiki/Linear_equation#Two-point_form) and thankfully find this diagram: ${tex`\left| \begin{array} { l l l } { x } & { y } & { 1} \\ { x _ { 1} } & { y _ { 1} } & { 1} \\ { x _ { 2} } & { y _ { 2} } & { 1} \end{array} \right| = 0`} For the determinant form of the two-point equation of a line. This matches up with ${tex`\left| \begin{array} { r l } { \nu _ { i ,1} ( \epsilon ) } & { \nu _ { i, 2 } ( \epsilon ) } & { 1} \\ { \nu , _ 1} & { \nu _ { j ,2} } & { 1} \\ { \nu _ { k ,1} } & { \nu _ { k ,2} } & { 1} \end{array} \right|`} From Simulation of Simplicity: it’s the core of their solution the point-in-polygon problem. So both are 3x3 matrices, but the one in Simulation of Simplicity has a little bit of extra magic in it. Anyway, since I need to relearn this: you’d transform the first example back into this: ${tex`\begin{vmatrix} x & y & 1\\x_1 & y_1 & 1\\x_2 & y_2 & 1 \end{vmatrix} = x\,\begin{vmatrix} y_1 & 1\\y_2 & 1 \end{vmatrix} - y\,\begin{vmatrix} x_1 & 1\\x_2 & 1 \end{vmatrix} + 1\,\begin{vmatrix} x_1 & y_1\\x_2 & y_2 \end{vmatrix} = xy_1 - xy_2 - x_1y + x_1y_2 + x_2y - x_2y_1`} Regrouping the terms:`*/}39:[] 30:{"id":"70965ad0643d4783","slug":"understanding-point-in-polygon","title":"Understanding point-in-polygon","thumbnail":"b15c7aba1f686865f36a4248c59bd5bd2307d7ae94404bb1f400a42344886f21","access_level":"public","create_time":"$D2018-03-08T00:49:48.971Z","update_time":"$D2021-02-12T21:23:02.909Z","view_time":null,"first_public_time":"$D2018-03-09T01:46:03.170Z","has_unpublished_edits":null,"fork":false,"stars":41,"starred":false,"owner":{"type":"team","id":"b47c3705cfcd7c1b","login":"tmcw","name":"Tom MacWright","avatar_url":"https://avatars-next.observableusercontent.com/669bdd774da8a1b9c1b7706ecbc1ea5136b9932c3d5e86b753aa5a098005f674","tier":"starter_2024"},"version":763,"publish_version":763,"role":null,"authored":false,"authors":[{"type":"individual","id":"0bf14b8d9c6eb88a","login":"tmcw","name":"Tom MacWright","avatar_url":"https://avatars-next.observableusercontent.com/669bdd774da8a1b9c1b7706ecbc1ea5136b9932c3d5e86b753aa5a098005f674","tier":null}],"members":[{"type":"individual","id":"0bf14b8d9c6eb88a","login":"tmcw","name":"Tom MacWright","avatar_url":"https://avatars-next.observableusercontent.com/669bdd774da8a1b9c1b7706ecbc1ea5136b9932c3d5e86b753aa5a098005f674","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":11,"collection":{"id":"b0b8183de839aa32","slug":"algorithms","title":"Algorithms","owner":{"id":"f35c755083683fe5","login":"observablehq","name":"Observable"}},"type":"notebook","body":{"titling":"auto","title":"Understanding point-in-polygon","cells":[{"id":0,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":15,"value":"{\n let length = 10;\n let values = [];\n for (var i = 0, j = length - 1; i < length; j = i++) {\n values.push(`i = ${i} j = ${j}`);\n }\n return html`

${values.join('\\n')}

`;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":34,"value":"md`Continuing on, the author describes _the method_ as having two tests:\n\n1. Is the point in the half-plane to the left of the extended edge? and\n2. Is the point's Y coordinate within the edge's Y-range?\n\nPicking apart the code, these refer to:\n\n1. \\`x < (xj - xi) * (y - yi) / (yj - yi) + xi\\`\n2. \\`yi > y != yj > y\\`\n\nAt this point, a graphic will be useful.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":60,"value":"svg``","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":158,"value":"md`## The first condition\n\n\\`yi > y != yj > y\\` is asking:\n\nGiven that yi and yj are the y-coordinates of the edge’s ends, is the point between them? It tests that by ask whether the answers to “is this below” are the opposite. If one is below and the other one is above, then the point is between. And if this condition fails, then the edge has no effect on the status of the point being inside or outside.\n\n_Side note: there’s no part of this algorithm that really cares whether you angle the line along the \\`x\\` or \\`y\\` axis. It could be a vertical line instead - people just tend to use horizontal lines._\n\nSo essentially we’re _eliminating_ edges from consideration at this step: in the shape above, only these edges matter, because they are intersect with a line leveled at the point of interest’s y coordinate:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":162,"value":"{\n let relevantEdges = [];\n for (let i = 0, j = shape.length - 1; i < shape.length; j = i++) {\n if ((shape[i][1] > point[1]) != (shape[j][1] > point[1])) {\n relevantEdges.push([shape[i], shape[j]]);\n }\n }\n return svg``\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":240,"value":"md`Already we can run through the algorithm in our minds and see that there are 3 red lines to the right of the red dot. 3 is an odd number, so the dot is inside the polygon. But for the algorithm to make that observation, it will need to do some more work.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":201,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":243,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":444,"value":"md`So, getting back to it, the centerpiece of the paper is this algorithm for determining whether a half-line intersects with an edge. The paper uses the term _half-line_ in the same way that the Wikipedia page uses a _ray_.\n\nPicking this apart:\n\n x < (xj - xi) * (y - yi) / (yj - yi) + xi\n\nWe’re testing whether \\`x\\` is less than the value, so we can assume that, essentially, the quantity on the right is supposed to be the x-coordinate of a point on the edge at the same y coordinate as the intersection point. I’ll test that theory by adding red dots at all the possible values of the right side of this inequality.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":473,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":501,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":665,"value":"$37","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":742,"value":"viewof restart = html``","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":134,"value":"point = [300, 120]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":127,"value":"shape = [\n [100, 110],\n [130, 50],\n [160, 150],\n [180, 140],\n [190, 30],\n [240, 250],\n [250, 240],\n [290, 20],\n [350, 50],\n [370, 180],\n [380, 190],\n [390, 10],\n [420, 10],\n [430, 150],\n [400, 200],\n [360, 200],\n [300, 260],\n [250, 280],\n [220, 240],\n [200, 260],\n [120, 240]\n ]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":651,"value":"$38","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W39","stdlib":"1"}} 31:"$30:body"