气泡词云图空间利用率最大化解决方案
之前有一个需求，需要把珊瑚图表的气泡词云图空间利用率最大化，一开始没有思路，研究的比较久，最后勉强实现了。涉及的知识点也比较多，这个算是我写代码以来，用到的数学、物理知识最深的一次了，所以在这里记录一下。
1 需求描述
气泡词云图之前使用d3-hierarchy实现的，该实现是一个圆形，如果在一个长宽不等的矩形容器中，空间利用率较低，会出现大片留白，希望实现为尽可能填满容器，将空间利用率最大化。
2 问题本质
之前没有遇到过类似的问题，没什么思路，经过不懈的努力，终于发现，这个问题的本质其实是--不等圆packing，不等圆packing问题是如何将N个任意半径的圆形物体互不相嵌的放入一个圆形容器中，找出一种布局使得这个圆形容器的半径最小。
不等圆packing问题具有非常广泛的应用场景，如裁缝如何最大效率的利用布料、设计航天器时需要把若干个半径给定的圆柱形部件放入一个尽可能小的圆柱形腔体内，此外，在无线通信、汽车工业材料切割等工业领域也有着广泛的应用。
该问题是一个典型的NP难度问题，使用传统的数学方法很难求解，迄今为止尚未设计出既严格（确保能找到全局最优解）又迅速（多项式计算复杂度 ）的求解算法。现有的求解算法几乎均为启发式算法，主要包括两大类，一是基于构造规则的构造算法，如占角算法（最大穴度算法）、PERM、Beam Search、Algorithm2等，二是基于演化规则的演化算法，如拟人拟物算法、禁忌搜索算法、粒子群算法、遗传算法、模拟退火算法、以及混合算法SATS等。
3 算法分析
算法主要分两种，构造型和演化型，本文选取其中的占角算法、拟人拟物算法、分子动力学模拟方法进行分析。
3.1 占角算法（最大穴度算法）
在不考虑n个圆的各种排列情况下首先将第一个圆放进与大圆R相切的一个位置，接着把第二个圆放到与前两个圆相切的位置之一，然后把第三个圆也放置到与前三个圆中的两个圆相切的位置之一，就这样将下一个圆放到前i个圆中任意两个圆相切的位置，如果这个位置导致了新放入的圆与其他圆有重叠则将新圆放置在满足与其他两个圆相切的另一位置，如果发现所有位置都不满足不重叠的条件，则回溯到之前的一个圆，并重新选取之前圆的位置。按照这样的规则满足条件则继续添加新圆，如果不满足则回溯到上一个圆放置的位置，直到所有的圆都放入到了大圆R中，则找到问题的解，或者遍历所有情况仍旧找不到解，则认为该问题没有解。
d3的打包图就是使用占角算法，首先将第一个圆放到原点，然后将第二个圆与第一个圆相切，将其居中：
// Place the first circle.
a = circles[0], a.x = 0, a.y = 0;
if (!(n > 1)) return a.r;
// Place the second circle.
b = circles[1], a.x = -b.r, b.x = a.r, b.y = 0;
if (!(n > 2)) return a.r + b.r;然后将第三个圆放到与两个圆同时相切的位置（占角），可用余弦定理或者两点间的距离公式求得第三个圆的坐标：
function place(b, a, c) {
var dx = b.x - a.x, x, a2,
dy = b.y - a.y, y, b2,
d2 = dx * dx + dy * dy;
if (d2) {
a2 = a.r + c.r, a2 *= a2;
b2 = b.r + c.r, b2 *= b2;
if (a2 > b2) {
x = (d2 + b2 - a2) / (2 * d2);
y = Math.sqrt(Math.max(0, b2 / d2 - x * x));
c.x = b.x - x * dx - y * dy;
c.y = b.y - x * dy + y * dx;
} else {
x = (d2 + a2 - b2) / (2 * d2);
y = Math.sqrt(Math.max(0, a2 / d2 - x * x));
c.x = a.x + x * dx - y * dy;
c.y = a.y + x * dy + y * dx;
}
} else {
c.x = a.x + c.r;
c.y = a.y;
}
}计算过程中防止重叠，还需要判断两圆之间的碰撞，碰撞检测使用两圆圆心的距离小于等于两圆的半径之和：
function intersects(a, b) {
var dr = a.r + b.r - 1e-6, dx = b.x - a.x, dy = b.y - a.y;
return dr > 0 && dr * dr > dx * dx + dy * dy;
}全部小圆放入之后再根据容器较短边等比例缩放，然后同时移动至容器中心。
效果如下：
这个过程一直是寻找最大凹点，也就是最大穴度，所以最终会是一个像圆的形状，要控制矩形边界的话不太容易，经过多次尝试，还是不能控制矩形边界。
3.2 拟人拟物算法
这是一个求解不等圆packing 问题效率较高的方法。拟物算法将大圆R看做封闭刚性不可改变空的圆形容器，将每个小圆看做光滑弹性的小球，开始将这n个小球强行放到容器中，我们定义小球之间因为弹性挤压产生的能量为弹性势能（U=U1+U2+···+Un),接着小球由于相互间以及和容器间的作用力不断运动，每个小球都向受挤压的方向运动，直至所有小球受力到达平衡或者不再受力。拟物方法就是模拟这样一个过程。
拟人方法是为了解决拟物方法中可能出现的一种“死锁”的情况，即n个小球在局部最小解的情况中反复移动无法跳出的情况（这时的弹性势能U不再变化且不为0),“人为地”从容器中取出最拥挤的小球重新放入容器内进行拟物计算，这样就解决了拟物算法可能遇到的僵持平衡的状态。整个过程直到所有的小球都不受到弹性力，弹性势能为零（U≈0)则认为找到了解。
代码实现：
// 开始计算数据
startSetData(circles, width, height, padding0, padding1) {
let step = 0.1;
let oldmaxPE = circles[0];
let maxPE = circles[1];
let minPE = circles[1];
let newPE;
let oldPE = 0;
let t = 0;
const num = circles.length;
while (true) {
if (step >= 0.001) {
// 计算所有圆的总势能，找出最大势能圆和最小势能圆
const peObj = this.calcPE(circles, maxPE, minPE);
newPE = peObj.PE;
maxPE = peObj.maxPE;
minPE = peObj.minPE;
// 总势能约等于0，结束循环
if (newPE <= 0.0001) {
break;
} else {
if (newPE >= oldPE) {
step = 0.8 * step;
}
oldPE = newPE;
for (let i = 0; i < num; i++) {
// 计算每个圆新的位置
circles[i].x -= step * circles[i].dx;
circles[i].y -= step * circles[i].dy;
// 边界检测
if (circles[i].x + circles[i].r >= width - padding0) {
circles[i].x = width - circles[i].r - padding0;
}
if (circles[i].x - circles[i].r <= padding0) {
circles[i].x = circles[i].r + padding0;
}
if (circles[i].y + circles[i].r >= height - padding1) {
circles[i].y = height - circles[i].r - padding1;
}
if (circles[i].y - circles[i].r <= padding1) {
circles[i].y = circles[i].r + padding1;
}
}
}
}
else {
if (maxPE == oldmaxPE) {
t = t + 1;
}
const NextDouble = Math.random();
if (t < 1) {
// 在圆盘上随机选点，重新摆放势能最大的圆
// 注意此时的maxPE 和 oldmaxPE均指向old_Circles中的对象。
maxPE.x = width * NextDouble;
maxPE.y = height * NextDouble;
oldmaxPE = maxPE;
step = 0.1;
}
else if (t == 1) {
minPE.x = width * NextDouble;
minPE.y = height * NextDouble;
t = 0;
oldmaxPE = circles[0];
step = 0.1;
}
}
}
}
// 计算势能
calcPE(circles, maxPE, minPE) {
let PE = 0;
const num = circles.length;
// 计算前，将以前的计算结果清零！！！
for (let i = 0; i < num; i++) {
circles[i].dx = 0;
circles[i].dy = 0;
circles[i].PE = 0;
}
// 逐个计算圆的势能以及移动量
for (let i = 0; i < num; i++) {
this.moveDirection(circles[i], num, circles, i);
}
// 累加势能
for (let i = 0; i < num; i++) {
PE += circles[i].PE;
}
// 找出最大势圆和最小势能圆
maxPE = minPE = circles[0];
for (let i = 1; i < num; i++) {
maxPE = maxPE.PE < circles[i].PE ? circles[i] : maxPE;
minPE = minPE.PE > circles[i].PE ? circles[i] : minPE;
}
return { PE, maxPE, minPE };
}
// 计算单个圆势能、移动量
moveDirection(circle, num, circles, numth) {
for (let i = 0; i < num; i++)
{
// 计算当前圆和其他圆之间的势能和移动量
if (numth != i) {
const dij = Math.sqrt(Math.pow(circle.x - circles[i].x, 2.0) + Math.pow(circle.y - circles[i].y, 2.0));
if (dij < circle.r + circles[i].r) {
if (circles[i].x === circle.x) {
circle.dx += circle.r + circles[i].r - dij;
} else {
circle.dx += (circles[i].x - circle.x) / dij * (circle.r + circles[i].r - dij);
}
if (circles[i].y === circle.y) {
circle.dy += circle.r + circles[i].r - dij
} else {
circle.dy += (circles[i].y - circle.y) / dij * (circle.r + circles[i].r - dij);
}
circle.PE += Math.pow(circle.r + circles[i].r - dij, 2.0);
}
}
}
}效果如下：
这种方法也有一些缺点，经过多次测试，发现有些数据可能出现计算很久的情况，还有一些数据无解，计算太久会卡住。
3.3 分子动力学模拟方法
以上两种方法均有不足之处，为了更加完美的解决问题，还可以使用分子动力学模拟的方法进行仿真模拟，它假设任意单位时间步长 Δt = 1，所有的粒子的单位质量常量 m = 1。作用在每个粒子上的合力 F 相当于在单位时间 Δt 内的恒定加速度 a。并且可以简单的通过为每个粒子添加速度并计算粒子的位置来模拟仿真。
d3-force提供了一个力模型，可以施加一些力，使圆之间不重叠，但是这个模型没有边界检测，需要新创建一种力，用于边界检测。
实现方式如下：
const x0 = boxPadding[0],
y0 = boxPadding[1],
x1 = svgWidth - boxPadding[0],
y1 = svgHieght - boxPadding[1] * 2;
const simulation = d3
.forceSimulation(rootData)
.force(
"collide",
d3
.forceCollide()
.radius(function (d) {
return d.r + 5;
})
.strength(2)
.iterations(2)
)
.force('charge', d3.forceManyBody().strength(20))
.force('boundary', forceBoundary(x0, y0, x1, y1));
simulation.tick(300);其中，用于边界检测的力模型如下：
function constant(x) {
return function() {
return x;
};
}
function forceBoundary(x0, y0, x1, y1) {
var strength = constant(0.1),
hardBoundary = true,
border = constant( Math.min((x1 - x0)/2, (y1 - y0)/2) ),
nodes,
strengthsX,
strengthsY,
x0z, x1z,
y0z, y1z,
borderz,
halfX, halfY;
if (typeof x0 !== "function") x0 = constant(x0 == null ? -100 : +x0);
if (typeof x1 !== "function") x1 = constant(x1 == null ? 100 : +x1);
if (typeof y0 !== "function") y0 = constant(y0 == null ? -100 : +y0);
if (typeof y1 !== "function") y1 = constant(y1 == null ? 100 : +y1);
function getVx(halfX, x, strengthX, border, alpha) {
return (halfX - x) * Math.min(2, Math.abs( halfX - x) / halfX) * strengthX * alpha;
}
function force(alpha) {
for (var i = 0, n = nodes.length, node; i < n; ++i) {
node = nodes[i];
if ((node.x - node.r <= (x0z[i] + borderz[i]) || node.x + node.r >= (x1z[i] - borderz[i])) ||
(node.y - node.r <= (y0z[i] + borderz[i]) || node.y + node.r >= (y1z[i] - borderz[i])) ) {
node.vx += getVx(halfX[i], node.x, strengthsX[i], borderz[i], alpha);
node.vy += getVx(halfY[i], node.y, strengthsY[i], borderz[i], alpha);
} else {
node.vx = 0;
node.vy = 0;
}
if (hardBoundary) {
if (node.x + node.r >= x1z[i]) node.vx += x1z[i] - node.x - node.r;
if (node.x - node.r <= x0z[i]) node.vx += x0z[i] - node.x + node.r;
if (node.y + node.r >= y1z[i]) node.vy += y1z[i] - node.y - node.r;
if (node.y - node.r <= y0z[i]) node.vy += y0z[i] - node.y + node.r;
}
}
}
function initialize() {
if (!nodes) return;
var i, n = nodes.length;
strengthsX = new Array(n);
strengthsY = new Array(n);
x0z = new Array(n);
y0z = new Array(n);
x1z = new Array(n);
y1z = new Array(n);
halfY = new Array(n);
halfX = new Array(n);
borderz = new Array(n);
for (i = 0; i < n; ++i) {
strengthsX[i] = (isNaN(x0z[i] = +x0(nodes[i], i, nodes)) ||
isNaN(x1z[i] = +x1(nodes[i], i, nodes))) ? 0 : +strength(nodes[i], i, nodes);
strengthsY[i] = (isNaN(y0z[i] = +y0(nodes[i], i, nodes)) ||
isNaN(y1z[i] = +y1(nodes[i], i, nodes))) ? 0 : +strength(nodes[i], i, nodes);
halfX[i] = x0z[i] + (x1z[i] - x0z[i])/2,
halfY[i] = y0z[i] + (y1z[i] - y0z[i])/2;
borderz[i] = +border(nodes[i], i, nodes)
}
}
force.initialize = function(_) {
nodes = _;
initialize();
};
force.x0 = function(_) {
return arguments.length ? (x0 = typeof _ === "function" ? _ : constant(+_), initialize(), force) : x0;
};
force.x1 = function(_) {
return arguments.length ? (x1 = typeof _ === "function" ? _ : constant(+_), initialize(), force) : x1;
};
force.y0 = function(_) {
return arguments.length ? (y0 = typeof _ === "function" ? _ : constant(+_), initialize(), force) : y0;
};
force.y1 = function(_) {
return arguments.length ? (y1 = typeof _ === "function" ? _ : constant(+_), initialize(), force) : y1;
};
force.strength = function(_) {
return arguments.length ? (strength = typeof _ === "function" ? _ : constant(+_), initialize(), force) : strength;
};
force.border = function(_) {
return arguments.length ? (border = typeof _ === "function" ? _ : constant(+_), initialize(), force) : border;
};
force.hardBoundary = function(_) {
return arguments.length ? (hardBoundary = _, force) : hardBoundary;
};
return force;
}这种方法实现效果和拟人拟物效果相似，计算速度更快，也不会出现无解的情况，算是比较完美的一种，所以最终选用了这种方式勉强完成了本次需求。
4 总结
以前在学校学习数学、物理知识的时候总感觉好像用不到，但是从这次需求中发现，没有这些基本的知识，要解决问题是比较困难的，特别是在数据可视化中，之前一直做一些比较基础的前端工作，会产生一种错觉，前端就是新时代的搬砖工，只能做基础的搬砖工作，接触到图表之后才发现，这里面水很深，还是有很多算法没有js版本，前端在算法领域成熟度还是有待提高。虽然在平时使用不到，但是真正用的时，如果不能理解算法底层原理，要解决这种类似的问题是比较困难的，甚至无从下手。