CCPD Snowden

Using a Capacity Constrained Point Distribution to display a picture of an American hero.

This is a variant of the original algorithm by Michael Balzer, Thomas Schlömer & Oliver Deussen (University of Konstanz, Germany, 2009), in which I use a Voronoi Diagram to create a topology of the current sites, and only swap the points between neighbouring sites (and neighbours of neighbours). It appears to be much faster than the original algorithm.

See also https://seenthis.net/messages/619204 for an application in cartography.

Note to self: do not edit on BlockBuilder

<!DOCTYPE html>
<meta charset="utf-8">
<style>
    .polygons {
        fill: none;
        stroke: #333;
    }
    
    canvas {
        display: none;
    }
</style>
<script src="https://d3js.org/d3.v4.min.js"></script>

function draw(img) {

var canvas = d3.select('body')
            .append('canvas')
            .attr('width', img.width)
            .attr('height', img.height)
            .node();

var ctx = canvas.getContext('2d');
        ctx.drawImage(img, 0, 0);
        var imageData = ctx.getImageData(0, 0, canvas.width, canvas.height);
        var data = imageData.data;

var height = 500,
            width = 960;

var svg = d3.select("body")
            .append('svg')
            .attr('width', width)
            .attr('height', height);

var points = [];
        var factor = height / canvas.height;
        for (var i = 0; i < data.length; i += 4) {
            if (Math.random() > (data[i] + data[i + 1] + data[i + 2]) / (3 * 256)) {
                var y = Math.floor(i / 4 / canvas.width),
                    x = (i / 4) - y * canvas.width;
                points.push([x * factor, y * factor]);
            }
        }

const capacity = 14;
        const npoints = Math.floor(points.length / capacity);

var sites = d3.range(npoints)
            .map(function (d) {
                var p = points[Math.floor(Math.random() * points.length)];
                p.index = d;
                return p;
            });

var voronoi = d3.voronoi().size([width, height]),
            diagram = voronoi(sites);

var polygon = svg.append("g")
            .attr("class", "polygons")
            .selectAll("path")
            .data(diagram.polygons())
            .enter().append("path")
            .attr('stroke-opacity', 0.1);

function distance2(a, b) {
            var dx = a[0] - b[0],
                dy = a[1] - b[1];
            return /*Math.sqrt*/ (dx * dx + dy * dy);
        }

var calc = 0,
            iterations = 0;

function iterate() {
            var swaps = 0;
            iterations++;

var links = new Array(sites.length);
            diagram.links().forEach(function (l) {
                var ext = d3.extent([l.source.index, l.target.index]),
                    i = ext[0],
                    j = ext[1];
                if (!links[i]) links[i] = [j];
                else links[i].push(j);
            });

for (var i in links) {
                var l = links[i];
                /* if (false)  */
                links[i] = d3.merge([links[i]].concat(links[i].map(function (j) {
                    return links[j] || [];
                })));

l.forEach(function (j) {
                    var Hi = [],
                        Hj = [],
                        k, ki, kj;
                    for (k = 0; k < capacity; k++) {
                        Hi.push(distance2(points[i * capacity + k], sites[j]) - distance2(points[i * capacity + k], sites[i]));
                        Hj.push(distance2(points[j * capacity + k], sites[i]) - distance2(points[j * capacity + k], sites[j]));

calc++;
                    }

while (Hi.length > 0 && Hj.length > 0 && ((ki = d3.scan(Hi)) || true) && ((kj = d3.scan(Hj)) || true) && (Hi[ki] + Hj[kj] < 0)) {
                        swaps++;
                        var temp = points[i * capacity + ki];
                        points[i * capacity + ki] = points[j * capacity + kj];
                        points[j * capacity + kj] = temp;
                        Hi = Hi.slice(0, ki).concat(Hi.slice(ki + 1));
                        Hj = Hj.slice(0, kj).concat(Hj.slice(kj + 1));
                    }
                });
            }

return swaps;
        }

var site = svg.selectAll('circle')
            .data(sites)
            .enter()
            .append('circle')
            .attr('r', 1.5)
            .attr('fill', '#333');

var lastswaps = null;
        var interval = d3.interval(function () {
            var swaps = iterate();

svg.selectAll('circle')
                .data(sites)
                .attr('transform', function (d) {
                    return 'translate(' + d + ')';
                });

svg.selectAll('rect')
                .data(points)
                .attr('transform', function (d) {
                    return 'translate(' + d + ')';
                });

diagram = voronoi(sites);

polygon = polygon.data(diagram.polygons())
                .attr("d", function (d) {
                    return d ? "M" + d.join("L") + "Z" : null;
                });

sites = sites.map(function (site, i) {
                var pts = points.slice(i * capacity, i * capacity + capacity);
                site[0] = d3.mean(pts.map(function (d) {
                    return d[0];
                }));
                site[1] = d3.mean(pts.map(function (d) {
                    return d[1];
                }));
                return site;
            });

console.log("" + swaps + " swaps, " + calc + " distances computed");
            if (swaps == lastswaps && swaps < 300) {
                console.log("stabilized after " + iterations + " iterations.")
                interval.stop();
                polygon
                    .transition()
                    .duration(2000)
                    .attr('stroke-opacity', 0.05);
                site
                    .transition()
                    .duration(2000)
                    .attr('fill', 'black');
            }
            lastswaps = swaps;
        })

}
</script>

https://d3js.org/d3.v4.min.js