Voronoï playground: Voronoï map (study 2)

This block is a Work In Progress, a continuation of a previous block, and tries to implements a Voronoï map (i.e., one-level treemap) based on the algorithm presented in Computing Voronoi Treemaps - Faster, Simpler, and Resolution-independent, and augmented with heuristics found on a Java implementation done by Arlind Nocaj, one of the co-authors of the scientific paper. Because this algorithm requires the computation of weighted Voronoï, this block uses the d3-weighted-voronoi plugin.

Without going into details (cf. the scientific paper), the overall algorithm doesn't differs from the one used in the previous block. Hence, at each iteration :

the weighted voronoi diagram is computed based on each weighted sites
then, each site is re-position at the center of its influence area
then, another weighted voronoi diagram is computed based on each weighted relocated sites
then, each site is re-weighted in order to fit its expected influence area

The algorithm stops when each site is no longer moved or re-weighted (more exactly, when each site is moved or is re-weighted below a certain treshold), or if a certain number of iterations is reached.

Differences from the algorithm used in the previous block are :

rework the limitation of weights when re-positioning and when re-weighting sites; both handled in function handleOverweigthed, which prevent any site to be overweighted (and hence prevent any site to produce an empty cell)
heuristics to mitigate near-zero weight (if active, min allowed weight = (max weight / 1000)); near-zero weights are difficult to handle because re-positioning systematically overweighted near-zero sites, making stabilization difficult to reach
heuristics to mitigate flickering (if active, if area error < 5% of total area, the more the flickering count, the less sites are re-positioned and re-weigted)

Notes :

only one level, no nested hierarchy
compared to the previous block, the algorithm used in this block performs better : in term of representativity of initial weights (reaching error below 1% is common, compared to 15%-20% in the previous block), and in term of number of iterations needed to reach the stable state
experience shows that stabilization is not reached when a lightweighted site is trapped between 2 heavyweighted sites and the edges of the encosing polygon; this is because (i) heavyweighted sites have reached their targeted area and don't need to be re-positioned or re-weighted by their own, (ii) lightweighted site has much larger area (than the one targeted) but can no longer lowers its weight (it has reached the minimim weight), and (iii) the others sites haven't enought energy to grow and reach their targeted area (which would push heavyweighedt sites on the lightweighted one)

User interactions :

you can choose to draw the Weighted Voronoï Diagram (default) or the weights (visualized as circles).
you can hide/show sites
you can choose among different rendering (greyscale, radial rainbow, or conical rainbow (default, having hard-&-fun time to implement it because canvas don't provides conical gradient)).

Acknowledgments to :

D3.js (v.4)
d3-weighted-voronoi plugin
blockbuilder.org
(part of) https://github.com/ArlindNocaj/Voronoi-Treemap-Library for a Java implementation

<html lang="en">
  <head>
    <meta charset="utf-8" />
    <title>Voronoï playground: Voronoï map (study 2)</title>
    <meta name="description" content="Voronoi map in D3.js (study 2)">
    <script src="https://d3js.org/d3.v4.min.js"></script>
    <script src="https://raw.githack.com/kcnarf/d3-weighted-voronoi/master/build/d3-weighted-voronoi.js"></script>
    <style>
      #layouter {
        text-align: center;
        position: relative;
      }
      
      #wip {
        display: none;
        position: absolute;
        top: 200px;
        left: 330px;
        font-size: 40px;
        text-align: center;
      }

.control {
        position: absolute;
      }
      .control.top{
        top: 5px;
      }
      .control.bottom {
        bottom: 5px;
      }
      .control.left{
        left: 5px;
      }
      .control.right {
        right: 5px;
      }
      .control.right div{
        text-align: right;
      }
      .control.left div{
        text-align: left;
      }
      .control .separator {
        height: 5px;
      }

canvas {
        margin: 1px;
        border-radius: 1000px;
        box-shadow: 2px 2px 6px grey;
      }
      canvas#background-image, canvas#alpha {
        display: none;
      }
    </style>
  </head>
  <body>
    <div id="layouter">
      <canvas id="background-image"></canvas>
      <canvas id="alpha"></canvas>
      <canvas id="colored"></canvas>

<div id="control0" class="control top left">
        <div>
          <input id="cellsOrWeights" type="radio" name="cellsOrWeights" value="cells" checked onchange="cellsOrWeightsUpdated('cells')"> Weighted Voronoï 
        </div>
        <div>
          <input id="cellsOrWeights" type="radio" name="cellsOrWeights" value="circles" onchange="cellsOrWeightsUpdated('weights')"> Weights
        </div>
        <div>
          <input id="cellsOrWeights" type="radio" name="cellsOrWeights" value="circles" onchange="cellsOrWeightsUpdated('cellsAndWeights')"> Both
        </div>
      </div>
      
      <div id="control1" class="control bottom left">
        <div>
          <input id="showSites" type="checkbox" name="showSites" onchange="siteVisibilityUpdated()"> Show sites
        </div>
      </div>
      
      <div id="control2" class="control bottom right">
        <div>
          Grey <input id="bgImgGrey" type="radio" name="bgImg" onchange="bgImgUpdated('grey')">
        </div>
        <div>
          Radial rainbow <input id="bgImgRadialRainbow" type="radio" name="bgImg" onchange="bgImgUpdated('radialRainbow')">
        </div>
        <div>
          Canonical rainbow <input id="bgImgCanonicalRainbow" type="radio" name="bgImg" checked onchange="bgImgUpdated('canonicalRainbow')">
        </div>
      </div>

<div id="wip">
        Work in progress ...
      </div>
    </div>
    
    <script>
      var _PI = Math.PI,
          _2PI = 2*Math.PI,
          sqrt = Math.sqrt,
          sqr = function(d) { return Math.pow(d,2); }
      		epsilon = 1;
      
      //begin: layout conf.
      var totalWidth = 550,
          totalHeight = 500,
          controlsHeight = 0,
          canvasbw = 1, // canvas border width
          canvasbs = 8, // canvas box-shadow
          radius = (totalHeight-controlsHeight-canvasbs-2*canvasbw)/2,
          width = 2*radius,
          height = 2*radius,
          halfRadius = radius/2
          halfWidth = halfRadius,
          halfHeight = halfRadius,
          quarterRadius = radius/4;
          quarterWidth = quarterRadius,
          quarterHeight = quarterRadius;
      //end: layout conf.
      
      //begin: drawing conf.
      var drawSites = false,
          bgType = "canonicalRainbow",
      		drawCellsOrWeights = "cells",
          bgImgCanvas, alphaCanvas, coloredCanvas,
          bgImgContext, alphaContext, coloredContext,
      		radialGradient;
      //end: drawing conf.
      
      //begin: init layout
      initLayout()
      //end: init layout
      
      //begin: voronoi treemap conf.
      var clippingPolygon = computeClippingPolygon(),
      		weightedVoronoi = d3.weightedVoronoi().clip(clippingPolygon),
          siteCount = 20,
          baseWeight = 2000,
          outlierCount = siteCount/10,
          outlierWeight = 3*baseWeight,
          convergenceTreshold = 0.01, // 0.01 means 1% error
      		totalArea = Math.abs(d3.polygonArea(clippingPolygon)),
      		areaErrorTreshold = convergenceTreshold*totalArea,
          areaErrorHistory = [], // used to detect flickering
          areaErrorHistoryLength = 10,
          maxIterationCount = 200,
          shouldBreakOnMaxIteration = true;
      var totalWeight, avgWeight;
      //end: voronoi treemap conf.
      
      //begin: algorithm conf.
      var shouldComputeVoronoiAfterReposition = true,
      		handleOverweightedVariant = 1,
          shouldRotateHandleOverweighted = false, // mixing severall heuristics seems to performs better (limit flickering), but sometimes freezes (infinite loop in handleOverweighted1)
          shouldMinimizeWeight = false, // when activated, not flickering, but stabilization at higher iterations
          shouldHandleNearZeroWeights = true,
          nearZeroWeightRatio = 0.01, // 0.01 means min allowed weight = 1% of max weight
          adaptPlacementsVariant = 1,
          adaptWeightsVariant = 1;
      var handleOverweight;
      //end: algorithm conf.
			
      function computeClippingPolygon() {
        var circlingPolygon = [];
        for (a=0; a<_2PI; a+=_2PI/60) {
          circlingPolygon.push(
            [radius + (radius-1)*Math.cos(a), radius + (radius-1)*Math.sin(a)]
          )
        }
        
      	return circlingPolygon;
      };
      
      //begin: user interaction handlers
      function siteVisibilityUpdated() {
        drawSites = d3.select("#showSites").node().checked;
      };
      
      function bgImgUpdated(newType) {
        bgType = newType;
        resetBackgroundImage() ;
      };
      
      function cellsOrWeightsUpdated(type) {
        drawCellsOrWeights = type;
      };
      //end: user interaction handlers
      
      function adapt(polygons, iterationCount) {
        var converged, adaptedTreemapPoints;
        
        if (shouldRotateHandleOverweighted) {
          rotateHandleOverweighted();
        }
        
        adaptPlacements(polygons);
        if (shouldComputeVoronoiAfterReposition) {
          adaptedTreemapPoints = polygons.map(function(p) { return p.site.originalObject; });
          polygons = weightedVoronoi(adaptedTreemapPoints);
          if (polygons.length<siteCount) {
            console.log("at least 1 site has no area, which is not supposed to arise");
            debugger;
          }
        }
        
        adaptWeights(polygons);
        adaptedTreemapPoints = polygons.map(function(p) { return p.site.originalObject; });
        polygons = weightedVoronoi(adaptedTreemapPoints);
        if (polygons.length<siteCount) {
          console.log("at least 1 site has no area, which is not supposed to arise");
          debugger;
        }
        
        redraw(adaptedTreemapPoints, polygons);
        
        converged = overallConvergence(polygons);
				
        if (shouldBreakOnMaxIteration && iterationCount===maxIterationCount) {
          console.log("Max iteration reached")
          setTimeout(reset, 1750);
        } else {
          if (converged) {
            console.log("Stopped at iteration "+iterationCount);
            finalize(adaptedTreemapPoints, polygons, 20);
          } else {
            setTimeout(function(){
              adapt(polygons, iterationCount+1);
            }, 50);
          }
        }
      };
      
			function finalize(treemapPoints, polygons, countDown) {
        redraw(treemapPoints, polygons);
        
        if (countDown === 0) {
          setTimeout(reset, 1750);
        } else {
          setTimeout(function(){
            finalize(treemapPoints, polygons, countDown-1);
          }, 50);
        }
      }
      
      function adaptPlacements0(polygons) {
        var newTreemapPoints = [];
        var polygon, treemapPoint, centroid;
        
        for(var i=0; i<siteCount; i++) {
          polygon = polygons[i];
          treemapPoint = polygon.site.originalObject;
          centroid = d3.polygonCentroid(polygon);
          
          treemapPoint.x = centroid[0];
          treemapPoint.y = centroid[1];
          
          newTreemapPoints.push(treemapPoint);
        }
        
        handleOverweighted(newTreemapPoints);
      };
      
      // flickering mitigation
      function adaptPlacements1(polygons) {
        var newTreemapPoints = [];
        var polygon, treemapPoint, centroid, flickeringInfluence;
        
        flickeringInfluence = 0.5*flickeringMitigationRatio(polygons);
        for(var i=0; i<siteCount; i++) {
          polygon = polygons[i];
          treemapPoint = polygon.site.originalObject;
          centroid = d3.polygonCentroid(polygon);
          
          dx = centroid[0] - treemapPoint.x;
          dy = centroid[1] - treemapPoint.y;
          
          //begin: handle excessive change;
          dx *= (1-flickeringInfluence);
          dy *= (1-flickeringInfluence);
          //end: handle excessive change;
          
          
          treemapPoint.x += dx;
          treemapPoint.y += dy;
          
          newTreemapPoints.push(treemapPoint);
        }
        
        handleOverweighted(newTreemapPoints);
      };
      
      function adaptWeights0(polygons) {
        var newTreemapPoints = [];
        var polygon, treemapPoint, currentArea, adaptRatio;
        
        for(var i=0; i<siteCount; i++) {
          polygon = polygons[i];
          treemapPoint = polygon.site.originalObject;
          currentArea = d3.polygonArea(polygon);
          adaptRatio = treemapPoint.targetedArea/currentArea;
          
          //begin: handle excessive change;
          adaptRatio = Math.max(adaptRatio, 0.9);
          adaptRatio = Math.min(adaptRatio, 1.1);
          //end: handle excessive change;
          
          adaptedWeight = treemapPoint.weight*adaptRatio;
          adaptedWeight = Math.max(adaptedWeight, epsilon);
					
          treemapPoint.weight = adaptedWeight;
          
          newTreemapPoints.push(treemapPoint);
        }
        
        handleOverweighted(newTreemapPoints);
      };
      
      // flickering mitigation
      function adaptWeights1(polygons) {
        var newTreemapPoints = [];
        var polygon, treemapPoint, currentArea, adaptRatio, flickeringInfluence;
        
        flickeringInfluence = 0.1*flickeringMitigationRatio(polygons);
        for(var i=0; i<siteCount; i++) {
          polygon = polygons[i];
          treemapPoint = polygon.site.originalObject;
          currentArea = d3.polygonArea(polygon);
          adaptRatio = treemapPoint.targetedArea/currentArea;
          
          //begin: handle excessive change;
          adaptRatio = Math.max(adaptRatio, 0.9+flickeringInfluence);
          adaptRatio = Math.min(adaptRatio, 1.1-flickeringInfluence);
          //end: handle excessive change;
          
          adaptedWeight = treemapPoint.weight*adaptRatio;
          adaptedWeight = Math.max(adaptedWeight, epsilon);
					
          treemapPoint.weight = adaptedWeight;
          
          newTreemapPoints.push(treemapPoint);
        }
        
        handleOverweighted(newTreemapPoints);
      };
      
      // heuristics: lower heavy weights
      function handleOverweighted0(treemapPoints) {
      	var fixCount = 0;
        var fixApplied, tpi, tpj, weightest, lightest, sqrD, adaptedWeight;
        do {
          fixApplied = false;
          for(var i=0; i<siteCount; i++) {
            tpi = treemapPoints[i];
            for(var j=i+1; j<siteCount; j++) {
              tpj = treemapPoints[j];
              if (tpi.weight > tpj.weight) {
                weightest = tpi;
                lightest = tpj;
              } else {
                weightest = tpj;
                lightest = tpi;
              }
              sqrD = squaredDistance(tpi, tpj);
              if (sqrD < weightest.weight-lightest.weight) {
                // adaptedWeight = sqrD - epsilon; // as in ArlindNocaj/Voronoi-Treemap-Library
                // adaptedWeight = sqrD + lightest.weight - epsilon; // works, but below loc performs better (less flickering)
                adaptedWeight = sqrD + lightest.weight/2;
                adaptedWeight = Math.max(adaptedWeight, epsilon);
                weightest.weight = adaptedWeight;
                fixApplied = true;
                fixCount++;
                break;
              }
            }
            if (fixApplied) { break; }
          }
        } while (fixApplied)
        
        if (fixCount>0) {
          if (shouldMinimizeWeight) {
          	minimizeWeight(treemapPoints);
          }
        	console.log("# fix: "+fixCount);
        }
      }
      
      // heuristics: increase light weights
      function handleOverweighted1(treemapPoints) {
      	var fixCount = 0;
        var fixApplied, tpi, tpj, weightest, lightest, sqrD, overweight;
        do {
          fixApplied = false;
          for(var i=0; i<siteCount; i++) {
            tpi = treemapPoints[i];
            for(var j=i+1; j<siteCount; j++) {
              tpj = treemapPoints[j];
              if (tpi.weight > tpj.weight) {
                weightest = tpi;
                lightest = tpj;
              } else {
                weightest = tpj;
                lightest = tpi;
              }
              sqrD = squaredDistance(tpi, tpj);
              if (sqrD < weightest.weight-lightest.weight) {
                overweight = weightest.weight - lightest.weight - sqrD
                lightest.weight += overweight + epsilon;
                fixApplied = true;
                fixCount++;
                break;
              }
            }
            if (fixApplied) { break; }
          }
        } while (fixApplied)
        
        if (fixCount>0) {
          if (shouldMinimizeWeight) {
          	minimizeWeight(treemapPoints);
          }
        	console.log("# fix: "+fixCount);
        }
      }
      
      function minimizeWeight(treemapPoints) {
        var minWeight = treemapPoints[0].weight;
        
        for (var i=1; i<siteCount; i++) {
          minWeight = Math.min(minWeight, treemapPoints[i].weight);
        }
        minWeight -= epsilon;
        for (var i=0; i<siteCount; i++) {
          treemapPoints[i].weight -= minWeight;
        }
      }
      
      function squaredDistance(s0, s1) {
        return sqr(s1.x - s0.x) + sqr(s1.y - s0.y);
      };

};
      
      function handleNearZeorWeights(basePoints) {
        var maxWeight = basePoints.reduce(function(max, bp){
          return Math.max(max, bp.weight);
        }, -Infinity);
        var minAllowedWeight = maxWeight*nearZeroWeightRatio,
            nearZeroCount = 0;
        
        basePoints.forEach(function(bp) {
          if (bp.weight<minAllowedWeight) {
            bp.weight = minAllowedWeight;
            nearZeroCount++;
          }
        })
        
        if (nearZeroCount>0) {
          console.log("# near-zero weights: "+nearZeroCount);
        }
      }
      
      function createTreemapPoints(basePoints, initialPlacementPolicy) {
        var totalWeight = basePoints.reduce(function(acc, bp){ return acc+=bp.weight; }, 0);
        var avgWeight = totalWeight/siteCount;
        		avgArea = totalArea/siteCount,
        		defaultWeight = avgArea/2;
        
          
        if (initialPlacementPolicy === 'sortedPie') {
          // sortedPie ensures :
          //	- a gradient from light weights to heavy weights
          //  - i.e., light weights will not be confine between heavy weights
          //	- i.e., light weights will move/re-weight more easily
          var sortedBasePoints = basePoints.sort(function(bp0, bp1){
                return bp0.weight < bp1.weight;
              });
          
          return createTreemapPoints(sortedBasePoints, 'pie');
        }
        else if (initialPlacementPolicy === 'pie') {
          var deltaRad = _2PI/siteCount;
          var rad;
          
					return basePoints.map(function(bp, i) {
            rad = deltaRad*i;

return {
              index: bp.index,
              targetedArea: totalArea*bp.weight/totalWeight,
              data: bp,
              x: radius+halfRadius*Math.cos(rad)+Math.random(),
              y: radius+halfRadius*Math.sin(rad)+Math.random(),
              weight: defaultWeight
            }
          })
        } else {
          var xExtent = d3.extent(clippingPolygon.map(function(p){return p[0];})),
          		yExtent = d3.extent(clippingPolygon.map(function(p){return p[1];})),
              dx = xExtent[1]-xExtent[0],
              dy = yExtent[1]-yExtent[0];
          var x,y;
          
          return basePoints.map(function(bp) {
            //use (x,y) instead of (r,a) for a better uniform placement of sites (ie. less centered)
            x = xExtent[0]+dx*Math.random();
            y = yExtent[0]+dy*Math.random();
            while (!d3.polygonContains(clippingPolygon, [x, y])) { 
              x = xExtent[0]+dx*Math.random();
              y = yExtent[0]+dy*Math.random();
            }

return {
              index: bp.index,
              targetedArea: totalArea*bp.weight/totalWeight,
              data: bp,
              x: x,
              y: y,
              weight: defaultWeight
            }
          })
        }
      }
      
      setAdaptPlacements();
      setAdaptWeights();
      setHandleOverweighted();
      reset();
      
      /********************************/
      /*      Drawing functions       */
      /*   Playing with canvas :-)    */
      /*                              */
      /* Experiment to draw           */
      /*   with a uniform color,      */
      /*   or with a radial gradient, */
      /*   or over a background image */
			/********************************/
      
      function initLayout() {
        d3.select("#layouter").style("width", totalWidth+"px").style("height", totalHeight+"px");
        d3.selectAll("canvas").attr("width", width).attr("height", height);

bgImgCanvas = document.querySelector("canvas#background-image");
        bgImgContext = bgImgCanvas.getContext("2d");
        alphaCanvas = document.querySelector("canvas#alpha");
        alphaContext = alphaCanvas.getContext("2d");
        coloredCanvas = document.querySelector("canvas#colored");
        coloredContext = coloredCanvas.getContext("2d");

//begin: set a radial rainbow
        radialGradient = coloredContext.createRadialGradient(radius, radius, 0, radius, radius, radius);
        var gradientStopNumber = 10,
            stopDelta = 0.9/gradientStopNumber;
            hueDelta = 360/gradientStopNumber,
            stop = 0.1,
            hue = 0;
        while (hue<360) {
          radialGradient.addColorStop(stop, d3.hsl(Math.abs(hue+160), 1, 0.45));
          stop += stopDelta;
          hue += hueDelta;
        }
        //end: set a radial rainbow
        
        //begin: set the background image
        resetBackgroundImage();
        //end:  set the initial background image
      }
      
      function resetBackgroundImage() {
        if (bgType==="canonicalRainbow") {
          //begin: make conical rainbow gradient
          var imageData = bgImgContext.getImageData(0, 0, 2*radius, 2*radius);
          
          var i = -radius,
              j = -radius,
              pixel = 0,
              radToDeg = 180/Math.PI;
          var aRad, aDeg, rgb;
          while (i<radius) {
            j = -radius;
            while (j<radius) {
              aRad = Math.atan2(j, i);
              aDeg = aRad*radToDeg;
              rgb = d3.hsl(aDeg, 1, 0.45).rgb();
              
              imageData.data[pixel++] = rgb.r;
              imageData.data[pixel++] = rgb.g;
              imageData.data[pixel++] = rgb.b;
              imageData.data[pixel++] = 255;
              
              j++;
            }
            i++;
          }
          bgImgContext.putImageData(imageData, 0, 0);
          //end: make conical rainbow gradient
        } else if (bgType==="radialRainbow") {
          bgImgContext.fillStyle = radialGradient;
          bgImgContext.fillRect(0, 0, width, height);
        } else {
          bgImgContext.fillStyle = "grey";
        	bgImgContext.fillRect(0, 0, width, height);
        }
      }

Updated missing url https://raw.githack.com/Kcnarf/d3-weighted-voronoi/master/build/d3-weighted-voronoi.js to https://raw.githack.com/kcnarf/d3-weighted-voronoi/master/build/d3-weighted-voronoi.js

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

https://raw.githack.com/Kcnarf/d3-weighted-voronoi/master/build/d3-weighted-voronoi.js