Хорошо, я пришел сюда в поисках ответа и не нашел ничего простого и понятного, поэтому я пошел дальше и сделал глупую, но эффективную (и относительно простую) вещь: оптимизацию Монте-Карло.
Проще говоря, алгоритм выглядит следующим образом: случайным образом возмущает вашу матрицу проекции, пока она не проецирует ваши известные трехмерные координаты на ваши известные двухмерные координаты.
Вот еще фото от двигателя танка Томаса:
Допустим, мы используем GIMP, чтобы найти двухмерные координаты того, что мы считаем квадратом на плоскости земли (действительно ли это квадрат, зависит от вашего суждения о глубине):
Я получаю четыре очка на 2D-изображении: (318, 247), (326, 312), (418, 241) и (452, 303).
По договоренности мы говорим, что эти точки должны соответствовать трехмерным точкам: (0, 0, 0), (0, 0, 1), (1, 0, 0) и (1, 0, 1). Другими словами, единичный квадрат в плоскости y = 0.
Проецирование каждой из этих трехмерных координат в 2D выполняется путем умножения вектора 4D [x, y, z, 1] на матрицу проекции 4x4, а затем деления компонентов x и y на z, чтобы фактически получить коррекцию перспективы. Это более или менее то, что делает gluProject () , за исключением того, что gluProject() также учитывает текущий видовой экран и учитывает отдельную матрицу вида модели (мы можем просто предположить, что матрица вида модели - это единичная матрица). Очень удобно смотреть на документацию gluProject(), потому что я на самом деле хочу решение, которое работает для OpenGL, но учтите, что в документации отсутствует деление на z в формуле.
Помните, алгоритм должен начинаться с некоторой матрицы проекций и случайным образом возмущать ее, пока она не даст желаемую проекцию. Итак, что мы собираемся сделать, это спроектировать каждую из четырех трехмерных точек и посмотреть, насколько мы приблизились к желаемым 2D точкам. Если наши случайные возмущения приводят к тому, что проецируемые двумерные точки становятся ближе к тем, которые мы отмечали выше, мы сохраняем эту матрицу как улучшение по сравнению с нашим первоначальным (или предыдущим) предположением.
Давайте определим наши точки:
# Known 2D coordinates of our rectangle
i0 = Point2(318, 247)
i1 = Point2(326, 312)
i2 = Point2(418, 241)
i3 = Point2(452, 303)
# 3D coordinates corresponding to i0, i1, i2, i3
r0 = Point3(0, 0, 0)
r1 = Point3(0, 0, 1)
r2 = Point3(1, 0, 0)
r3 = Point3(1, 0, 1)
Нам нужно начать с некоторой матрицы, единичная матрица кажется естественным выбором:
mat = [
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1],
]
Нам нужно на самом деле реализовать проекцию (которая в основном является умножением матрицы):
def project(p, mat):
x = mat[0][0] * p.x + mat[0][1] * p.y + mat[0][2] * p.z + mat[0][3] * 1
y = mat[1][0] * p.x + mat[1][1] * p.y + mat[1][2] * p.z + mat[1][3] * 1
w = mat[3][0] * p.x + mat[3][1] * p.y + mat[3][2] * p.z + mat[3][3] * 1
return Point(720 * (x / w + 1) / 2., 576 - 576 * (y / w + 1) / 2.)
Это в основном то, что gluProject() делает, 720 и 576 - это ширина и высота изображения соответственно (то есть область просмотра), и мы вычитаем из 576, чтобы считать тот факт, что мы посчитали координаты y сверху, в то время как OpenGL обычно считает их снизу. Вы заметите, что мы не вычисляем z, это потому, что нам это здесь не нужно (хотя это может быть удобно, чтобы убедиться, что оно попадает в диапазон, который OpenGL использует для буфера глубины).
Теперь нам нужна функция для оценки того, насколько мы близки к правильному решению. Значение, возвращаемое этой функцией - это то, что мы будем использовать, чтобы проверить, является ли одна матрица лучше другой. Я решил пройти по сумме квадратов расстояний, т.е.
# The squared distance between two points a and b
def norm2(a, b):
dx = b.x - a.x
dy = b.y - a.y
return dx * dx + dy * dy
def evaluate(mat):
c0 = project(r0, mat)
c1 = project(r1, mat)
c2 = project(r2, mat)
c3 = project(r3, mat)
return norm2(i0, c0) + norm2(i1, c1) + norm2(i2, c2) + norm2(i3, c3)
Чтобы возмущать матрицу, мы просто выбираем элемент для возмущения на случайную величину в некотором диапазоне:
def perturb(amount):
from copy import deepcopy
from random import randrange, uniform
mat2 = deepcopy(mat)
mat2[randrange(4)][randrange(4)] += uniform(-amount, amount)
(Стоит отметить, что наша функция project() на самом деле вообще не использует mat[2], поскольку мы не вычисляем z, а все наши координаты y равны 0, значения mat[*][1] также не имеют значения Мы могли бы использовать этот факт и никогда не пытаться нарушать те значения, которые могли бы привести к небольшому ускорению, но это оставлено в качестве упражнения ...)
Для удобства давайте добавим функцию, которая выполняет основную часть аппроксимации, вызывая perturb() снова и снова для определения наилучшей из найденных нами матриц:
def approximate(mat, amount, n=100000):
est = evaluate(mat)
for i in xrange(n):
mat2 = perturb(mat, amount)
est2 = evaluate(mat2)
if est2 < est:
mat = mat2
est = est2
return mat, est
Теперь осталось только запустить ...:
for i in xrange(100):
mat = approximate(mat, 1)
mat = approximate(mat, .1)
Я считаю, что это уже дает довольно точный ответ. Через некоторое время я нашел следующую матрицу:
[
[1.0836000765696232, 0, 0.16272110011060575, -0.44811064935115597],
[0.09339193527789781, 1, -0.7990570384334473, 0.539087345090207 ],
[0, 0, 1, 0 ],
[0.06700844759602216, 0, -0.8333379578853196, 3.875290562060915 ],
]
с ошибкой около 2.6e-5. (Обратите внимание, что элементы, которые мы сказали, не использовались в вычислениях, на самом деле не были изменены по сравнению с нашей исходной матрицей; это потому, что изменение этих записей не приведет к изменению результата оценки, и поэтому изменение никогда не будет перенесено.)
Мы можем передать матрицу в OpenGL, используя glLoadMatrix() (но не забудьте сначала транспонировать ее, и не забудьте загрузить матрицу вида модели с матрицей идентификаторов):
def transpose(m):
return [
[m[0][0], m[1][0], m[2][0], m[3][0]],
[m[0][1], m[1][1], m[2][1], m[3][1]],
[m[0][2], m[1][2], m[2][2], m[3][2]],
[m[0][3], m[1][3], m[2][3], m[3][3]],
]
glLoadMatrixf(transpose(mat))
Теперь мы можем, например, перевести вдоль оси z, чтобы получить разные позиции вдоль дорожек:
glTranslate(0, 0, frame)
frame = frame + 1
glBegin(GL_QUADS)
glVertex3f(0, 0, 0)
glVertex3f(0, 0, 1)
glVertex3f(1, 0, 1)
glVertex3f(1, 0, 0)
glEnd()
Конечно, это не очень элегантно с математической точки зрения; вы не получите уравнение в замкнутой форме, в которое вы можете просто вставить свои числа и получить прямой (и точный) ответ. ОДНАКО, это позволяет вам добавлять дополнительные ограничения, не беспокоясь о усложнении ваших уравнений; например, если мы хотим также указать высоту, мы можем использовать этот угол дома и сказать (в нашей функции оценки), что расстояние от земли до крыши должно быть примерно одинаковым, и снова запустить алгоритм. Так что да, это своего рода грубая сила, но она работает и работает хорошо.
