如何在一组球体上最有效地执行碰撞检测

假设我有一个有几个内核的CPU，我想找到哪些球体是接触的。任何一组球体，其中每个球体相连。它们都至少接触集合中的一个球体)被称为"群"，并被组织成一个向量，在下面的例子中称为"group_members"。为了实现这一点，我目前使用了一个相当昂贵的操作，在概念上看起来像这样:

vector<Sphere*> unallocated_spheres = all_spheres; // start with a copy of all spheres
vector<vector<Sphere*>> group_sequence; // groups will be collected here
while (unallocated_spheres.size() > 0U) // each iteration of this will represent the creation of a new group
{
    std::vector<Sphere*> group_members; // this will store all members of the current group
    group_members.push_back(unallocated_spheres.back()); // start with the last sphere (pop_back requires less resources than erase)
    unallocated_spheres.pop_back(); // it has been allocated to a group so remove it from the unallocated list
    // compare each sphere in the new group to every other sphere, and continue to do so until no more spheres are added to the current group
    for (size_t i = 0U; i != group_members.size(); ++i) // iterators would be unsuitable in this case
    {
        Sphere const * const sphere = group_members[i]; // the sphere to which all others will be compared to to check if they should be added to the group
        auto it = unallocated_spheres.begin();
        while (it != unallocated_spheres.end())
        {
            // check if the iterator sphere belongs to the same group
            if ((*it)->IsTouching(sphere))
            {
                // it does belong to the same group; add it and remove it from the unallocated_spheres vector and repair iterators
                group_members.push_back(*it);
                it = unallocated_spheres.erase(it); // repair the iterator
            }
            else ++it; // if no others were found, increment iterator manually
        }
    }
    group_sequence.push_back(group_members);
}

有没有人对提高这段代码在wall time方面的效率有什么建议?我的程序花费了相当一部分时间来运行这些循环，任何关于如何在结构上改变它以使其更有效的建议都将不胜感激。

请注意，由于这些是球体，"IsTouching()"是一个非常快速的浮点操作(比较两个球体的位置和半径)。它看起来像这样(注意x,y和z是球体在欧几里得维中的位置):

// input whether this cell is touching the input cell (or if they are the same cell; both return true)
bool const Sphere::IsTouching(Sphere const * const that) const
{
    // Apply pythagoras' theorem in 3 dimensions
    double const dx = this->x - that->x;
    double const dy = this->y - that->y;
    double const dz = this->z - that->z;
    // get the sum of the radii of the two cells
    double const rad_sum = this->radius + that->radius;
    // to avoid taking the square root to get actual distances, we instead compare 
    // the square of the pythagorean distance with the square of the radii sum
    return dx*dx + dy*dy + dz*dz < rad_sum*rad_sum;
}