Аффинное 3D-преобразование в Python

Вот решение с использованием numpy и scipy. scipy в основном используется для генерации случайных вращений, за исключением scipy.linalg.norm, который легко написать самостоятельно. Основные вещи, используемые из numpy, — это перекрестное произведение и умножение матриц, которые также легко кодировать самостоятельно.

Основная идея такова: по трем неколлинеарным точкам x1,x2,x3 можно найти ортогональную тройку векторов (осей) v1,v2,v3, причем v1 направлена ​​в направлении x2-x1, v2 в плоскости натянутый (x2-x1) и (x3-x1), и v3 завершает тройку.

Точки y1,y2,y3 поворачиваются и перемещаются относительно x1,x2,x3. Оси w1,w2,w3, сгенерированные из y1,y2,y3, вращаются (т. е. без перемещения) из v1,v2,v3. Каждый из этих двух наборов троек ортогонален, поэтому по ним легко найти поворот: R = W * transpose(V)

Когда у нас есть вращение, найти перевод просто: y1 = R*x + t, поэтому t = y1 - R*x. Возможно, было бы лучше использовать решатель наименьших квадратов и объединить все три точки, чтобы получить оценку t.

import numpy
import scipy.linalg


def rand_rot():
    """Return a random rotation

    Return a random orthogonal matrix with determinant 1"""
    q, _ = scipy.linalg.qr(numpy.random.randn(3, 3))
    if scipy.linalg.det(q) < 0:
        # does this ever happen?
        print "got a negative det"
        q[:, 0] = -q[:, 0]
    return q


def rand_noncollinear():
    """Return 3 random non-collinear vectors"""
    while True:
        b = numpy.random.randn(3, 3)
        sigma = scipy.linalg.svdvals(b)
        if sigma[2]/sigma[0] > 0.1:
            # "very" non-collinear
            break
        # "nearly" collinear; try again

    return b[:, 0], b[:, 1], b[:, 2]


def normalize(a):
    """Return argument normalized"""
    return a/scipy.linalg.norm(a)


def colstack(a1, a2, a3):
    """Stack three vectors as columns"""
    return numpy.hstack((a1[:, numpy.newaxis],
                         a2[:, numpy.newaxis],
                         a3[:, numpy.newaxis]))


def get_axes(a1, a2, a3):
    """Generate orthogonal axes from three non-collinear points"""
    # I tried to do this with QR, but something didn't work
    b1 = normalize(a2-a1)
    b2 = normalize(a3-a1)
    b3 = normalize(numpy.cross(b1, b2))
    b4 = normalize(numpy.cross(b3, b1))
    return b1, b4, b3

# random rotation and translation
r = rand_rot()
t = numpy.random.randn(3)

# three non-collinear points
x1, x2, x3 = rand_noncollinear()
# some other point
x4 = numpy.random.randn(3)

# the images of the above in the transformation.
# y4 is for checking only -- won't be used to estimate r or t
y1, y2, y3, y4 = [numpy.dot(r, x) + t
                  for x in x1, x2, x3, x4]


v1, v2, v3 = get_axes(x1, x2, x3)
w1, w2, w3 = get_axes(y1, y2, y3)

V = colstack(v1, v2, v3)
W = colstack(w1, w2, w3)

# W = R V, so R = W * inverse(V); but V orthogonal, so inverse(V) is
# transpose(V):
rfound = numpy.dot(W, V.T)

# y1 = R x1 + t, so...
tfound = y1-numpy.dot(r, x1)

# get error on images of x2 and x3, just in case

y2err = scipy.linalg.norm(numpy.dot(rfound, x2) + tfound - y2)
y3err = scipy.linalg.norm(numpy.dot(rfound, x3) + tfound - y3)

# and check error image of x4 -- getting an estimate of y4 is the
# point of all of this
y4err = scipy.linalg.norm(numpy.dot(rfound, x4) + tfound - y4)

print "y2 error: ", y2err
print "y3 error: ", y3err
print "y4 error: ", y4err

И описание, и ваш код сбивают с толку. Описание немного расплывчато, а в примерах кода отсутствуют важные фрагменты. Итак, как я понимаю вопрос:

Зная три точки в двух пространствах, как построить преобразование из пространства A в пространство B?

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

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

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

транс = q0 - p0

Это оставляет вам вращение, которое было объяснено в нескольких сообщениях:

• python + maya: повернуть ось Y вдоль вектор
• Как преобразовать трехмерный вектор в вращение Эйлера в таких программах, как Maya, с помощью python

После этого вы также можете рассчитать коэффициент масштабирования.

я понял это

p0p1 = p1-p0
p0p2 = p2-p0
p0p3 = p3-p0

q0q1 = q1-q0
q0q2 = q2-q0
q0q3 = q3-q0

before = dt.Matrix([p0.x, p0.y, p0.z, 0],[p1.x, p1.y, p1.z, 0],[p2.x, p2.y, p2.z, 0], [0,0,0,1]);
after = dt.Matrix([q0.x, q0.y, q0.z, 0],[q1.x, q1.y, q1.z, 0],[q2.x, q2.y, q2.z, 0], [0,0,0,1]);

normal = p0p1.cross(p0p2).normal()
dist = p0p3.dot(normal)
q3 = p3 - dist*normal

transformMatrix = before.inverse()*after
solve = dt.Matrix(q3.x, q3.y, q3.z, 1)*transformMatrix

q3 =  dt.Vector(solve[0][0], solve[0][1], solve[0][2])

newNormal = q0q1.cross(q0q2).normal()
q3 = q3 + newNormal*dist

pm.move(list[3], q3, r=False)

Матрица преобразования работала только для точек, находящихся в пределах плоскости p0p1p2. Поэтому я решил это, преобразовав спроецированную точку p3, а затем переместив ее от плоскости на такое же расстояние.

Если у вас есть решение, которое включает только матрицу, не стесняйтесь поделиться, оно все равно может мне помочь! :)