We use the Javascript framework GPU.js to compile custom WebGL shaders to render and animate the heatmap together with normal Wolfram Language graphics primitives. Based on the blog here.

All dependencies are kept within this notebook and do not require installation or compilation.

Building common module (optional)

We build a simple module, which will hook up a CommonJS module to a global variable so that we can continue experimenting with compute shader code without needing to rebuild the whole thing.

.esm    
window.gpuJS = require('./js/gpu.js');

Please, create a folder in your notebook directory js/ and place there a fresh build of gpu.js or download an a self-extracting notebook from the original blog-page.

Now we have an easy access to this framework.

Compute shaders

Here we use GPU.js to write a compute shader in Javascript, which use privded array of 2D points and renders their density to a canvas. We shall support the possibility to update the image in real-time and use PlotRange, ImageSize parameters similar to Plot function

.js
function getColor(t) {     
    if (t >= 0 && t <= 0.25)       
        return [255.0 + (- 255.0) * t / 0.25, 255.0 + ( - 255.0) * t / 0.25, 255.0];     
    if (t >= 0.25 && t <= 0.55)       
        return [0.0, 255.0 * (t - 0.25) / 0.3, 255.0 + (- 255.0) * (t - 0.25) / 0.3];     
    if (t >= 0.55 && t <= 0.85)       
        return [255.0 * (t - 0.55) / 0.3, 255.0, 0.0];     
    if (t >= 0.85 && t <= 1)       
        return [255.0, 255.0 + ( - 255.0) * (t - 0.85) / 0.15, 0.0];     
    return [255.0, 255.0, 255.0];   
}

core.HeatMap = async (args, env) => {
  const gpu = new window.gpuJS.GPU();
  const options = await core._getRules(args, env);

  env.local.gpu = gpu;

  let [width, height] = [320, 240];
  let plotRange = [[-1,1], [-1,1]];

  if (options.ImageSize) {
    if (Array.isArray(options.ImageSize))
      [width, height] = options.ImageSize;
    else
      [width, height] = [Math.round(options.ImageSize), Math.round(options.ImageSize/1.6180339)];
  }

  if (options.PlotRange) {
    plotRange = options.PlotRange;
  }

  plotRange = [
    plotRange[0][0],
    width / (plotRange[0][5] - plotRange[0][0]),

    plotRange[1][0],
    height / (plotRange[1][6] - plotRange[1][0])
  ];


  const points = await interpretate(args[0], env);
  let alpha = 10000.0;
  if (args.length - Object.keys(options).length > 1) {
    alpha = await interpretate(args[1], env);
    env.local.alpha_has = true;
  }

  env.local.alpha = alpha;

  env.local.kernel = gpu.createKernel(function(points, pointCount, alpha) {
    let intensity = 0;
    for (let i = 0; i < pointCount; i++) {
        const pointX = (points[i * 2] - this.constants.plotRange[0]) * this.constants.plotRange[1];
        const pointY = (points[i * 2 + 1] - this.constants.plotRange[2]) * this.constants.plotRange[3];
        const dx = (this.thread.x - pointX)/this.constants.width;
        const dy = (this.thread.y - pointY)/this.constants.height;
        const dist = dx * dx + dy * dy;
        intensity += 1.0/(dist * alpha + 1.0); // Blurry effect
    }

    intensity = Math.min(intensity, 1); // Cap intensity
    const baseColor = getColor(intensity);

    let blendedColor = [
        baseColor[0] * intensity + (1.0 - intensity)*255,
        baseColor[1] * intensity + (1.0 - intensity)*255,
        baseColor[2] * intensity + (1.0 - intensity)*255
    ];
    this.color(blendedColor[0] / 255, blendedColor[1] / 255, blendedColor[2] / 255);
  }, {
    dynamicArguments: true,
    constants: {
      plotRange: plotRange,
      width: width,
      height: height
    }
  })
  .setOutput([width, height])
  .setGraphical(true)
  .setFunctions([getColor]);  

  env.element.appendChild(env.local.kernel.canvas);

  if (points instanceof NumericArrayObject) {
    env.local.kernel(new Float32Array(points.buffer), points.dims[0], alpha);
    return;
  }

  env.local.kernel(points.flat(), points.length, alpha);
}

core.HeatMap.virtual = true

core.HeatMap.destroy = (args, env) => {
  env.local.kernel.destroy();
  env.local.gpu.destroy();
  delete env.local.gpu;
  delete env.local.kernel;
}

core.HeatMap.update = async (args, env) => {
  let alpha = env.local.alpha;
  if (env.local.alpha_has) {
    alpha = await interpretate(args[1], env);
  }

  const data = await interpretate(args[0], env);
  if (data instanceof NumericArrayObject) {
    env.local.kernel(new Float32Array(data.buffer), data.dims[0], alpha);
    return;
  }

  env.local.kernel(data.flat(), data.length, alpha);
}

this.ondestroy(() => {
  gpu.destroy();
});

return 'If you see it, it works'

Output form

To display the heatmap, we need to define an output form in Wolfram Language, that will use Javascript to render the widget

HeatMap /: MakeBoxes[h_HeatMap, StandardForm] := With[{},
  ViewBox[h, h]
]

HeatMap /: MakeBoxes[h_HeatMap, StandardForm] := With[{
  o = CreateFrontEndObject[h]
},
  MakeBoxes[o, StandardForm]
] /; ByteCount[h] > 1024

Options[HeatMap] = {ImageSize->{320, 240}, PlotRange->{{-1,1}, {-1,1}}};

Test

Create a bunch of points and render it

HeatMap[Exp[-Norm[#]]# &/@ RandomReal[{-Pi,Pi}, {30,2}]]

Integration with Graphics

Using basic Inset we can place it under any graphics

range = {{-1,1}, {-1,1}}/2.0;
dataset = Exp[-Norm[#]]# &/@ RandomReal[{-Pi,Pi}, {50,2}];

Graphics[{
  ColorData[97][8], PointSize[0.03], 
  Inset[
    HeatMap[dataset, 4000, PlotRange->range, ImageSize->{344,244}],
    (range[[All,1]]+range[[All,2]])/2.0
  ],
  Point[dataset]
}, 
  PlotRange->range, ImageSize->{420,320}, ImagePadding->30,
  Frame->True, Axes->True
]

Make it dynamic

We should take a full advantage of the WebGL and perform smooth animation of the data points. For this we define helper function to guide our points though the calculated path

fixedPoint[args__, maxN_:Infinity] := NestWhile[Function[input,
  With[{
    pos = input[[1]], 
    target = input[[2]], 
    velocity = input[[3]], 
    iteration = input[[4]]
  }, {
    pos + 0.01 velocity,
    target,
    velocity 0.9 + 0.1 (Normalize /@ ( target - pos)),
    iteration + 1
  }]
], {args}, (Max[Abs[#[[1]]-#[[2]]]]>0.1 && #[[4]] <= maxN)&][[{1,4}]]

Now we generate the data, find the number of interations steps until all points converge the target positions and assign a slider to control the process

(* input data *)
range = {{-1,1}, {-1,1}};
points = RandomPoint[RegionDifference[Disk[{0,0}, 0.6], Disk[{0,0}, 0.5]], 100];

(* place points outside the viewport *)
offscreenPos = With[{corners = range // Transpose}, 
  RandomPoint[
    RegionDifference[
      Rectangle @@ (2 corners),
      Rectangle @@ corners
    ], Length[points]
  ]
];

(* procedurally animate each particle *)
With[{
  vel = 0.1 RandomReal[{-1,1}, {Length[points], 2}],
  target = points
},
  points = offscreenPos;

  With[{
    maxIterations = fixedPoint[offscreenPos, target, vel, 0, Infinity][[2]]
  },
    EventHandler[InputRange[0, maxIterations, 1, 0], Function[i,
      points = fixedPoint[offscreenPos, target, vel, 0, i][[1]]
    ]]
  ]
]

(* render *)
Graphics[{
  ColorData[97][10], PointSize[0.03], 
  Inset[
    HeatMap[points // Offload, 4000, PlotRange->range, ImageSize->{344,244}],
    (range[[All,1]]+range[[All,2]])/2.0
  ],
  Point[points // Offload]
}, 
  PlotRange->range, ImageSize->{420,320}, ImagePadding->30,
  Frame->True, Axes->True, TransitionType->None, "Control"->False
]