在数据可视化中，我们往往会使用图来表达数据中所蕴含的信息。而图布局算法可以使散乱的信息 (信息多以点线的关系承载) 通过一种清晰的方式呈现出来，并符合相应的美学标准。在图布局算法模型中，其建立在粒子物理理论的基础上，将节点模拟成为原子，通过原子间的引力和斥力来得到节点的速度与加速度，计算其移动方位与距离，最终达到一个稳定平衡的状态，从而完成布局。以下就是由 d3 实现的力引导布局:
在 d3 的实现中，为了达到性能与效果的平衡，节点与节点间模拟同种电荷相互排斥，并将节点存入四叉树中，利用 Barnes–Hut 近似来减少节点间电荷斥力的计算量。同时连线间的节点模拟弹簧牵引力，节点的速度综合斥力引力得出，并发生阻尼衰减，最终达到整图平衡。
在本文中，我们将对 d3 实现的力导引布局进行一步步分解，详细剖析其实现过程与背后的原理。在此之前，读者们自行阅读上图实现源码，熟悉一下 d3-force API。
节点的处理与优化
初始化导入节点
首先，需要将节点形成所符合的数据结构导入 d3 当中，对于每一个节点进行预处理，节点按一定的半径和旋转角度环绕起来，vx 与 vy 分别表示节点在 x 轴和 y 轴方向上的速度分量。
| 1 2 3 4 5 6 7 8 9 10 11 | var initialRadius = 10, initialAngle = Math.PI * (3 - Math.sqrt(5)); for (var i = 0, n = nodes.length, node; i < n; ++i) { node = nodes[i], node.index = i; ... var radius = initialRadius * Math.sqrt(i), angle = i * initialAngle; node.x = radius * Math.cos(angle); node.y = radius * Math.sin(angle); ... node.vx = node.vy = 0; } |
建立节点四叉树
四叉树 (Q-Tree) 是一种树形数据结构。它的每个节点下至多可以有四个子节点，通常把一部分二维空间细分为四个象限并把该区域里的相关信息存入到四叉树节点中。四叉树的每一个节点代表一个矩形区域，每一个矩形区域又可划分为四个小矩形区域，这四个小矩形区域作为四个子节点所代表的矩形区域，正如下图的矩阵虚线划分:
遍历导入所有节点数据，求出节点最小 x 值 x0x0 ，最大 x 值 x1x1，最小 y 值 y0y0，最大 y 值 y1y1，将 (x0,y0)(x0,y0) 点坐标以 22 倍增，使 x0∗2n≥x1,y0∗2n≥y1x0∗2n≥x1,y0∗2n≥y1，那么新 (x1,y1)(x1,y1) 为 (x0∗2n,y0∗2n)(x0∗2n,y0∗2n)。以 22 倍增符合四叉树在一维上不断对半划分的特点。
一边添加节点，一边对四叉树进行划分，因此一般来说四叉树是不完全的。在添加时会遇到如下 4 种情况:
case1: 当前四叉树为空，直接将当前节点作为树的根节点。
case2: 当前查询节点为真实节点(对象节点)，即与添加节点所处范围矩阵一致，需要建立数组索引，即再次划分该矩阵，直到查询节点与添加节点分处不同矩阵。
case3: 当前查询节点为索引节点(数组节点)，且添加节点所处范围矩阵为空，直接添加。
case4: 当前查询节点为真实节点，且添加节点的坐标与其完全相等，直接下挂。
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | function add(tree, x, y, d) { var node = tree._root, leaf = {data: d}, x0 = tree._x0, y0 = tree._y0, x1 = tree._x1, y1 = tree._y1, ... // 当前tree._root为空，对于case1 // If the tree is empty, initialize the root as a leaf. if (!node) return tree._root = leaf, tree; // 类似与二分查找，自顶向下搜索 // Find the existing leaf for the new point, or add it. while (node.length) { if (right = x >= (xm = (x0 + x1) / 2)) x0 = xm; else x1 = xm; if (bottom = y >= (ym = (y0 + y1) / 2)) y0 = ym; else y1 = ym; // 判断当前添加节点所处象限是否为空，对应case3 if (parent = node, !(node = node[i = bottom << 1 | right])) return parent[i] = leaf, tree; } // node为查询节点 // Is the new point is exactly coincident with the existing point? xp = +tree._x.call(null, node.data); yp = +tree._y.call(null, node.data); // 特殊情况，若当前需要加入的节点与父节点完全重合，对应case4 if (x === xp && y === yp) return leaf.next = node, parent ? parent[i] = leaf : tree._root = leaf, tree; // 不停分割，直至处于不同象限，对应case2 // Otherwise, split the leaf node until the old and new point are separated. do { parent = parent ? parent[i] = new Array(4) : tree._root = new Array(4); if (right = x >= (xm = (x0 + x1) / 2)) x0 = xm; else x1 = xm; if (bottom = y >= (ym = (y0 + y1) / 2)) y0 = ym; else y1 = ym; } while ((i = bottom << 1 | right) === (j = (yp >= ym) << 1 | (xp >= xm))); return parent[j] = node, parent[i] = leaf, tree; } |
四叉树的构建完成后，自四叉树由下而上，求所有节点的合坐标与合静电电荷量。若当前节点 quadquad 为索引节点，其中 strength 为静电电荷量，默认值均为 -30，可由开发者自定义。
合坐标公式: (quad.x=∑i=03(quad[i].value∗quad[i].x)∑i=03(quad[i].value),quad.y=∑i=03(quad[i].value∗quad[i].y)∑i=03(quad[i].value))(quad.x=∑i=03(quad[i].value∗quad[i].x)∑i=03(quad[i].value),quad.y=∑i=03(quad[i].value∗quad[i].y)∑i=03(quad[i].value))
合电荷量公式: ∑i=03(quad[i].value)∑i=03(quad[i].value)
若当前节点为真实节点，坐标即为当前真实节点坐标，静电电荷量则为真实节点及其下挂真实节点 strength 总和。
斥力的优化求解
节点间的关键就在于电荷斥力的求解。在原来的朴素算法中，需要每个节点对其他所有节点求解斥力，算法复杂度为 O(n2)O(n2)，而目前采用四叉树与 Barnes-Hut 近似，算法复杂度 O(nlogn)O(nlogn)，其核心思想在于，当前节点 (node) 计算远处节点的斥力对速度位移影响时，可将其周围相近的点产生的相同斥力效果做整合处理，形成”质心”，而这个”周围”的大小就由四叉树的矩阵与 Barnes-Hut 近似精度 theta 所决定，theta 默认值为 (0.9)2(0.9)2。
当前象限区域面积 S1S1 为 (x1−x2)2(x1−x2)2，当前节点 (node) 与象限合节点 (quad) 形成矩阵面积 S2S2 为 (quad.x−node.x)∗(quad.y−node.y)(quad.x−node.x)∗(quad.y−node.y)，当 S1S2<thetaS1S2<theta 时，计算当前节点与象限合节点的相互作用力，将其转化为 node 速度变化量:
根据 Velocity Verlet 算法， v(t+Δt)=v(t)+f(t+Δt)+f(t)2mΔtv(t+Δt)=v(t)+f(t+Δt)+f(t)2mΔt，d3 进行了简化处理，设 m=1,Δt=1m=1,Δt=1:
Δv=v(t+Δt)−v(Δt)=f(t+Δt)+f(t)2mΔt≈f(t)Δv=v(t+Δt)−v(Δt)=f(t+Δt)+f(t)2mΔt≈f(t)
而在 d3 源码中，node.vx+=(x∗quad.valuex2+y2)∗alphanode.vx+=(x∗quad.valuex2+y2)∗alpha，alpha 是阻尼衰减系数。暂时无法理解这个求解过程，按照电场力求解公式 F=kQqr2F=kQqr2 推导，再将力分解到 x 轴，与现有求解不符，已提 issue 给代码的作者，希望能得到他的解答。
具体代码如下，apply 函数第一个参数 quad 为四叉树索引节点，内有索引下属子节点的合坐标 (quad.x,quad.y)(quad.x,quad.y) 和合电荷量。其返回 true 意味着，当前节点及其子节点已完成计算，否则需要继续向下遍历节点。
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | var distanceMin2 = 1, distanceMax2 = Infinity, theta2 = 0.81; function apply(quad, x1, _, x2) { if (!quad.value) return true; var x = quad.x - node.x, y = quad.y - node.y, w = x2 - x1, l = x * x + y * y; // 是否可采用 Barnes-Hut 近似 // Apply the Barnes-Hut approximation if possible. // Limit forces for very close nodes; randomize direction if coincident. if (w * w / theta2 < l) { if (l < distanceMax2) { if (x === 0) x = jiggle(), l += x * x; if (y === 0) y = jiggle(), l += y * y; if (l < distanceMin2) l = Math.sqrt(distanceMin2 * l); node.vx += x * quad.value * alpha / l; node.vy += y * quad.value * alpha / l; } return true; } // 无法采用 Barnes-Hut 近似且 quad 有节点，或 l 大于距离上限，需要继续向下遍历 // Otherwise, process points directly. else if (quad.length || l >= distanceMax2) return; // 排除自身对自身影响 // Limit forces for very close nodes; randomize direction if coincident. if (quad.data !== node || quad.next) { if (x === 0) x = jiggle(), l += x * x; if (y === 0) y = jiggle(), l += y * y; if (l < distanceMin2) l = Math.sqrt(distanceMin2 * l); } do if (quad.data !== node) { w = strengths[quad.data.index] * alpha / l; node.vx += x * w; node.vy += y * w; } while (quad = quad.next); } |
节点连线的处理
相比节点的处理而言，节点的连线简单了许多，没有使用任何优化，先初始化连线，统计每个节点的度，求每一条边的起点 (source) 度的占比，使 bias = 起点度 / (起点度 + 终点度)。每条边的默认长度 (distance) 为30，默认弹簧劲度系数 (strength) 为 1 / min(起点度, 终点度)，这是为了减小对于度较大节点的引力，提高稳定性。
遍历所有连线，计算施加在连线两端节点的引力，最终推导出速度的变化:
target.vx−=k∗Δx∗cosθ∗alpha∗biastarget.vx−=k∗Δx∗cosθ∗alpha∗bias
=strengths[i]∗(l−distances[i])∗xx2+y2√∗alpha∗bias=strengths[i]∗(l−distances[i])∗xx2+y2∗alpha∗bias
target.vy、source.vx、source.vy 同理。
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | var iterations = 1; function force(alpha) { for (var k = 0, n = links.length; k < iterations; ++k) { for (var i = 0, link, source, target, x, y, l, b; i < n; ++i) { link = links[i], source = link.source, target = link.target; x = target.x + target.vx - source.x - source.vx || jiggle(); y = target.y + target.vy - source.y - source.vy || jiggle(); l = Math.sqrt(x * x + y * y); l = (l - distances[i]) / l * alpha * strengths[i]; x *= l, y *= l; target.vx -= x * (b = bias[i]); target.vy -= y * b; source.vx += x * (b = 1 - b); source.vy += y * b; } } } |
布局的形成
布局的形成主要依靠不断的迭代计算，每处理一次节点与节点连线称为一步，循环往复，通过 alpha 逐渐衰减与 (fx, fy) 来进行控制。
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 | var simulation, alpha = 1, alphaMin = 0.001, // alpha衰减率 alphaDecay = 1 - Math.pow(alphaMin, 1 / 300), alphaTarget = 0, // 速度衰减 velocityDecay = 0.6, stepper = timer(step), // tick事件与end事件 event = dispatch("tick", "end"); function step() { tick(); event.call("tick", simulation); if (alpha < alphaMin) { stepper.stop(); event.call("end", simulation); } } function tick() { // alpha不断衰减 alpha += (alphaTarget - alpha) * alphaDecay; // 不停迭代 forces.each(function(force) { force(alpha); }); // 速度转化为距离改变 for (i = 0; i < n; ++i) { node = nodes[i]; if (node.fx == null) node.x += node.vx *= velocityDecay; // 具有fx，说明当前节点被控制，不需要受到力的影响，速度置为0 else node.x = node.fx, node.vx = 0; if (node.fy == null) node.y += node.vy *= velocityDecay; else node.y = node.fy, node.vy = 0; } } |
小结
d3-force 的实现与传统图力导引布局中的 FR 算法具有完全相同的算法思路，均将迭代分为三部分，先计算节点之间相互的排斥力，然后是计算图中有边连接的节点之间相互的吸引力，最后综合吸引力和排斥力，来得到每个节点的速度，而不是加速度。同时采用类似于模拟退火的衰减方案，使布局趋于稳定。
在力的求解方式上，d3 的实现与 FR 的求解有所不同，对于 FR 中引力和斥力的公式推导也暂时没有想明白，这是一个优化点。此外，d3 中包括节点的静电电荷量，连线弹簧的长度、劲度系数，连线引力的迭代次数目前都采用了默认值，可否进行调整，这也是一个优化点。
参考资料
d3-force doc
d3-quadtree doc
Velocity_Verlet wiki
A Multilevel Algorithm for Force-Directed Graph Drawing
图布局力导引算法研究与实现
基于Barnes Hut算法的N-body问题模拟
出处：
914
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
