Подтвердить что ты не робот

Обнаруживать, пересекаются ли куб и конус друг с другом?

Рассмотрим два геометрических объекта в 3D:

  • куб, выровненный с осями и определяемый положением его центра и его протяженностью (длина края)
  • конус, не выровненный с осями и определяемый положением его вершины, расположением центра его основания и полууголом в вершине

Вот небольшой код для определения этих объектов в С++:

// Preprocessor
#include <iostream>
#include <cmath>
#include <array>

// 3D cube from the position of its center and the side extent
class cube
{ 
    public:
        cube(const std::array<double, 3>& pos, const double ext)
        : _position(pos), _extent(ext) 
        {;}
        double center(const unsigned int idim) 
            {return _position[idim];}
        double min(const unsigned int idim)
            {return _position[idim]-_extent/2;}
        double max(const unsigned int idim)
            {return _position[idim]+_extent/2;}
        double extent()
            {return _extent;}
        double volume()
            {return std::pow(_extent, 3);}
    protected:
        std::array<double, 3> _position;
        double _extent;
};

// 3d cone from the position of its vertex, the base center, and the angle
class cone
{
    public:
        cone(const std::array<double, 3>& vert, 
             const std::array<double, 3>& bas, 
             const double ang)
        : _vertex(vert), _base(bas), _angle(ang)
        {;}
        double vertex(const unsigned int idim)
            {return _vertex[idim];}
        double base(const unsigned int idim)
            {return _base[idim];}
        double angle()
            {return _angle;}
        double height()
            {return std::sqrt(std::pow(_vertex[0]-_base[0], 2)+std::pow(
            _vertex[1]-_base[1], 2)+std::pow(_vertex[2]-_base[2], 2));}
        double radius()
            {return std::tan(_angle)*height();}
        double circle()
            {return 4*std::atan(1)*std::pow(radius(), 2);}
        double volume()
            {return circle()*height()/3;}
    protected:
        std::array<double, 3> _vertex;
        std::array<double, 3> _base;
        double _angle;
};

Я хотел бы написать функцию, чтобы определить, пусто или нет пересечение куба и конуса:

// Detect whether the intersection between a 3d cube and a 3d cone is not null
bool intersection(const cube& x, const cone& y)
{
    // Function that returns false if the intersection of x and y is empty
    // and true otherwise
}

Вот иллюстрация проблемы (иллюстрация в 2D, но моя проблема в 3D): Cube and cone intersection

Как это сделать эффективно (я ищу алгоритм, поэтому ответ может быть в C, С++ или Python)?

Примечание. Здесь пересечение определяется как: он содержит ненулевой объем 3D, который находится в кубе и в конусе (если куб находится внутри конуса, или если конус находится внутри куба, они пересекаются).

4b9b3361

Ответ 1

Для кода

Этот ответ будет немного более общим, чем ваша проблема (например, вместо куба я использую поле). Адаптация к вашему делу должна быть очень простой.

Некоторые определения

/*
    Here is the cone in cone space:

            +        ^
           /|\       |
          /*| \      | H
         /  |  \     |
        /       \    |
       +---------+   v

    * = alpha (angle from edge to axis)
*/
struct Cone // In cone space (important)
{
    double H;
    double alpha;
};

/*
    A 3d plane
      v
      ^----------+
      |          |
      |          |
      +----------> u
      P
*/
struct Plane
{
    double u;
    double v;
    Vector3D P;
};

// Now, a box.
// It is assumed that the values are coherent (that only for this answer).
// On each plane, the coordinates are between 0 and 1 to be inside the face.
struct Box
{
    Plane faces[6];
};

Пересечение линии - конус

Теперь давайте вычислим пересечение между отрезком и нашим конусом. Обратите внимание, что я буду делать вычисления в пространстве конуса. Также обратите внимание, что я беру ось Z как вертикальную. Изменение его на Y остается как упражнение для читателя. Предполагается, что линия находится в пространстве конуса. Направление сегмента не нормируется; вместо этого отрезок имеет длину вектора направления и начинается в точке P:

/*
    The segment is points M where PM = P + t * dir, and 0 <= t <= 1
    For the cone, we have 0 <= Z <= cone.H
*/
bool intersect(Cone cone, Vector3D dir, Vector3D P)
{
    // Beware, indigest formulaes !
    double sqTA = tan(cone.alpha) * tan(cone.alpha);
    double A = dir.X * dir.X + dir.Y * dir.Y - dir.Z * dir.Z * sqTA;
    double B = 2 * P.X * dir.X +2 * P.Y * dir.Y - 2 * (cone.H - P.Z) * dir.Z * sqTA;
    double C = P.X * P.X + P.Y * P.Y - (cone.H - P.Z) * (cone.H - P.Z) * sqTA;

    // Now, we solve the polynom At² + Bt + C = 0
    double delta = B * B - 4 * A * C;
    if(delta < 0)
        return false; // No intersection between the cone and the line
    else if(A != 0)
    {
        // Check the two solutions (there might be only one, but that does not change a lot of things)
        double t1 = (-B + sqrt(delta)) / (2 * A);
        double z1 = P.Z + t1 * dir.Z;
        bool t1_intersect = (t1 >= 0 && t1 <= 1 && z1 >= 0 && z1 <= cone.H);

        double t2 = (-B - sqrt(delta)) / (2 * A);
        double z2 = P.Z + t2 * dir.Z;
        bool t2_intersect = (t2 >= 0 && t2 <= 1 && z2 >= 0 && z2 <= cone.H);

        return t1_intersect || t2_intersect;
    }
    else if(B != 0)
    {
        double t = -C / B;
        double z = P.Z + t * dir.Z;
        return t >= 0 && t <= 1 && z >= 0 && z <= cone.H;
    }
    else return C == 0;
}

Прямоугольное пересечение

Теперь мы можем проверить, пересекает ли прямоугольная часть плана конус (это будет использоваться для проверки того, пересекает ли кусок куба конус). Все еще в пространстве конуса. План передается таким образом, который поможет нам: 2 вектора и точка. Для упрощения вычислений векторы не нормируются.

/*
    A point M in the plan 'rect' is defined by:
        M = rect.P + a * rect.u + b * rect.v, where (a, b) are in [0;1]²
*/
bool intersect(Cone cone, Plane rect)
{
    bool intersection = intersect(cone, rect.u, rect.P)
                     || intersect(cone, rect.u, rect.P + rect.v)
                     || intersect(cone, rect.v, rect.P)
                     || intersect(cone, rect.v, rect.P + rect.u);

    if(!intersection)
    {
        // It is possible that either the part of the plan lie
        // entirely in the cone, or the inverse. We need to check.
        Vector3D center = P + (u + v) / 2;

        // Is the face inside the cone (<=> center is inside the cone) ?
        if(center.Z >= 0 && center.Z <= cone.H)
        {
            double r = (H - center.Z) * tan(cone.alpha);
            if(center.X * center.X + center.Y * center.Y <= r)
                intersection = true;
        }

        // Is the cone inside the face (this one is more tricky) ?
        // It can be resolved by finding whether the axis of the cone crosses the face.
        // First, find the plane coefficient (descartes equation)
        Vector3D n = rect.u.crossProduct(rect.v);
        double d = -(rect.P.X * n.X + rect.P.Y * n.Y + rect.P.Z * n.Z);

        // Now, being in the face (ie, coordinates in (u, v) are between 0 and 1)
        // can be verified through scalar product
        if(n.Z != 0)
        {
            Vector3D M(0, 0, -d/n.Z);
            Vector3D MP = M - rect.P;
            if(MP.scalar(rect.u) >= 0
               || MP.scalar(rect.u) <= 1
               || MP.scalar(rect.v) >= 0
               || MP.scalar(rect.v) <= 1)
                intersection = true;
        }
    }
    return intersection;
}

Коробчатое пересечение коней

Теперь, окончательная часть: весь куб:

bool intersect(Cone cone, Box box)
{
    return intersect(cone, box.faces[0])
        || intersect(cone, box.faces[1])
        || intersect(cone, box.faces[2])
        || intersect(cone, box.faces[3])
        || intersect(cone, box.faces[4])
        || intersect(cone, box.faces[5]);
}

Для математики

Все еще в пространстве конуса уравнения конуса:

// 0 is the base, the vertex is at z = H
x² + y² = (H - z)² * tan²(alpha)
0 <= z <= H

Теперь параметрическое уравнение линии в 3D:

x = u + at
y = v + bt
z = w + ct

Вектор направления (a, b, c), а точка (u, v, w) лежит на прямой.

Теперь давайте сведем уравнения вместе:

(u + at)² + (v + bt)² = (H - w - ct)² * tan²(alpha)

Затем, после разра- ботки и рефакторизации этого уравнения, получим следующее:

At² + Bt + C = 0

где A, B и C показаны в первой функции пересечения. Просто разрешите это и проверьте граничные условия на z и t.

Ответ 2

  • представить 2 бесконечные строки

    • ось конуса
    • проходящей через точку P (центр куба для стартеров), которая перпендикулярна оси конуса.

    Вам известна коническая ось, так что это легко, вторая строка определяется как

    P+t*(perpendicular vector to cone axis)
    

    Этот вектор может быть получен поперечным произведением вектора оси конуса и вектора, перпендикулярного вашему изображению (с учетом оси Z). t - это параметр скалярного значения...

  • вычислить пересечение этих двух линий/осей

    если вы не знаете, какие уравнения вывести их или google. Пусть точка пересечения Q

  • , если точка пересечения Q не лежит внутри конуса

    (между вершиной и базой), то точка P не пересекает конус. Из уравнений пересечения вы получите параметры t1 и t2

    • let t1 для оси P оси
    • и t2 для линии оси конуса

    если вектор направления оси линии также является длиной конуса, тогда пересечение внутри конуса, если t2 = <0,1>

  • , если P не находится внутри треугольника (разрезанный конус к плоскости, порожденной этими двумя осями)

    вам также легко знать положение внутреннего конуса Q (t2), чтобы вы знали, что конус находится в P -axis от Q до расстояния R*t2, где R радиус основания конуса. Таким образом, вы можете вычислить |P-Q| и проверить, есть ли это <=R*t2 или использовать напрямую t1 (если P вектор направления оси - единица).

    если расстояние больше, чем R*t2 точка P не пересекает конус.

  • , если # 3 и # 4 положительны, то P пересекает конус

    cone

    • надеюсь, что вы не возражаете, вот ваш образ с несколькими вещами, добавленными для ясности.

[примечания]

Теперь в жесткой части есть граничные случаи, когда вершина куба не пересекает конус, но сам куб все равно пересекает конус. Это может произойти, когда ||P-Q|-R*t2| = <0,half cube size> В этом случае вы должны проверить больше точек, а затем просто кубировать вершины вдоль ближайшей грани куба.

Другой подход:

  • создать матрицу преобразования для конуса

    Где:

    • его вершина как начало
    • его ось как ось +Z
    • и XY плоскость параллельна своей базе

    поэтому любая точка внутри конуса, если

    • Z = <0,h>
    • X*X + Y*Y <= (R*Z/h)^2 или X*X + Y*Y <= (R*Z*tan(angle))^2
  • конвертировать вершины куба в пространство конуса

    и проверьте, не верна ли какая-либо вершина внутри конуса, вы также можете проверить все граничные линии куба с условиями из # 1 (алгебраически) или использовать больше точек вдоль кубических граней, как в предыдущем методе.

Обсуждение чата: http://chat.stackoverflow.com/rooms/48756/discussion-between-spektre-and-joojaa

Ответ 3

Информация: Я не знаю, является ли эта идея уже запатентованной интеллектуальной собственностью (в вашем регионе), или нет, или как узнать, или что-то еще. Я делаю это для удовольствия.:)

Но вот говядина:

  • Шаг 1: Аппроксимация: Для эффективности обрабатывайте оба объекта как сферы (используйте внешние сферы). Вычислите их distance (между их двумя центральными точками), чтобы выяснить, достаточно ли они, чтобы пересечься по все. Быстрое возвращение false, если они не могут пересекаться (потому что их расстояние больше суммы радиуса обеих сфер).

  • Шаг 2: Точные вычисления: Это простой способ сделать это: интерпретировать конус как пакет трехмерных пикселей, называемый voxels (или legos): выберите любое разрешение (гранулярность), которое вы считаете приемлемым (возможно, 0,01). Создайте вектор, указывающий от (0,0,0) до любой воксельной точки внутри тома вашего конуса (начиная с точки, которую вы уже назвали "вершиной" ). Верно true, если координата этого воксела существует внутри вашего куба. Повторите это для каждого воксела, который вы можете вычислить для вашего объекта конуса, на основе выбранной гранулярности.

  • Шаг 3: Если ничего не соответствует, верните false.

Оптимизация. Функция, которая определяет, может ли какая-либо заданная трехмерная точка внутри куба быть оптимизирована, рассмотрим внутреннюю сферу куба. Смысл, любая данная 3-D точка находится внутри куба, если она достаточно близко, чтобы центр внутренней сферы куба находился внутри этой сферы. Веселье начинается, когда вы начинаете заполнять пустые углы куба дополнительными сферами для большей оптимизации (это совершенно необязательно).

Я уверен, что для шага 2 возможны дальнейшие оптимизации. Однако приятная вещь об этом подходе заключается в том, что вы можете свободно регулировать детализацию, точную настройку между временем вычисления и точностью вычислений.

Вы также можете создать решатель, который автоматически уменьшает детализацию по нескольким итерациям. Значение точности увеличивается с течением времени (для лучшего разрешения).

Ответ 4

Несколько лет назад я создал эффективный алгоритм для проверки пересечения между конусом и aabb для некоторого кода рендеринга. Я недавно нуждался в этом тесте для того, над чем я сейчас работаю, поэтому я пересмотрел его и сделал его еще более эффективным. Я потратил больше времени, чем я хотел бы признать работу над этой проблемой, и поэтому решил избавить вас от страдания и опубликовать код. Таким образом, это AFAIK полностью уникальное решение и не будет найдено в учебнике (т.е. Решение Дэвида Эберли).

BIG EDIT: мой предыдущий алгоритм обрабатывал большинство ситуаций достаточно хорошо, что я не заметил никаких проблем - это до тех пор, пока я не получу блок времени, достаточный для тщательного тестирования всех ситуаций. В одном случае он имел небольшую погрешность, и в одном случае был довольно большой погрешность при тестировании на 6-плоскостях грубой силы. Тест на сегмент-ящик, который я делал в конце, не мог объяснить все возможные странные ситуации, когда расширяющийся конус проникал в ящик. Это выглядело хорошо на моих 2-х и плохо-нарисованных 3D-тестах, но не на практике. Моя новая идея несколько дороже, но она пуленепробиваемая.

Шаги, которые выполняет новая версия, следующие:

1.) Раньше, если вершина находится в ограничивающей рамке.

2.) Определите лица, к которым может коснуться конус (Макс. 3), и добавьте их вершины в массив.

3.) Проектировать все вершины в массиве в "пространство конуса".

4.) Раньше, если проецируемые вершины находятся внутри конуса

5.) Проведите круг-полигон-тест против массива вершин и радиус конуса, чтобы поймать пересечения ребер с полигоном на границе конуса

bool Intersect(const Cone& pCone) const {
    Vector3 pFaceVerts[12];
    U32 uVertCount;
    int piClipSigns[3];
    U32 uClipCount = GetClipInfo(pCone.GetApex(), piClipSigns);

    switch (uClipCount) {
        // If the clip count is zero, the apex is fully contained in the box
        xcase 0: {
            return true;
        }
        // 1) Clips single face, 4 vertices, guaranteed to not touch any other faces
        xcase 1: {
            int iFacet = piClipSigns[0] != 0 ? 0 : (piClipSigns[1] != 0 ? 1 : 2);
            GetFacetVertices(iFacet, piClipSigns[iFacet], pFaceVerts);
            uVertCount = 4;
        }
        // 2) Clips an edge joining two candidate faces, 6 vertices
        // 3) Clips a vertex joining three candidate faces, 7 vertices
        xcase 2:
        acase 3: {
            uVertCount = 0;
            for (U32 iFacet = 0; iFacet < 3; iFacet++) {
                if (piClipSigns[iFacet] != 0) {
                    GetFacetVertices(iFacet, piClipSigns[iFacet], pFaceVerts + uVertCount);
                    uVertCount += 4;
                }
            }
            FixVertices(pFaceVerts, uVertCount);
        }
    }

    // Project vertices into cone-space
    F32 fConeRadiusSquared = Square(pCone.GetRadius());
    F32 pfLengthAlongAxis[6];
    bool bOutside = true;
    for (U32 i = 0; i < uVertCount; i++) {
        pfLengthAlongAxis[i] = Dot(pCone.GetAxis(), pFaceVerts[i] - pCone.GetApex());
        bOutside &= Clamp1(pfLengthAlongAxis[i], LargeEpsilon, pCone.GetHeight() - LargeEpsilon);
    }
    // Outside the cone axis length-wise
    if (bOutside) {
        return false;
    }
    for (U32 i = 0; i < uVertCount; i++) {
        Vector3 vPosOnAxis = pCone.GetApex() + pCone.GetAxis() * pfLengthAlongAxis[i];
        Vector3 vDirFromAxis = pFaceVerts[i] - vPosOnAxis;

        F32 fScale = (pCone.GetHeight() / pfLengthAlongAxis[i]);
        F32 x = fScale * Dot(vDirFromAxis, pCone.GetBaseRight());
        F32 y = fScale * Dot(vDirFromAxis, pCone.GetBaseUp());
        // Intersects if any projected points are inside the cone
        if (Square(x) + Square(y) <= fConeRadiusSquared) {
            return true;
        }
        pFaceVerts[i] = Vector2(x, y);
    }
    // Finally do a polygon circle intersection with circle center at origin
    return PolygonCircleIntersect(pFaceVerts, uVertCount, pCone.GetRadius());
}

GetClipInfo:

inline U32 GetClipInfo(const Vector3& P, int piClipSigns[3]) const {
    U32 N = 0;
    for (U32 i = 0; i < 3; i++) {
        if (P[i] < m_vMin[i]) {
            piClipSigns[i] = -1;
            N++;
        } else if (P[i] > m_vMax[i]) {
            piClipSigns[i] = +1;
            N++;
        } else {
            piClipSigns[i] = 0;
        }
    }
    return N;
}

GetFacetVertices и FixVertices немного взломаны на данный момент, но они получают вершины на лице и фиксируют вершины как выпуклые и ccw-упорядоченные соответственно.

Альтернативой было бы просто проецировать все вершины в коническое пространство без какой-либо логики, но мне нужно, чтобы я был как можно быстрее, поэтому я разбил его на несколько случаев.

Я пробовал различные другие подходы, в частности испытание разделительной оси между осью конуса и коробкой, и использовал самую большую разделительную ось, чтобы получить самую близкую поверхность для проверки с помощью PolygonCircleIntersect, но я определил случай сбоя и поэтому выбросил его.