Stippling of 2D Scalar Fields
This notebook shows an example of using the Linde–Buzo–Gray stippling algorithm to abstract two-dimensional scalar fields. Please note that this is not the actual implementation that this is a simplified version of the actual algorithm—a more detailed description of our method can be found in the paper.
For demonstration purposes, we will use the Himmelblau function, which is often used to benchmark optimization algortihms.
himmelblau = (x, y) => Math.pow(x * x + y - 11, 2) + Math.pow(x + y * y - 7, 2);
We can put the values in a grid and visualize the function using color.
himmelblau_grid = new Grid([400,400], himmelblau, [-5,5], [-5,5], [0, 800]);
For and it looks like this:
viewof canvas = { return himmelblau_grid.to_context().canvas }
Now, our goal is to create a stippled abstraction of this function. For this we will use an adaption of the Linde–Buzo-Gray stippling algorithm.
We are now ready to iterate the algorithm.
stipples_new = { const lbg_stippling = new Stippling(himmelblau_grid); lbg_stippling.with_adaptive_radius(1,5); let status = lbg_stippling.iterate(); const svg = d3.select(DOM.svg(himmelblau_grid.width, himmelblau_grid.height)); while(status.splitted.length != 0 || status.deleted.length != 0) { console.log(status); status = lbg_stippling.iterate(); const circles = svg.selectAll("circle") .data(lbg_stippling.stipples); circles.enter().append("circle") .attr("r", d => d.radius) .attr("cx", d => d[0]) .attr("cy", d => d[1]); circles.merge(circles) .attr("r", d => d.radius) .attr("cx", d => d[0]) .attr("cy", d => d[1]); circles.exit().remove(); yield Promises.delay(100, svg.node()); } yield svg.node(); }
html`<hr>`
Below, we define a Grid type and the stippling algorithm, which we use to create the examples above.
class Grid { // `func` has to be a 2d function constructor([n, m], func, [xmin, xmax], [ymin, ymax], [vmin, vmax]) { this.width = n; this.height = m; this.x = d3.scaleLinear().domain([0.5, this.width]).range([xmin, xmax]); this.y = d3.scaleLinear().domain([0.5, this.height]).range([ymax, ymin]); // Turn upside down this.color = d3.scaleSequentialSqrt(d3.interpolateYlGnBu).domain([vmin, vmax]); this.values = new Array(this.width * this.height); // populate array with values for (var j = 0.5, k = 0; j < this.height; ++j) { for (var i = 0.5; i < this.width; ++i, ++k) { this.values[k] = func(this.x(i), this.y(j)); } } } get(i, j) { return this.values[i + j * this.width]; } pos_at_index(idx) { const x = idx % this.width; const y = (idx - x) / this.width; return [x + 0.5, y + 0.5]; // We are looking at centers of the grid. } to_context() { let ctx = DOM.context2d(this.width, this.height); ctx.canvas.width = this.width; ctx.canvas.height = this.height; for (let u of d3.range(this.width)) { for (let v of d3.range(this.height)) { ctx.fillStyle = this.color(this.get(u,v)); ctx.fillRect(u,v,1,1); } } ctx.canvas.id = "canvas" return ctx; } }
class Stippling { constructor(grid) { this.grid = grid; this.radius_extent = [3, 3]; this.threshold = 0.4; this.delta_threshold = 0.01; // we initialize the algorithm with random stipples this.stipples = this.random_stipples(100); } // enables adaptive stipple sizes, depending on the underlying density. with_adaptive_radius(min_radius, max_radius) { this.radius_extent = [min_radius, max_radius]; return this; } // returns the thresholds for split and delete thresholds(radius) { const area = Math.PI * radius * radius; return [(1.0 + this.threshold / 2.0) * area, (1.0 - this.threshold / 2.0) * area]; } // generates `num` random stipples on top of the grid. random_stipples(num) { const stipples = new Array(num); const xran = d3.randomUniform(0, this.grid.width); const yran = d3.randomUniform(0, this.grid.height); for (var i = 0; i < num; ++i) { stipples[i] = [xran(), yran()]; stipples[i].radius = this.radius_extent[0]; } return stipples; } // performs one iteration of the stippling algorithm iterate() { const voronoi = d3.voronoi() .extent([[0, 0], [this.grid.width, this.grid.height]]) .x(d => d[0]) .y(d => d[1]); const diagram = voronoi(this.stipples); // initialize densities diagram.cells.forEach(cell => { cell.site.data.density = 0; cell.site.data.moment10 = 0; cell.site.data.moment01 = 0; cell.site.data.moment11 = 0; cell.site.data.moment20 = 0; cell.site.data.moment02 = 0; }); // compute the density and the weighted centroid of each cell this.grid.values.forEach((val, idx) => { const coord = this.grid.pos_at_index(idx); var closest_cell = diagram.find(coord[0], coord[1]); closest_cell.data.density += val; // Moment00 closest_cell.data.moment10 += coord[0] * val; closest_cell.data.moment01 += coord[1] * val; closest_cell.data.moment11 += coord[0] * coord[1] * val; closest_cell.data.moment20 += coord[0] * coord[0] * val; closest_cell.data.moment02 += coord[1] * coord[1] * val; }); const grid_extent = d3.extent(this.grid.values); const density_to_radius = d3.scaleLinear().domain(grid_extent).range(this.radius_extent); let deleted = []; let splitted = []; let relaxed = []; diagram.polygons().forEach(polygon => { const density = polygon.data.density; // determine point size based on average intensity const area = d3.polygonArea(polygon); const avg_density = density / area; const radius = density_to_radius(avg_density); const [split_threshold, delete_threshold] = this.thresholds(radius); if (density < delete_threshold * grid_extent[1]) { // Delete deleted.push[polygon.data]; // nop } else if (density > split_threshold * grid_extent[1]) { // Split const centroid = d3.polygonCentroid(polygon); const dist = Math.sqrt(d3.polygonArea(polygon) / Math.PI) / 2.0; const x = polygon.data.moment20 / polygon.data.density - centroid[0] * centroid[0]; const y = 2 * (polygon.data.moment11 / polygon.data.density - centroid[0] * centroid[1]); const z = polygon.data.moment02 / polygon.data.density - centroid[1] * centroid[1]; var orientation = Math.atan2(y, x - z) / 2.0; var deltaX = dist * Math.cos(orientation); var deltaY = dist * Math.sin(orientation); splitted.push([centroid[0] + deltaX, centroid[1] + deltaY]); splitted.push([centroid[0] - deltaX, centroid[1] - deltaY]); splitted[splitted.length - 2].radius = radius; splitted[splitted.length - 1].radius = radius; } else { // Relax const centroid = [polygon.data.moment10 / polygon.data.density, polygon.data.moment01 / polygon.data.density]; relaxed.push(centroid); relaxed[relaxed.length - 1].radius = radius; } }); // we increase the threshold with each iteration for faster convergence this.threshold = this.threshold + this.delta_threshold; this.stipples = [...relaxed, ...splitted]; return {deleted, relaxed, splitted}; } }
Below we just include d3 for convenience:
d3 = require('d3');