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

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

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

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

// 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): Пересечение куба и конуса

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

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

4 Solutions collect form web for “Обнаруживать, пересекаются ли куб и конус друг с другом?”

Для кода

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

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

 /* 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's 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 * PX * dir.X +2 * PY * dir.Y - 2 * (cone.H - PZ) * dir.Z * sqTA; double C = PX * PX + PY * PY - (cone.H - PZ) * (cone.H - PZ) * 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 = PZ + t1 * dir.Z; bool t1_intersect = (t1 >= 0 && t1 <= 1 && z1 >= 0 && z1 <= cone.H); double t2 = (-B - sqrt(delta)) / (2 * A); double z2 = PZ + 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 = PZ + 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.PX * nX + rect.PY * nY + rect.PZ * nZ); // Now, being in the face (ie, coordinates in (u, v) are between 0 and 1) // can be verified through scalar product if(nZ != 0) { Vector3D M(0, 0, -d/nZ); 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.

  1. представить 2 бесконечных линии

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

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

     P+t*(perpendicular vector to cone axis) 

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

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

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

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

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

    • пусть t1 для линии оси P
    • и t2 для линии оси конуса

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

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

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

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

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

    конус

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

[заметки]

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

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

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

    Где:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

БОЛЬШОЙ РЕДАКТОР. Мой предыдущий алгоритм достаточно хорошо справлялся с большинством ситуаций, что я не заметил никаких проблем – это до тех пор, пока я не получу блок времени, достаточный для тщательного тестирования всех ситуаций. В одном случае он имел небольшую погрешность, и в одном случае был довольно большой погрешность при тестировании на 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, но я определил случай сбоя и поэтому выбросил его.

  • Является ли `scipy.misc.comb` быстрее, чем вычисление биномиального ad-hoc?
  • Евклидовой алгоритм (GCD) с несколькими номерами?
  • Интегрирование многомерного интеграла в scipy
  • Теорема Парсеваля в Python
  • Программно программные решения
  • отрицательный pow в python
  • Странное поведение деления в python
  • Как я рисую линейную регрессию
  •  
    Interesting Posts for Van-Lav

    Проксирование на другой веб-сервис с помощью Flask

    Есть ли встроенная функция Python для генерации 100 чисел от 0 до 1?

    Я пытаюсь выбрать структуру для продукта, который я собираюсь построить, и до сих пор я склоняюсь к Нагаре … Любые мысли?

    Разница в цепочке ES6 Promises и PEP3148 Фьючерсы

    Неблокирование чтения на подпроцессе.PIPE в python

    Приложение Flask вызывает ошибку 500 без исключения

    Импорт csv-файла в матрицу / массив в Python

    Threading в приложении PyQt: используйте потоки Qt или потоки Python?

    Возможно ли иметь статические утверждения типа в PyCharm?

    помощь с помощью argparse без дублирования ALLCAPS

    Python / Matplotlib – изменить относительный размер подзаголовка

    Можно ли мариновать объект cookiejar?

    Как создать экспоненциально масштабированную ось?

    Как выполнить проверку подлинности с помощью модуля запросов с использованием хранилища доверенных сертификатов?

    OpenCV – чтение 16-битного изображения в оттенках серого

    Python - лучший язык программирования в мире.