31:T431,We are working towards enabling more flexible georeference transformations such as **Thin Plate Spline transformations** in the [Allmaps](https://allmaps.org/) project. To make this happen we need to do three things: 1. Write an algorithms that break our map images (more specifically the area inside image mask) up into small **triangles**. This is a specific case of 'polygon triangulation'. 👈 **This is what this notebook is about.** 2. Write some functions that, given a set of ground controle points, can **transform** any image coordinate to it's corresponding map coordinate. (Currently only the standard 'affine' transformation is implemented in Allmaps, which translates, rotates and shears the image but does not allow more flexible reprojections). For this, see our [next notebook](https://observablehq.com/d/0b57d3b587542794) in this series. 3. Implement both if these in our **render** code (tell WebGL to draw each triangle of the map image in a specific location), and add this functionality to the viewer. So, let's make a triangulation algorithm!32:T850,### Polygon triangulation, you say? We are looking for an algorithm that divides the masked area of a map (i.e. a polygon) into a set of (small enough) triangles. Such polygon triangulation is a specific case of ['mesh generation'](https://en.wikipedia.org/wiki/Mesh_generation): for 2D shapes and using triangles for the mesh (as opposed to using squares, hexagons, ...). For the similar problem of drawing triangles between *points* the well-known ['Delaunay triangulation'](https://en.wikipedia.org/wiki/Delaunay_triangulation) algorithm exists. It's a beautiful algorithm and has really nice properties. One being that it maximises the minimum angle. In practice this means it delivers relatively well-conditioned triangles (i.e. they tend to look more like isosceles triangles then like fine wedges). We, however, are starting from a *polygon* (or more specifically a planar straight-line graph) which is more complicated, as we want the resulting triangles not only to respect the polygon vertices but also all its edges. Also, we want small enough triangles. So solve this problem we have to turn to generalisations of Delaunay triangulation. One option is to use algorithms for triangulations between points that allow for setting conditions on the resulting triangles, such as respecting certain edges (the polygon edges) or points (extra points we could add in the middle to ensure the triangles are small enough). This is called ['constrained Delaunay triangulation'](https://en.wikipedia.org/wiki/Constrained_Delaunay_triangulation)). Another option would be to use algorithms specifically for polygons that automatically add extra points as long as the triangles don't meet the input requirements on e.g. their size or maximum angle. This are called ['conforming Delaunay triangulation'](https://en.wikipedia.org/wiki/Delaunay_refinement#Ruppert's_algorithm) algorithms, such as Ruppert's algorithm. Finally, there are algorithms that follow other methods, like ear slicing (and constrained versions of those). They can be very fast, but their triangles are typically poorly conditioned: they are long and thin.33:Tb84,### Starting from the right existing algorithm We are looking for an algorithm with the following properties: - Implements **constrained Delaunay triangulation** of points set (and write our own code around it 📝) or **conforming Delaunay triangulation** of polygons (suites our needs out of the box ✨) - Written in **JavaScript (or TypeScript)** and on NPM - **Small** (<500KB including dependencies) - **Robust enough** (also works with weird polygons, but our inputs should be quite simple) - **Fast enough** (our input polygons are small and few, so we don't need record speeds) Here's an overview of some algorithms we found out there. | Name | Algorithm | Language | Size (gziped) | Robustness | Speed | | ----------------------------------------------------- | --------------------------------------- | ----------------------- | -------------------------- | ---------- | ----- | | [Delaunator](https://github.com/mapbox/delaunator) | Delaunay ❌ | JS | 8.1 KB | ++ | ++ | | [Earcut](https://github.com/mapbox/earcut) | Constrained Ear-slicing triangulation | JS | 6.5 KB | + | +++ | | [libtess](https://github.com/brendankenny/libtess.js) | Polygon triangulation (not conforming?) ❌ | JS (port from OpenGL) | 15 KB | +++ | + | | [cdt-js](https://github.com/savithru-j/cdt-js) | Constrained Delaunay triangulation | JS (demo site, not NPM) ❌ | 723 KB (unzipped) | ? | ? | | [Triangle WASM](https://github.com/brunoimbrizi/triangle-wasm) | Conforming Delaunay triangulation | JS (WebAssembly from C) 📦 | 12.8 KB | +++ | ++ | | [poly2tri.js](https://github.com/r3mi/poly2tri.js) | Constrained Delaunay triangulation | TS (port from C++) | 19.3 KB | + | + | (See the earcut readme for [performance benchmarks](https://github.com/mapbox/earcut) on standard Delaunay triangulation) There's also a [Triangulations](https://github.com/mkacz91/Triangulations) project on GitHub, an interactive website to explain and demo some triangulation algorithms, including conforming Delaunay triangulation - nice! The code used in the demos is written in JS but not structured as a library, and said to have sub-obtimal robustness. Maybe it can serve as a good start we want to take a different approach later. Since we want a triangulation in with relatively well conditioned triangles (not long and thin) we won't pick the Earcut algorithm. Since we prefer to work with plain JS/TS code over WebAssembly, we don't pick Triangle WASM. 👉 So **we'll build our polygon triangulation algorithms around poly2tri.js**! 34:T5a7,### What our algorithm does Goal: > Given a polygon and a grid size, we want to get a set of triangles, smaller then the grid size and triangulating the polygon. What our code does: > It takes the incoming polygon, breaks it's edges into smaller pieces (of the 'grid size' length), builds a grid in the bounding box of the polygon, and gives poly2tri the interpolated polygon edges and the grid points that which are inside the polygon, as Steiner points, to assure the resulting triangles are small enough. Then poly2tri.js does the heavy lifting of the actual triangulation. Many thanks to it's original poly2tri developers and the community who ported it! Possible improvements: - Assure it work for all (or most) possible polygons (map masks) we give it. poly2tri is not very good at handeling weird polygons: it does not support repeated vertex points within epsilon, polygons can not be self-intersecting and Steiner points can not be on the existing lines ('collinear') (!). One way to deal with this could be to move the vertex points around by a very small epsilon value in x and y directions. Another way could be to use something like [robust-point-in-polygon](https://github.com/mikolalysenko/robust-point-in-polygon) (5KB instead of the <1KB [point-in-polygon](https://www.npmjs.com/package/point-in-polygon). - Support holes. We don't need them for the Allmaps project, but poly2tri supports them so we could easily add this.35:[] 2f:{"title":"A simple polygon triangulation algorithm","cells":[{"id":0,"value":"# A simple polygon triangulation algorithm","pinned":false,"hidden":false,"mode":"md"},{"id":873,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":889,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":879,"value":"$33","pinned":false,"hidden":false,"mode":"md"},{"id":909,"value":"$34","pinned":false,"hidden":false,"mode":"md"},{"id":884,"value":"## Try it out!","pinned":false,"hidden":false,"mode":"md"},{"id":854,"value":"# ✎👇","pinned":false,"hidden":false,"mode":"md"},{"id":683,"value":"viewof polygon = polygonInput({\n title: \"Draw a polygon\",\n description:\n \"Drawing three or more points creates a polygon. It's triangulation is then shown below.\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1012,"value":"polygon","pinned":true,"hidden":false,"mode":"ojs"},{"id":13,"value":"{\n let width = 400;\n let height = Math.round((210 / 400) * width);\n\n let ctx = DOM.context2d(width, height);\n\n function render() {\n ctx.clearRect(0, 0, width, width);\n if (useGrid) {\n drawGrid(ctx, grid);\n }\n if (polygon.length > 0) {\n drawPolygon(ctx, polygon);\n }\n if (polygon.length > 2) {\n drawTriangles(ctx, triangles);\n }\n }\n\n render();\n\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":800,"value":"Play around with these settings to better understand what's happening","pinned":false,"hidden":false,"mode":"md"},{"id":803,"value":"viewof useGrid = Inputs.toggle({ label: \"Use Grid\", value: true })","pinned":false,"hidden":false,"mode":"ojs"},{"id":805,"value":"viewof usePolygonInterpolation = Inputs.toggle({\n label: \"Use Polygon Interpolation\",\n value: true\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"viewof gridSize = Inputs.range([5, 80], {\n value: 20,\n step: 1,\n label: \"Grid size\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":837,"value":"The triangulation library 'poly2tri.js' used in this code does not allow grid points ('Steiner points') on edges ('collinear'). This may be circumvented my moving the grid points by a small epsilon, or by using the more robust `robust-point-in-polygon` checker.","pinned":false,"hidden":false,"mode":"md"},{"id":834,"value":"viewof useOffset = Inputs.toggle({\n label: \"Use Offset\",\n value: false\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":953,"value":"viewof useRobustPointInPolygon = Inputs.toggle({\n label: \"Use robust-point-in-polygon\",\n value: true\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":818,"value":"---","pinned":false,"hidden":false,"mode":"md"},{"id":666,"value":"### Polygon","pinned":false,"hidden":false,"mode":"md"},{"id":679,"value":"import { polygonInput } from \"@harrystevens/polygon-input\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":61,"value":"drawPolygon = function (ctx, polygon, f) {\n ctx.fillStyle = \"white\";\n for (let i = 0; i + 1 < polygon.length; i++) {\n ctx.beginPath();\n ctx.arc(polygon[i][0], polygon[i][1], 4, 0, 2 * Math.PI, false);\n ctx.moveTo(polygon[i][0], polygon[i][1]);\n ctx.lineTo(polygon[i + 1][0], polygon[i + 1][1]);\n ctx.stroke(); // Render the path\n ctx.fill();\n }\n ctx.beginPath();\n var i = polygon.length - 1;\n ctx.arc(polygon[i][0], polygon[i][1], 4, 0, 2 * Math.PI, false);\n ctx.moveTo(polygon[i][0], polygon[i][1]);\n ctx.lineTo(polygon[0][0], polygon[0][1]);\n ctx.stroke(); // Render the path\n ctx.fill();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":671,"value":"### Grid","pinned":false,"hidden":false,"mode":"md"},{"id":773,"value":"getBbox = function (polygon) {\n // polygon = [[3,5], [7,10], [20, 1]]\n const result = [Infinity, Infinity, -Infinity, -Infinity];\n polygon.forEach((point) => {\n if (result[0] > point[0]) {\n result[0] = point[0];\n }\n if (result[1] > point[1]) {\n result[1] = point[1];\n }\n if (result[2] < point[0]) {\n result[2] = point[0];\n }\n if (result[3] < point[1]) {\n result[3] = point[1];\n }\n });\n return result;\n // result = [minX, minY, maxX, maxY]\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":9,"value":"grid = {\n let grid = [];\n let bbox = getBbox(polygon);\n let epsilon = useOffset ? Math.pow(10, -10) : 0;\n for (let x = bbox[0] + gridSize, i = 0; x <= bbox[2]; i++, x += gridSize) {\n for (let y = bbox[1] + gridSize, j = 0; y <= bbox[3]; j++, y += gridSize) {\n grid.push([x + epsilon, y + epsilon]);\n }\n }\n return grid;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"drawGrid = function (ctx, grid, f) {\n for (let i = 0; i < grid.length; i++) {\n ctx.fillStyle = \"grey\";\n ctx.beginPath();\n ctx.arc(grid[i][0], grid[i][1], 2, 0, 2 * Math.PI, false);\n ctx.fill();\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":584,"value":"### Triangulation imports","pinned":false,"hidden":false,"mode":"md"},{"id":83,"value":"poly2tri = import('https://cdn.skypack.dev/poly2tri@1.5.0?min')","pinned":false,"hidden":false,"mode":"ojs"},{"id":277,"value":"pointInPolygon = require('https://bundle.run/point-in-polygon@1.1.0')","pinned":false,"hidden":false,"mode":"ojs"},{"id":932,"value":"robustPointInPolygon = import(\n \"https://cdn.skypack.dev/robust-point-in-polygon@1.0.3\"\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":576,"value":"### Triangulation helper functions","pinned":false,"hidden":false,"mode":"md"},{"id":324,"value":"getSegmentDistance = function (segment) {\n // segment = [[x0, y0], [x1, y1]]\n // return lenght of segment\n return Math.sqrt(\n Math.pow(segment[1][0] - segment[0][0], 2) +\n Math.pow(segment[1][1] - segment[0][1], 2)\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":339,"value":"getSegmentAngle = function (segment) {\n // segment = [[x0, y0], [x1, y1]]\n // return angle of segment (in radians, signed)\n return Math.atan2(\n segment[1][1] - segment[0][1],\n segment[1][0] - segment[0][0]\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":329,"value":"getNextPoint = function (point, distance, angle) {\n // point = [x0, y0], distance = number, angle = number in radians, signed\n // compute next point from point along angle over distance\n return [\n point[0] + Math.cos(angle) * distance,\n point[1] + Math.sin(angle) * distance\n ];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":298,"value":"interpolateSegment = function (segment, distance) {\n // segment = [[3,5], [7,10]]\n // distance = 1\n\n var currentPoint = segment[0];\n var result = [currentPoint];\n\n while (getSegmentDistance([currentPoint, segment[1]]) > distance) {\n var nextPoint = getNextPoint(\n currentPoint,\n distance,\n getSegmentAngle(segment)\n );\n result.push(nextPoint);\n currentPoint = nextPoint;\n }\n // segment[1] is not pushed\n return result;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":516,"value":"interpolatePolygon = function (polygon, distance) {\n // polygon = [[3,5], [7,10], [20, 1]]\n // distance = 1\n\n // close polygon\n polygon = [...polygon, polygon[0]];\n\n var result = [];\n for (let i = 0; i < polygon.length - 1; i++) {\n result = result.concat(\n interpolateSegment([polygon[i], polygon[i + 1]], distance)\n );\n }\n\n return result;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":607,"value":"### Triangulation","pinned":false,"hidden":false,"mode":"md"},{"id":488,"value":"TRIANG = {\n function TRIANG(polygon, points, distance) {\n // initialise Constrained Delaunay Triangulation with polygon\n let polygonUsed = usePolygonInterpolation\n ? interpolatePolygon(polygon, distance)\n : polygon;\n var swctx = new poly2tri.SweepContext(\n polygonUsed.map((p) => new poly2tri.Point(p[0], p[1]))\n );\n if (useGrid) {\n // add grid points as Steiner points if they are inside the polygon\n for (var i = 0; i < points.length; i++) {\n if (useRobustPointInPolygon) {\n if (robustPointInPolygon.default(polygon, points[i]) == -1) {\n swctx.addPoint(new poly2tri.Point(points[i][0], points[i][1]));\n }\n } else {\n if (pointInPolygon(points[i], polygon)) {\n swctx.addPoint(new poly2tri.Point(points[i][0], points[i][1]));\n }\n }\n }\n }\n // triangulate\n swctx.triangulate();\n let triangles = swctx.getTriangles();\n return triangles;\n }\n return TRIANG;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":478,"value":"triangles = polygon.length > 2 ? TRIANG(polygon, grid, gridSize) : []","pinned":false,"hidden":false,"mode":"ojs"},{"id":988,"value":"To get the result as a list of polygons:","pinned":false,"hidden":false,"mode":"md"},{"id":964,"value":"triangles.map((t) => t.getPoints().map((p) => [p.x, p.y]))","pinned":true,"hidden":false,"mode":"ojs"},{"id":115,"value":"drawTriangles = function (ctx, triangles) {\n triangles.forEach(function (t) {\n ctx.fillStyle = \"hsl(\" + 360 * Math.random() + \", 50%, 50%, 50%)\";\n ctx.beginPath();\n ctx.moveTo(t.getPoint(0).x, t.getPoint(0).y);\n ctx.lineTo(t.getPoint(1).x, t.getPoint(1).y);\n ctx.lineTo(t.getPoint(2).x, t.getPoint(2).y);\n ctx.fill();\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1024,"value":"test = function () {\n const pixelTriangleIndex = [];\n for (let index = 0; index < triangles.length; index++) {\n pixelTriangleIndex[3 * index] = index;\n pixelTriangleIndex[3 * index + 1] = index;\n pixelTriangleIndex[3 * index + 2] = index;\n }\n return pixelTriangleIndex;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1029,"value":"test()","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}