Почему мы должны гомогенизировать эти векторы? - PullRequest
0 голосов
/ 08 мая 2019

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

Я наткнулся на этот пример здесь:

Numpy и пересечения линий

def get_intersect(a1, a2, b1, b2):
    s = np.vstack([a1,a2,b1,b2])        # s for stacked
    h = np.hstack((s, np.ones((4, 1)))) # h for homogeneous
    l1 = np.cross(h[0], h[1])           # get first line
    l2 = np.cross(h[2], h[3])           # get second line
    x, y, z = np.cross(l1, l2)          # point of intersection
    if z == 0:                          # lines are parallel
        return (float('inf'), float('inf'))
    return (x/z, y/z)

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

  1. Почему вектор должен быть однородным (часть, в которой мы заполняем столбец единицами)?
  2. Какгомогенное решение отличается от неоднородного (если оно вообще есть)?
  3. Почему мы проверяем результат только на параллельность вдоль оси Z, а не на X и Y?

Я чувствую, что упускаю что-то очень очевидное, но не могу понять, что это такое.

Ответы [ 2 ]

3 голосов
/ 09 мая 2019

Для справки, уравнение из Википедии:

[]

Пусть a1 = (x1, y1), a2 = (x2, y2), b1 = (x3, y3), b2 = (x4, y4).


Вычислительная интерпретация

Соблюдайтепервые два перекрестных произведения в связанном ответе:

l1 = np.cross(h[0], h[1]) = (x1, y1, 1) ^ (x2, y2, 1)
   = (y1 - y2, x2 - x1, x1*y2 - x2*y1)

l2 = np.cross(h[2], h[3]) = (x3, y3, 1) ^ (x4, y4, 1)
   = (y3 - y4, x4 - x3, x3*y4 - x4*y3)

Эти две строки - все, что требуется для вычисления 6 различных членов в приведенном выше уравнении.И последнее перекрестное произведение:

x, y, z = np.cross(l1, l2)
      x = (x2 - x1) * (x3*y4 - x4*y3) - (x4 - x3) * (x1*y2 - x2*y1)
-->   y = (y3 - y4) * (x1*y2 - x2*y1) - (y1 - x2) * (x3*y4 - x4*y3)
      z = (y1 - y2) * (x4 - y3) - (y3 - y4) * (x2 - x1)

Эти числа в точности совпадают с числителями и знаменателем в уравнении Википедии.

Такое довольно сложное выражение потребовало бы десятков инструкций FPU для вычислениятермин-на-перспектива.

Гомогенизация векторов позволяет использовать этот метод перекрестных произведений, который можно оптимизировать под несколько инструкций SIMD - гораздо более эффективно.


Геометрическая интерпретация

Предположим, вы рассматриваете гомогенизированные векторы как точки в трехмерном пространстве и соединяете каждую пару с началом координат, чтобы сделать два треугольника:

enter image description here

Все 4 точки лежат наплоскость Z = 1 (серая).

Линия L (зеленая) - это пересечение плоскостей двух треугольников (синий + красный), ипроходит через начало координат O и желаемую точку пересечения P (желтый).

нормальный треугольника задается крестом-продукт его боковых векторов.В этом случае боковые векторы просто задаются 4 точками, поскольку другая точка является началом координат.В коде нормали задаются как l1 и l2.

Одно из определений плоскости состоит в том, что все лежащие в ней прямые должны быть перпендикулярны ее нормали.Поскольку прямая L лежит в плоскостях обоих треугольников, она должна быть перпендикулярна l1 и l2, то есть ее направление задается как np.cross(l1, l2).

Гомогенизация позволяетумный последний шаг, который использует подобные треугольники для вычисления P :

enter image description here

if z == 0:                          # lines are parallel
    return (float('inf'), float('inf'))
return (x/z, y/z)                   # Px = x / z, Py = y / z
1 голос
/ 14 мая 2019

Как небольшая поправка к ответам, приведенным выше, использование трехмерных перекрестных произведений для вычисления пересечения двухмерных линий в значительной степени наименее эффективно. Перекрестные продукты не оптимизируют SIMD вообще (если только вы не пойдете на все, чтобы использовать структуру данных SOA, но даже если это чрезвычайно расточительный подход). Наилучший подход - использовать старые добрые уравнения одновременности.

Начиная с наших точек ввода, которые определяют линию:

line1: [(x1, y1), (x2, y2)]
line2: [(x3, y3), (x4, y4)]

Вычислить некоторые векторы направления:

// 1st direction
u1 = x2 - x1
v1 = y2 - y1
D1 = [u1, v1]

// 2nd direction
u2 = x4 - x3
v2 = y4 - y3
D2 = [u2, v2]

Теперь давайте переформулируем уравнения линий в виде лучей и составим уравнение для любой точки вдоль этих лучей:

// coords of a point 'd1' distance along the 1st ray
Px = x1 + d1*u1
Py = y1 + d1*v1

// coords of a point 'd2' distance along the 2nd ray
Px = x3 + d2*u2
Py = y3 + d2*v2

Давайте предположим, что линии пересекаются, что означает, что оба P будут одинаковыми, что позволяет нам утверждать:

x1 + d1*u1 = x3 + d2*u2
y1 + d1*v1 = y3 + d2*v2

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

 d1 = x3 + d2*u2 - x1
      ---------------
            u1

 d1 = y3 + d2*v2 - y1
      ---------------
            v1

Теперь у нас есть два уравнения для d1, поэтому давайте создадим другое уравнение для одновременного получения значения для d2:

x3 + d2*u2 - x1     y3 + d2*v2 - y1
---------------  =  ---------------
      u1                  v1

Переставить для изоляции d2:

 d2 = u1(y3 - y1) - v1(x3 - x1)
      -------------------------
           v1*u2 - u1*v2

Если (v1 * u2 - u1 * v2) окажется равным нулю, для этого уравнения нет решения (назовем его детерминантом, потому что это то, что есть!). Если определитель не равен нулю, то просто вычислите d2, используя приведенное выше уравнение, и вернитесь к одному из наших предыдущих уравнений, чтобы найти значение точки:

Px = x3 + d2*u2
Py = y3 + d2*v2 

Некоторый непроверенный код C ++:

bool computeIntersectPoint(
  float x1, float y1,
  float x2, float y2,
  float x3, float y3,
  float x4, float y4,
  float outPoint[2])
{
  // compute direction vectors
  // in some cases, it can be worth 
  // storing the lines as rays as an 
  // optimisation. (avoids 4 subs)
  const float u1 = x2 - x1;
  const float v1 = y2 - x1;
  const float u2 = x4 - x3;
  const float v2 = y4 - x3;

  // check to see if we have a solution
  // 1 mul, 1 fmsub
  float determinant = v1*u2 - u1*v1; 
  if(determinant == 0)
    return false;

  // 2 sub, 1 mul, 1 fmsub
  float temp = u1 * (y3 - y1) - v1 * (x3 - x1);

  // 1 div
  float intersectDistance = temp / determinant;

  // 2 fma
  outPoint[0] = intersectDistance * u2 + x3;
  outPoint[1] = intersectDistance * v2 + y3;

  return true;
}

Быстрое доказательство десмоса: https://www.desmos.com/calculator/gtlmnmzn6l

Здесь стоит сравнить количество команд между двумя подходами здесь. Для трехмерного перекрестного произведения требуется 3 мульт и 3 инструкции fmsub (или 6 мул + 3 суб, если FMA недоступен). Так как у нас есть 3 из них, мы до: 9 mul и 9 fmsub ops. Добавьте в 2 подразделения, и мы получим:

9 mul
9 fmsub
2 div

Подход, который я разместил, потребует:

1 div
6 sub
4 fma
2 mul

Хотя вы можете сохранить 4 из этих подводных лодок, если вам удастся сохранить линии в виде лучей.

...