Пересечение плоскостей, ограниченных лучами
Я пытаюсь написать трассировщик лучей в свободное время. В настоящее время пытаются сделать пересечение лучей с ограниченными плоскостями.
Моя программа уже работает с бесконечными плоскостями. Я пытаюсь вычислить математику для небесконечных плоскостей. Пробовал гуглить, но все ресурсы говорят только о бесконечных плоскостях.
Моя плоскость имеет угловую точку (называемую позицией), от которой отходят два вектора (u и v) (их длина соответствует длине сторон). Луч имеет начало и направление.
Сначала я вычисляю точку пересечения с бесконечной плоскостью по формуле
t = normal * (position - origin) / (normal * direction)
Нормаль вычисляется как перекрестное произведение u и v. Тогда по формуле
origin + direction * t
Я получаю саму точку пересечения.
Следующий шаг — проверить, находится ли эта точка в границах прямоугольника, и здесь у меня возникли проблемы.
Моя идея состояла в том, чтобы получить относительный вектор intersection - position, который простирается от угла плоскости до точки пересечения, затем преобразовать его в новый базис u, нормали и v, а затем проверить, короче ли длины преобразованных векторов, чем u и вектора v.
bool BoundedPlane::intersect(const Vec3f &origin, const Vec3f &direction, float &t) const {
t = normal * (position - origin) / (normal * direction);
Vec3f relative = (origin + direction * t) - position;
Mat3f transform{
Vec3f(u.x, normal.x, v.x),
Vec3f(u.y, normal.y, v.y),
Vec3f(u.z, normal.z, v.z)
};
Vec3f local = transform.mul(relative);
return t > 0 && local.x >= 0 && local.x <= u.x && local.z <= 0 && local.z <= v.z;
}
В конце я проверяю, больше ли t, чем 0, что означает, что пересечение находится перед камерой, и находятся ли длины векторов внутри границ. Это дает мне странную строку:
.
Плоскость должна появиться под сферами вот так:
(это использовало ручную проверку, чтобы увидеть, правильно ли они отображаются, если числа верны).
Я не уверен, что я делаю неправильно, и есть ли более простой способ проверить границы. Заранее спасибо.
Редактировать1:
Я переместил вычисления матрицы преобразования в конструктор, так что теперь проверка пересечения выглядит так:
bool BoundedPlane::intersect(const Vec3f &origin, const Vec3f &direction, float &t) const {
if (!InfinitePlane::intersect(origin, direction, t)) {
return false;
}
Vec3f local = transform.mul((origin + direction * t) - position);
return local.x >= 0 && local.x <= 1 && local.z >= 0 && local.z <= 1;
}
Элемент преобразования является инверсией матрицы преобразования.
Добавил картинку в пост. Что вы подразумеваете под случаем отказа? Когда луч вообще не пересекается (в этом случае функция просто возвращает false из-за условия t>0)?
Случай сбоя — это случай, когда функция дает неправильные результаты. То есть ограниченная плоскость и луч, которые пересекаются должен, но пересекаются не надо, или которые пересекаются не должен, но делать.
Я все еще не уверен, что ваш вопрос. Вы просите конкретное дело с цифрами?
Список возможных проблем: 1) Векторы, которые вы передаете в конструктор Mat3f, это строки или столбцы? - если это строки, то они должны быть тремя базисными векторами. 2) Они также должны быть нормализованы, иначе вам нужно будет учитывать масштабирование. 3) Поскольку local находится в... базисе местный плоскости, он должен быть связан (|u|, |v|) (или (width, length), если вы решите нормализовать u, v), а не (u.x, v.z), которые являются координатами в базисе Мир.
Ответы 1
Могу ли я предложить другой подход? Рассмотрим кадр с началом
position и базисные векторы
u = { u.x, u.y, u.z }
v = { v.x, v.y, v.z }
direction = { direction.x, direction.y, direction.z}
Шаг 1: Сформировать матрицу
M = {
{u.x, v.x, direction.x},
{u.y, v.y, direction.y},
{u.z, v.z, direction.z}
}
Шаг 2: Вычислить вектор w, который является решением системы линейных уравнений 3 x 3
M * w = origin - position, т.е.
w = inverse(M) * (origin - position);
Убедитесь, что direction не лежит в одной плоскости с u, v, иначе пересечения нет и inverse(M) не существует.
Шаг 3: если 0.0 <= w.x && w.x <= 1.0 && 0.0 <= w.y && w.y <= 1.0, то прямая пересекает параллелограмм, натянутый векторами u, v, и точка пересечения
w0 = { w.x, w.y , 0 };
intersection = position + M * w0;
иначе прямая не пересекает параллелограмм, натянутый векторами u, v
Идея этого алгоритма состоит в том, чтобы рассмотреть (неортонормированный) фрейм position, u, v, direction. Затем матрица M меняет все в координатах этого нового кадра. В этом кадре линия вертикальна, параллельна оси «z», точка origin имеет координаты w, а вертикальная линия, проходящая через w, пересекает плоскость в точке w0.
Редактировать 1: Вот шаблонная формула для обратной матрицы 3x3:
Если исходная матрица M
a b c
d e f
g h i
инверсия
(1 / det(M)) * {
{e*i - f*h, c*h - b*i, b*f - c*e},
{f*g - d*i, a*i - c*g, c*d - a*f},
{d*h - e*g, b*g - a*h, a*e - b*d},
}
куда
det(M) = a*(e*i - f*h) + b*(f*g - d*i) + c*(d*h - e*h)
является определителем М.
Таким образом, алгоритм обращения может быть следующим:
Дано
M = {
{a, b, c},
{d, e, f},
{g, h, i},
}
- Рассчитать
inv_M = {
{e*i - f*h, c*h - b*i, b*f - c*e},
{f*g - d*i, a*i - c*g, c*d - a*f},
{d*h - e*g, b*g - a*h, a*e - b*d},
};
- Рассчитать
det_M = a*inv_M[1][1] + b*inv_M[2][1] + c*inv_M[3][1];
- Возвратите обратную матрицу M
inv_M = (1/det_M) * inv_M;
Редактировать 2: Давайте попробуем другой подход, чтобы ускорить процесс.
Шаг 1: Для каждой плоскости, определяемой точкой position и двумя векторами u и v, предварительно вычислите следующие величины:
normal = cross(u, v);
u_dot_u = dot(u, u);
u_dot_v = dot(u, v);
v_dot_v = dot(v, v); // all these need to be computed only once for the u and v vectors
det = u_dot_u * v_dot_v - u_dot_v * u_dot_v; // again only once per u and v
Шаг 2: Теперь для заданной линии с точкой origin и направлением direction, как и раньше, вычисляем точку пересечения int_point с плоскостью, натянутой на u и v:
t = dot(normal, position - origin) / dot(normal, direction);
int_point = origin + t * direction;
rhs = int_point - position;
Шаг 3: Расчет
u_dot_rhs = dot(u, rhs);
v_dot_rhs = dot(v, rhs);
w1 = (v_dot_v * u_dot_rhs - u_dot_v * v_dot_rhs) / det;
w2 = (- u_dot_v * u_dot_rhs + u_dot_u * v_dot_rhs) / det;
Шаг 4:
if (0 < = w1 && w1 <= 1 && 0 < = w2 && w2 <= 1 ){
int_point is in the parallelogram;
}
else{
int_point is not in the parallelogram;
}
Итак, что я здесь делаю, так это нахожу точку пересечения линии origin, direction с плоскостью, заданной position, u, v, и ограничиваю себя плоскостью, что позволяет мне работать в 2D, а не в 3D. я представляю
int_point = position + w1 * u + w2 * v;
rhs = int_point - position = w1 * u + w2 * v
и нахождение w1 и w2 путем умножения этого векторного выражения на базисные векторы u и v, что приводит к линейной системе 2x2, которую я решаю напрямую.
Спасибо, в итоге у меня была исходная матрица преобразования, но ваш пост объяснил мне, что моя математика неверна, и мне пришлось инвертировать матрицу. Теперь это работает, но очень медленно по сравнению с тем, что было раньше, простая сцена в вопросе заняла у меня около 2 секунд, чтобы вычислить, с 4 плоскостями это занимает около 20-30 секунд.
@DanielSmith Просто из любопытства, замедляет ли это инверсия матрицы? Какой алгоритм вы используете для инверсии матриц (думаю, это часть библиотеки)?
Да, я думаю, что это инверсия, которая замедляет его. Я использую самописный заголовок линейной алгебры (я пробовал Eigen 3, но это сильно замедлило вычисления). Прямо сейчас это просто жестко закодированная инверсия 3x3, которую я написал, чтобы проверить, дает ли ваш ответ правильный результат. Я собираюсь переделать это.
@DanielSmith Технически вам не нужно инвертировать матрицу, вам просто нужно решить систему линейных уравнений на шаге 2. Может быть, мы можем придумать быстрый решатель для линейной системы вместо инвертирования.
До сих пор я старался избегать каких-либо библиотек, потому что они, как правило, слишком медленны для такого большого количества вычислений, и прошло несколько лет с тех пор, как у меня была линейная алгебра в uni, поэтому обратное звучит проще, чем решатель линейного уравнения. Хотя я открыт для всего.
@DanielSmith Может быть, мы можем просто написать обратную матрицу напрямую, используя формулу шаблона, или мы могли бы использовать алгоритм исключения Гаусса или Гаусса-Джордана: en.wikipedia.org/wiki/Gaussian_elimination
Я только что понял, что, поскольку я использую векторы u, normal и v для матрицы преобразования, мне не нужно вычислять инверсию для каждого теста пересечения, потому что матрица зависит только от плоскости. Я перенес расчеты в конструктор, поэтому проверка пересечения требует лишь нескольких незначительных операций (редактировать 1). Однако это не привело к увеличению времени выполнения.
@DanielSmith Я добавил предложенный алгоритм обращения матрицы. Я должен сказать, что с вашим подходом вы просчитываете много вещей, которые вы могли бы обойти. Вы вычисляете векторное произведение для нормали (которая составляет примерно 2/3 определителя), которая вам на самом деле не нужна, затем вычисляете точку пересечения, а затем вычисляете обратную матрицу 3x3. Именно поэтому я предложил метод выше, потому что я вычисляю точку пересечения и проверку на ее нахождение в параллелограмме, вычисляя только одну обратную матрицу.
@DanielSmith О, и чтобы сэкономить еще немного, чек может быть просто (w.x)*(w.x) <= 1 && (w.y)*(w.y) <= 1
Забыл упомянуть, функция инверсии написана именно так, как вы предложили. Я поищу более быстрые способы инвертирования функций, когда доберусь туда. Что касается других расчетов, то нормаль вычисляется один раз в конструкторе и в любом случае нужна для источников света, так что это не дополнительные вычисления. Я использовал свою версию (u, normal, v) базисных векторов, потому что они зависят только от самой плоскости, поэтому достаточно вычислить инверсию только один раз. Если вы поменяете нормаль на направление луча, вам нужно вычислить инверсию для каждого луча, что должно быть медленнее.
Спасибо за ваш ответ и идеи тоже!
@DanielSmith Я добавил Edit 2, который избегает матриц 3 x 3 и является гораздо более прямым, учитывая уже доступную информацию. Не стесняйтесь попробовать и дайте мне знать, если это быстрее. Надеюсь, я не перепутал математику, но если что-то не так, дайте мне знать.
Я попробовал, и оказалось, что никакого прироста скорости не было (ничего вообще). Это наводит меня на мысль, что замедление происходит где-то еще. Я, вероятно, закончу тестированием каждой части rayracer, чтобы локализовать проблему. Дайте мне знать, если у вас есть какие-либо мысли.
Что такое "странная линия"? У вас есть случай отказа?