N体受力问题（一）：四叉树

最近由于项目需要，研究了下力导向图的实现，发现了不少有意思的东西。如力导向图的弹性布局效果的实现是将节点看成带电粒子，粒子因为受到其他粒子的引力和斥力发生运动，最终受力平衡而稳定下来。

要计算连接粒子间的受力不算难，根据距离套用公式即可。但是单个粒子不只是受到与其连接的粒子的力，而是受到来自除该粒子外所有粒子的力，即任何两个粒子之间都将受到力的影响。这种力就叫N体力(Many-Body)。

基本思想

假如一个力导向图有100个节点，要计算每个粒子的N体力，就要进行100*100 = 10000次计算，这是极其耗时的。所以要快速地计算N体力，需要减少问题的规模，提升N体问题模拟算法的速度，一种非常重要的思想就是把相互接近的一组物体近似看成单独的一个物体。对于一组距离足够远的物体，我们可以将它的力作用近似看成是来自其质心的作用力。一组物体的质心是这种物体经过质量加权后的平均位置。例如，如果两个物体的位置分别是(x1,y1)和(x2,y2)，质量分别为m1和m2，那么它们的总质量和总质心(x,y)分别为：

将平面进行等分划分是一种将物体分组的好方法，在同一区块的物体被视作一组。而通常情况下，四等分平面便可以解决大部分问题，为了更细致地分组，我们可以进行递归四等分，如下图所示：

具体的划分规则是：每个物体都被划分到各自的 区块（block） 里，同样地，每个最小区块（即没有子区块的区块）最多只能有一个物体，若一个最小区块有了一个以上的物体，则这个区块要被再次四等分。

我们需要对这种理论操作找对对应的数据结构，显然这种结构就是四叉树（quad-tree）。有了上面的划分规则，我们很容易想到如何构建该四叉树：首先，若一个区块是最小区块，那么该区块里面的物体单独为一个四叉树叶子节点，否则，该区块就是一个非叶子节点。这种子区块的划分最终将建造一棵非完全的树。上图平面的划分对应的四叉树如下图所示：

如果要计算某个单独的物体所受到的合力，那么就从根开始遍历树中的节点。如果一个非叶子节点的质心离某个物体足够远，那么就将树中那个部分所包含的物体近似看成一个整体，其位置就是整组物体的质心，其质量就是整组物体的总质量。这个算法相当高效，因为我们无需逐一地单独检验某一组中的个别物体。关于利用四叉树计算N体力的具体方法这里先不细说。

解决N体力问题只是四叉树的一种应用，少量动态物体与大量静态物体的碰撞检测也可以用四叉树分块做粗检查阶段。

四叉树的构建

为了构建一棵四叉树，我们将一个一个地向树中插入节点。具体来说，当向一个由（根）节点Root 所表示的树中插入一个节点A时，我们就递归地执行如下步骤（注意这里我们给“根”字加了一个括号的意思是Root不一定是整棵树的根，它也可能是其中某个子树的根）：

如果Root为空节点（即没有任何物体），则把新的物体A作为Root。
如果Root是一个非叶子节点（即非最小区块），就更新Root的总质量和质心。递归地将物体A插入到Root的子区块中的某个，即A成为Root的四个孩子节点中的某个。
如果Root是一个叶子节点（即最小区块），它包含有一个物体B，那么也就是说在一个最小区块里包含有两个物体Root和B。那么便需要将该区块划分为四个子区块：新建一个非叶子节点C，把Root保存为A，将C作为Root，然后递归地插入物体A和B到Root中。最终更新Root的质心和总质量。

下面这个例子演示了构建一棵包含5个节点的四叉树的过程，我们按A，B，C，D，E的顺序插入节点。

其中，根节点中储存了全部5个物体的质心和总质量。另外一个非叶子节点则包含有B、C和D三个物体所构成的整体的总质量和质心。

四叉树的实现

首先，定义好四叉树的两种节点：叶子节点（最小区块）和非叶子节点（非最小区块）：

// 区块区间
export class Domain {
    // 区块左上角横坐标
    x: number;
    // 区块左上角纵坐标 
    y: number;
    // 区块宽度
    width: number;
    // 区块高度
    height: number;
}

// 子区块
export class ChildBlock {
    // north-west：左上角子区块
    nw: NodeType;
    // north-east：右上角子区块
    ne: NodeType;
    // south-east：右下角子区块
    se: NodeType;
    // south-west：左下角子区块
    sw: NodeType;
}

// 非最小区块节点
export class BlockNode {
    // 四个子区块
    public childBlock: ChildBlock;
    // 该区块的区间
    public domain: Domain;

    // 质量
    public mass: number;
    // 质心
    public centroid: number[];

    constructor(domain: Domain) {
        this.domain = domain;

        // 初始化子区间
        this.childBlock = {
            nw: null,
            ne: null,
            se: null,
            sw: null
        };
    }

    // 重新质点计算
    calcCentroid() {}

    // 重新计算质量
    calcMass() {}
}

// 最小区块节点（即叶节点，物体节点）
export class ItemNode {
    // 物体
    public item: any;
    // 该区块的区间
    public domain: Domain;
    // 质心
    public mass: number;
    // 质心
    public centroid: number[];
    // 受到的合力
    public force: number[];

    constructor(item) {
        this.item = item;

        this.mass = item.mass;
        this.centroid = item.centroid;
        this.force = [0, 0];
    }

    // 重新计算合外力
    calcForce() {}
}

// 节点类型：叶子节点和非叶子节点的联合类型
export type NodeType = BlockNode | ItemNode;

关于质心，质量和合力的计算这里先不实现。

四叉树类实现如下：

// 四叉树
export default class QuadTree {
    // 根节点
    private root: NodeType;
    // 最外层区块的区间
    private domain: Domain;

    constructor(domain: Domain) {
        // 初始化根节点
        this.root = null;
        this.domain = domain;
    }

    /**
     * 计算给定的点在哪个子区块中
     * @param domain 当前区块区间
     * @param pos 给定的点
     * @return {
     *     block: string 位置
     *     childDomain：Domian 子区块区间
     * }
     */
    pos2Block(domain: Domain, pos: number[]) {
        let x = domain.x,
            y = domain.y,
            w = domain.width,
            h = domain.height,
            // 位置（nw，ne，se，sw）
            block = '',
            // 子区块区间
            childDomain: Domain = new Domain();

        childDomain.width = w/2;
        childDomain.height = h/2;

        // 点在上方（north）
        if(pos[1] > y && pos[1] < y + h/2) {
            block += 'n';
            childDomain.y = y;
        }

        // 点在下方（south）
        if(pos[1] > y + h/2 && pos[1] < y + h) {
            block += 's';
            childDomain.y = y + h/2;
        }

        // 点在左方（west）
        if(pos[0] > x && pos[0] < x + w/2) {
            block += 'w';
            childDomain.x = x;
        }

        // 点在右方（east）
        if(pos[0] > x + w/2 && pos[0] < x + w) {
            block += 'e';
            childDomain.x = x + w/2;
        }

        return {
            block,
            childDomain
        };
    }

    /**
     * 插入节点
     * @param node 当前节点
     * @param itemNode 要插入的新节点
     */
    insert(node: NodeType, itemNode: ItemNode): NodeType {
        // 若当前节点为叶子节点
        if(node instanceof ItemNode) {
            // 保存当前叶子节点
            let t = node;

            // 将当前节点划分，新建非叶子节点取代叶子节点
            node = new BlockNode(t.domain);

            // 插入刚保存的节点到新建的非叶子节点
            node = this.insert(node, t);
            // 插入新节点到新建的非叶子节点
            node = this.insert(node, itemNode);
        }
        // 若当前节点为非叶子节点
        else if(node instanceof BlockNode) {
            // 计算新加入的节点该放到哪个子区块
            let blockInfo = this.pos2Block(node.domain, itemNode.centroid);

            // 更新该新节点的区块区间
            itemNode.domain = blockInfo.childDomain;
            // 将新节点插入
            node[blockInfo.block] = this.insert(node[blockInfo.block], itemNode);
        }
        // 若当前节点为空，新节点直接作为当前节点
        else {
            node = itemNode;
        }

        return node;
    }

    /**
     * 添加节点
     * @param item 物体
     */
    addNode(item) {
        // 新建一个叶子节点
        let itemNode = new ItemNode(item);

        // 若当前四叉树为空，则设置该叶子节点的区间为最外层区块的区间
        if(this.root === null) {
            itemNode.domain = this.domain;
        }

        // 插入新节点到四叉树
        this.root = this.insert(this.root, itemNode);
    }
}

最终效果如下，每次新加入物体都会动态划分区块：

---EOF---

phenomLi / Blog