|
<!DOCTYPE html> |
|
<h1>Original Points</h1> |
|
<svg id="originalpoints" width="960" height="673" stroke="#fff" stroke-width="0.5"></svg> |
|
<h1>Step 1: delaunay triagulation: find the smallest triangle of the closest 3 points</h1> |
|
<canvas id="s3" width="960" height="673"></canvas> |
|
<h1>Step 2: pick a regular grid and interpolate the grid point</h1> |
|
<p>grid point is interpolated from the triangle that contains the point using "Berycentric Interpolation" algorithm</p> |
|
<svg id="grid" width="960" height="673" stroke="#fff" stroke-width="0.5"></svg> |
|
<h1>Step 3: Contour plot using "marching squares" algorithm</h1> |
|
<svg id="contourplot" width="960" height="673" stroke="#fff" stroke-width="0.5"></svg> |
|
<script src="https://d3js.org/d3.v4.min.js"></script> |
|
<script src="https://d3js.org/d3-hsv.v0.1.min.js"></script> |
|
<script src="https://d3js.org/d3-contour.v1.min.js"></script> |
|
<script src="https://unpkg.com/[email protected]/delaunay.js"></script> |
|
<script src="https://unpkg.com/lodash/lodash.min.js"></script> |
|
<script> |
|
// grided points |
|
d3.json("ca55.json", function (error, aero) { |
|
if (error) throw error; |
|
var i0 = d3.interpolateHsvLong(d3.hsv(120, 1, 0.65), d3.hsv(60, 1, 0.90)), |
|
i1 = d3.interpolateHsvLong(d3.hsv(60, 1, 0.90), d3.hsv(0, 0, 0.95)), |
|
interpolateTerrain = function (t) { return t < 0.5 ? i0(t * 2) : i1((t - 0.5) * 2); }; |
|
let triangles; |
|
|
|
const xmin = _.minBy(aero, getx)[0]; |
|
const ymin = _.minBy(aero, gety)[1]; |
|
const zmin = _.minBy(aero, getz)[2]; |
|
const xmax = _.maxBy(aero, getx)[0]; |
|
const ymax = _.maxBy(aero, gety)[1]; |
|
const zmax = _.maxBy(aero, getz)[2]; |
|
|
|
// contour |
|
(function () { |
|
const svg = d3.select("#contourplot"), |
|
width = +svg.attr("width"), |
|
height = +svg.attr("height"); |
|
|
|
const x = d3.scaleLinear().range([0, width]).domain([xmin, xmax]), |
|
y = d3.scaleLinear().range([0, height]).domain([ymin, ymax]), |
|
color = d3.scaleSequential(interpolateTerrain).domain([zmin, zmax]); |
|
|
|
const gridCellSize = 15 // unit: pixel |
|
const levels = 10; // number of contour levels |
|
|
|
const ib = Math.trunc(width / gridCellSize); |
|
const jb = Math.trunc(height / gridCellSize); |
|
|
|
triangles = new Delaunay(aero).triangulate(); |
|
|
|
console.time("generateData"); |
|
// generate a regular grid |
|
const data = new Float64Array(ib * jb); |
|
const gridData = []; |
|
for (let i = 0; i < ib; i++) { |
|
for (let j = 0; j < jb; j++) { |
|
const p = [x.invert((i + 0.5) * gridCellSize), y.invert((j + 0.5) * gridCellSize)]; |
|
|
|
let value; |
|
for (let i = 0; i < triangles.length; i += 3) { |
|
const param = contains(triangles[i], triangles[i + 1], triangles[i + 2], p); |
|
if (param) { |
|
const [u, v] = param |
|
value = (1 - u - v) * triangles[i][2] + u * triangles[i][2] + v * triangles[i][2] |
|
break; |
|
} |
|
} |
|
if (value == null) { |
|
// FIXME: should try to find the closest point? |
|
value = zmin |
|
} |
|
data[j * ib + i] = value; |
|
gridData[j * ib + i] = p.concat(value); |
|
} |
|
} |
|
console.timeEnd("generateData"); |
|
console.log("data size", data.length); |
|
|
|
svg.selectAll("path") |
|
.data(d3.contours() |
|
.size([ib, jb]) |
|
.thresholds(d3.range(zmin, zmax, (zmax - zmin) / levels)) |
|
(data)) |
|
.enter().append("path") |
|
.attr("d", d3.geoPath(d3.geoIdentity().scale(width / ib))) |
|
.attr("fill", function (d) { return color(d.value); }); |
|
const svgGridPoints = d3.select("#grid"); |
|
svgGridPoints.selectAll("circle") |
|
.data(gridData) |
|
.enter().append("circle") |
|
.attr('cx', function (d) { return x(d[0]) }) |
|
.attr('cy', function (d) { return y(d[1]) }) |
|
.attr('r', 2) |
|
.attr('fill', function (d) { return color(d[2]) }); |
|
}()); |
|
// original points |
|
(function () { |
|
let svg = d3.select("#originalpoints"), |
|
width = +svg.attr("width"), |
|
height = +svg.attr("height"); |
|
|
|
const xmin = _.minBy(aero, getx)[0]; |
|
const ymin = _.minBy(aero, gety)[1]; |
|
const zmin = _.minBy(aero, getz)[2]; |
|
const xmax = _.maxBy(aero, getx)[0]; |
|
const ymax = _.maxBy(aero, gety)[1]; |
|
const zmax = _.maxBy(aero, getz)[2]; |
|
|
|
let x = d3.scaleLinear().range([0, width]).domain([xmin, xmax]), |
|
y = d3.scaleLinear().range([0, height]).domain([ymin, ymax]), |
|
color = d3.scaleSequential(interpolateTerrain).domain([zmin, zmax]); |
|
|
|
svg.selectAll("circle") |
|
.data(aero) |
|
.enter().append("circle") |
|
.attr('cx', function (d) { return x(d[0]) }) |
|
.attr('cy', function (d) { return y(d[1]) }) |
|
.attr('r', 2) |
|
.attr('fill', function (d) { return color(d[2]) }); |
|
}()); |
|
// triagulation |
|
(function () { |
|
var canvas = document.getElementById("s3"), |
|
ctx = canvas.getContext("2d"); |
|
const x = d3.scaleLinear().range([0, 960]).domain([xmin, xmax]), |
|
y = d3.scaleLinear().range([0, 673]).domain([ymin, ymax]); |
|
for (i = triangles.length; i;) { |
|
ctx.beginPath(); |
|
--i; ctx.moveTo(x(triangles[i][0]), y(triangles[i][1])); |
|
--i; ctx.lineTo(x(triangles[i][0]), y(triangles[i][1])); |
|
--i; ctx.lineTo(x(triangles[i][0]), y(triangles[i][1])); |
|
ctx.closePath(); |
|
ctx.stroke(); |
|
} |
|
}()) |
|
|
|
}) |
|
|
|
// https://github.com/ironwallaby/delaunay/blob/2c60fbd36662ae970b1e8402b18285317a60d0b6/delaunay.js#L203 |
|
const contains = function (p0, p1, p2, p) { |
|
/* Bounding box test first, for quick rejections. */ |
|
if ((p[0] < p0[0] && p[0] < p1[0] && p[0] < p2[0]) || |
|
(p[0] > p0[0] && p[0] > p1[0] && p[0] > p2[0]) || |
|
(p[1] < p0[1] && p[1] < p1[1] && p[1] < p2[1]) || |
|
(p[1] > p0[1] && p[1] > p1[1] && p[1] > p2[1])) |
|
return null; |
|
|
|
var a = p1[0] - p0[0], |
|
b = p2[0] - p0[0], |
|
c = p1[1] - p0[1], |
|
d = p2[1] - p0[1], |
|
i = a * d - b * c; |
|
|
|
/* Degenerate tri. */ |
|
if (i === 0.0) |
|
return null; |
|
|
|
var u = (d * (p[0] - p0[0]) - b * (p[1] - p0[1])) / i, |
|
v = (a * (p[1] - p0[1]) - c * (p[0] - p0[0])) / i; |
|
|
|
/* If we're outside the tri, fail. */ |
|
if (u <= 0.0 || v <= 0.0 || (u + v) >= 1.0) |
|
return null; |
|
|
|
return [u, v]; |
|
} |
|
|
|
// berycentric interpolate |
|
// http://answers.unity3d.com/questions/383804/calculate-uv-coordinates-of-3d-point-on-plane-of-m.html |
|
// https://math.stackexchange.com/questions/516219/finding-out-the-area-of-a-triangle-if-the-coordinates-of-the-three-vertices-are |
|
// p: [x,y,z?], tri: [p0,p1,p2] |
|
// turns out method "contains" already calculates this |
|
const interpolate = function (tri, p) { |
|
const A = (tri[1][0] - tri[0][0]) * (tri[2][1] - tri[0][1]) - (tri[2][0] - tri[0][0]) * (tri[1][1] - tri[0][1]) |
|
const A1 = (tri[1][0] - p[0]) * (tri[2][1] - p[1]) - (tri[2][0] - p[0]) * (tri[1][1] - p[1]) |
|
const A2 = (p[0] - tri[0][0]) * (tri[2][1] - tri[0][1]) - (tri[2][0] - tri[0][0]) * (p[1] - tri[0][1]) |
|
const A3 = (tri[1][0] - tri[0][0]) * (p[1] - tri[0][1]) - (p[0] - tri[0][0]) * (tri[1][1] - tri[0][1]) |
|
|
|
const u = A2 / A, v = A3 / A |
|
|
|
return (1 - u - v) * tri[0][2] + u * tri[1][2] + v * tri[2][2] |
|
} |
|
|
|
const getx = function (o) { return o[0]; } |
|
const gety = function (o) { return o[1]; } |
|
const getz = function (o) { return o[2]; } |
|
|
|
</script> |