c++ - 使用 Barnes-Hut 进行图形放置的优化问题

Question

我一直在尝试解决我的图形可视化应用程序中的力导向图/Barnes-Hut 问题。到目前为止，我已经检查了八叉树的创建，它看起来正确(树由方框表示，圆圈是我的图形节点): 我的 Quadtree 中的字段如下:

class Quadtree
{
    public:
        int level;
        Quadtree* trees[2][2][2];
        glm::vec3 vBoundriesBox[8];
        glm::vec3 center;
        bool leaf;
        float combined_weight = 0;
        std::vector<Element*> objects;
        //Addition methods/fields
    private:
    //Additional methods/fields
    protected:
}

这就是我向四叉树递归添加元素的方式:

#define MAX_LEVELS 5

void Quadtree::AddObject(Element* object)
{
    this->objects.push_back(object);
}

void Quadtree::Update()
{
    if(this->objects.size()<=1 || level > MAX_LEVELS)
    {
        for(Element* Element:this->objects)
        {
            Element->parent_group = this;
            this->combined_weight += Element->weight;
        }
        return;
    }

    if(leaf)
    {
        GenerateChildren();
        leaf = false;
    }

    while (!this->objects.empty())
    {
        Element* obj = this->objects.back();
        this->objects.pop_back();
        if(contains(trees[0][0][0],obj))
        {
            trees[0][0][0]->AddObject(obj);
            trees[0][0][0]->combined_weight += obj->weight;
        } else if(contains(trees[0][0][1],obj))
        {
            trees[0][0][1]->AddObject(obj);
            trees[0][0][1]->combined_weight += obj->weight;
        } else if(contains(trees[0][1][0],obj))
        {
            trees[0][1][0]->AddObject(obj);
            trees[0][1][0]->combined_weight += obj->weight;
        } else if(contains(trees[0][1][1],obj))
        {
            trees[0][1][1]->AddObject(obj);
            trees[0][1][1]->combined_weight += obj->weight;
        } else if(contains(trees[1][0][0],obj))
        {
            trees[1][0][0]->AddObject(obj);
            trees[1][0][0]->combined_weight += obj->weight;
        } else if(contains(trees[1][0][1],obj))
        {
            trees[1][0][1]->AddObject(obj);
            trees[1][0][1]->combined_weight += obj->weight;
        } else if(contains(trees[1][1][0],obj))
        {
            trees[1][1][0]->AddObject(obj);
            trees[1][1][0]->combined_weight += obj->weight;
        } else if(contains(trees[1][1][1],obj))
        {
            trees[1][1][1]->AddObject(obj);
            trees[1][1][1]->combined_weight += obj->weight;
        }
    }

    for(int i=0;i<2;i++)
    {
        for(int j=0;j<2;j++)
        {
            for(int k=0;k<2;k++)
            {
                trees[i][j][k]->Update();
            }
        }
    }
}

bool Quadtree::contains(Quadtree* child, Element* object)
{
    if(object->pos[0] >= child->vBoundriesBox[0][0] && object->pos[0] <= child->vBoundriesBox[1][0] &&
       object->pos[1] >= child->vBoundriesBox[4][1] && object->pos[1] <= child->vBoundriesBox[0][1] &&
       object->pos[2] >= child->vBoundriesBox[3][2] && object->pos[2] <= child->vBoundriesBox[0][2])
        return true;
    return false;
}

正如您在图片中看到的那样，节点非常密集。我一直在尝试找出解决排斥力计算的方法，但它仍然不起作用，结果还是一样。

那么我是如何计算它的:

首先在我的主文件中，我在所有图形节点中运行循环:

for(auto& n_el:graph->node_vector)
{
    tree->CheckNode(&n_el);
}

接下来在我的Qyadtree类中，(tree是这个类对象)，我有这个递归方法:

void Quadtree::CheckNode(Node* node)
{
    glm::vec3 diff = this->center - node->pos;

    double distance_sqr = (diff.x * diff.x) + (diff.y*diff.y) + (diff.z*diff.z);
    double width_sqr = (vBoundriesBox[1][0] - vBoundriesBox[0][0]) * (vBoundriesBox[1][0] - vBoundriesBox[0][0]);
    if(width_sqr/distance_sqr < 10.0f || leaf)
    {
        if(leaf)
        {
            for(auto& n: objects)
            {
                n->Repulse(&objects);
            }
        }
        else
        {
            node->RepulseWithGroup(this);
        }
    }
    else
    {
        for(int i=0; i<2; i++)
        {
            for(int j=0; j<2; j++)
            {
                for(int k=0; k<2; k++)
                {
                    trees[i][j][k]->CheckNode(node);
                }
            }
        }
    }
}

最后，我有两种方法计算排斥力，具体取决于排斥力是在组和节点之间还是在两个节点之间:

double Node::Repulse(std::vector<Node*>* nodes)
{
    double dx;
    double dy;
    double dz;
    double force = 0.0;
    double distance_between;
    double delta_weights;
    double temp;
    for(auto& element_node:*nodes)
    {
        if(this->name == element_node->name)
        {
            continue;
        }
        if(!element_node->use) continue;
        delta_weights = 0.5 + abs(this->weight - element_node->weight);
        dx = this->pos[0] - element_node->pos[0];
        dy = this->pos[1] - element_node->pos[1];
        dz = this->pos[2] - element_node->pos[2];
        distance_between = dx * dx + dy * dy + dz * dz;
        force = 0.19998 * delta_weights/(distance_between * distance_between);
        temp = std::min(1.0, force);
        if(temp<0.0001)
        {
            temp = 0;
        }
        double mx = temp * dx;
        double my = temp * dy;
        double mz = temp * dz;
        this->pos[0] += mx;
        this->pos[1] += my;
        this->pos[2] += mz;
        element_node->pos[0] -= mx;
        element_node->pos[1] -= my;
        element_node->pos[2] -= mz;
    }
}

void Node::RepulseWithGroup(Quadtree* tree)
{
    double dx;
    double dy;
    double dz;
    double force = 0.0;
    double distance_between;
    double delta_weights;
    double temp;

    delta_weights = 0.5 + abs(this->weight - tree->combined_weight);
    dx = this->pos[0] - tree->center.x;
    dy = this->pos[1] - tree->center.y;
    dz = this->pos[2] - tree->center.z;
    distance_between = dx * dx + dy * dy + dz * dz;
    force = 0.19998 * delta_weights/(distance_between * distance_between);
    temp = std::min(1.0, force);
    if(temp<0.0001)
    {
        temp = 0;
    }
    double mx = temp * dx;
    double my = temp * dy;
    double mz = temp * dz;
    this->pos[0] += mx + this->parent_group->repulsion_force.x;
    this->pos[1] += my + this->parent_group->repulsion_force.y;
    this->pos[2] += mz + this->parent_group->repulsion_force.z;
}

如果这个想法:

if(width_sqr/distance_sqr < 10.0f || leaf)
    {
        if(leaf)
        {
            for(auto& n: objects)
            {
                n->Repulse(&objects);
            }
        }
        else
        {
            node->RepulseWithGroup(this);
        }
    }

不清楚是因为我发现，一片树叶中实际上可能有多个元素。如果可能已经达到最大级别并且静止元素位于一个框中，则可能会发生这种情况。然后我还需要计算盒子内对内部节点的力。

更让我困扰的是这种方法的速度(这表明八叉树无法正常工作)是速度。这是表示时间/节点数的简单图: 据我所知，原始力导向图算法的复杂度为 O(n^2)，但对于 Barnes-Hut，它应该是 O(nlogn)。然而，情节甚至与 nlogn 相去甚远。

有人可以告诉我我在这里做错了什么吗？我已经查看这段代码很长时间了，但我看不出我遗漏了什么。

编辑:

根据@Ilmari Karonen 的回答，我对 MAX_LEVELS 5、20、50、100 进行了测试。结果如下。看起来没有什么有意义的区别我会说(不幸的是)

Answer 1

就在我的头顶，

#define MAX_LEVELS 5

看起来很低。您可能只是用尽了八叉树的深度，导致您的算法退化为 O(n²) 直接求和。您可能想尝试将 MAX_LEVELS 增加到一个明显更高的值(至少，比如 10 或 20)，看看这是否会提高性能。

我没有测试过你的代码，所以我不能确定这是真正的问题还是唯一的问题。但这绝对是我首先要检查的内容。

---

仔细查看您的代码，我还发现了其他几个潜在问题。严格来说，这些可能不会影响性能，但它们可能会影响结果的正确性。

首先，您的 Quadtree 类中有一个 center vector ，大概代表子树中节点的质心，但您似乎从未更新该 vector 将节点添加到树中时。由于您确实在计算中使用该 vector ，因此您可能会得到虚假的结果。

(事实上，由于您使用 center vector 的一件事是计算节点和子树之间的距离，因此决定是否深入子树，这也可能弄乱你的表现。)

此外，您似乎在遍历树时直接更新位置，这意味着您的算法生成的轨迹将取决于遍历节点和扩展树的顺序。为了获得更一致和可重现的结果，您可能需要首先计算算法当前迭代期间每个节点的位移，将其存储在单独的 vector 中，然后在节点上运行第二次传递以将位移添加到它们的位置(并为下一次迭代重置它)。

此外，我肯定不会是唯一一个发现您有一个名为 Quadtree 的类实现了一个 oc 树这一事实的人，我可以吗？ :)