it-swarm-es.tech

¿Cómo proyecto en reversa los puntos 2D en 3D?

Tengo 4 puntos 2D en el espacio de la pantalla, y necesito proyectarlos hacia atrás en el espacio 3D. Sé que cada uno de los 4 puntos es una esquina de un rectángulo rígido girado en 3D, y sé el tamaño del rectángulo. ¿Cómo puedo obtener coordenadas 3D de esto?

No estoy usando ninguna API en particular, y no tengo una matriz de proyección existente. Solo estoy buscando matemática básica para hacer esto. Por supuesto, no hay suficientes datos para convertir un solo punto 2D a 3D sin otra referencia, pero imagino que si tiene 4 puntos, sabe que están todos en ángulo recto entre sí en el mismo plano, y conoces la distancia entre ellos, deberías poder averiguarlo desde allí. Desafortunadamente, no puedo entender cómo.

Esto podría caer bajo el paraguas de la fotogrametría, pero las búsquedas en Google no me han llevado a ninguna información útil.

54
Joshua Carmody

Bien, vine aquí buscando una respuesta y no encontré algo simple y directo, así que seguí adelante e hice lo tonto pero efectivo (y relativamente simple): la optimización de Monte Carlo.

En pocas palabras, el algoritmo es el siguiente: perturbar aleatoriamente su matriz de proyección hasta que proyecte sus coordenadas 3D conocidas a sus coordenadas 2D conocidas.

Aquí hay una foto fija de Thomas the Tank Engine:

Thomas the Tank Engine

Digamos que usamos GIMP para encontrar las coordenadas 2D de lo que creemos que es un cuadrado en el plano del suelo (si realmente es un cuadrado o no depende de su juicio de la profundidad):

With an outline of the square

Obtengo cuatro puntos en la imagen 2D: (318, 247), (326, 312), (418, 241) y (452, 303).

Por convención, decimos que estos puntos deben corresponder a los puntos 3D: (0, 0, 0), (0, 0, 1), (1, 0, 0) y (1, 0, 1). En otras palabras, una unidad cuadrada en el plano y = 0.

La proyección de cada una de estas coordenadas 3D en 2D se realiza multiplicando el vector 4D [x, y, z, 1] con una matriz de proyección 4x4, luego dividiendo los componentes x e y por z para obtener la corrección de la perspectiva. Esto es más o menos lo que gluProject () hace, excepto que gluProject() también tiene en cuenta la ventana gráfica actual y tiene en cuenta una matriz de vista de modelo separada (podemos suponer que la matriz de vista de modelo es la matriz de identidad). Es muy útil mirar la documentación de gluProject() porque realmente quiero una solución que funcione para OpenGL, pero tenga en cuenta que a la documentación le falta la división por z en la fórmula.

Recuerde, el algoritmo es comenzar con alguna matriz de proyección y perturbarla aleatoriamente hasta obtener la proyección que deseamos. Entonces, lo que vamos a hacer es proyectar cada uno de los cuatro puntos 3D y ver qué tan cerca nos acercamos a los puntos 2D que queríamos. Si nuestras perturbaciones aleatorias hacen que los puntos 2D proyectados se acerquen a los que marcamos anteriormente, entonces mantenemos esa matriz como una mejora sobre nuestra suposición inicial (o anterior).

Definamos nuestros puntos:

# 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)

Necesitamos comenzar con alguna matriz, la matriz de identidad parece una elección natural:

mat = [
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, 1, 0],
    [0, 0, 0, 1],
]

Necesitamos implementar realmente la proyección (que es básicamente una multiplicación matricial):

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.)

Esto es básicamente lo que hace gluProject(), 720 y 576 son el ancho y la altura de la imagen, respectivamente (es decir, la ventana gráfica), y restamos de 576 para contar el hecho de que contamos las coordenadas y desde la parte superior, mientras que OpenGL generalmente las cuenta desde El fondo. Notarás que no estamos calculando z, eso es porque realmente no lo necesitamos aquí (aunque podría ser útil para garantizar que se encuentre dentro del rango que OpenGL usa para el búfer de profundidad).

Ahora necesitamos una función para evaluar qué tan cerca estamos de la solución correcta. El valor devuelto por esta función es lo que usaremos para verificar si una matriz es mejor que otra. Elegí ir por la suma de las distancias al cuadrado, es decir:

# 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)

Para perturbar la matriz, simplemente elegimos un elemento para perturbar por una cantidad aleatoria dentro de cierto rango:

def perturb(amount):
    from copy import deepcopy
    from random import randrange, uniform
    mat2 = deepcopy(mat)
    mat2[randrange(4)][randrange(4)] += uniform(-amount, amount)

(Vale la pena señalar que nuestra función project() en realidad no usa mat[2] en absoluto, ya que no calculamos z, y dado que todas nuestras coordenadas y son 0, los valores mat[*][1] tampoco son relevantes. podría usar este hecho y nunca tratar de perturbar esos valores, lo que daría una pequeña aceleración, pero eso se deja como un ejercicio ...)

Por conveniencia, agreguemos una función que haga la mayor parte de la aproximación llamando a perturb() una y otra vez en cuál es la mejor matriz que hemos encontrado hasta ahora:

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

Ahora todo lo que queda por hacer es ejecutarlo ...:

for i in xrange(100):
    mat = approximate(mat, 1)
    mat = approximate(mat, .1)

Creo que esto ya da una respuesta bastante precisa. Después de correr por un tiempo, la matriz que encontré fue:

[
    [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  ],
]

con un error de alrededor de 2.6e-5. (Observe cómo los elementos que dijimos que no se usaron en el cálculo en realidad no se han cambiado desde nuestra matriz inicial; eso es porque cambiar estas entradas no cambiaría el resultado de la evaluación y, por lo tanto, el cambio nunca se llevaría a cabo)

Podemos pasar la matriz a OpenGL usando glLoadMatrix() (pero recuerde transponerla primero, y recuerde cargar su matriz de vista de modelo con la matriz de identidad):

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))

Ahora podemos, por ejemplo, traducir a lo largo del eje z para obtener diferentes posiciones a lo largo de las pistas:

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()

With 3D translation

Por supuesto, esto no es muy elegante desde un punto de vista matemático; no obtienes una ecuación de forma cerrada en la que solo puedes conectar tus números y obtener una respuesta directa (y precisa). SIN EMBARGO, le permite agregar restricciones adicionales sin tener que preocuparse por complicar sus ecuaciones; Por ejemplo, si también quisiéramos incorporar la altura, podríamos usar esa esquina de la casa y decir (en nuestra función de evaluación) que la distancia desde el suelo hasta el techo debería ser regular, y ejecutar el algoritmo nuevamente. Entonces, sí, es una especie de fuerza bruta, pero funciona y funciona bien.

Choo choo!

62
Vegard

Este es el problema clásico para la realidad aumentada basada en marcadores.

Tiene un marcador cuadrado (código de barras 2D) y desea encontrar su Pose (traslación y rotación en relación con la cámara), después de encontrar los cuatro bordes del marcador. Descripción general

No estoy al tanto de las últimas contribuciones al campo, pero al menos hasta cierto punto (2009) se suponía que el RPP superaba a POSIT que se menciona anteriormente (y de hecho es un enfoque clásico para esto). Consulte los enlaces, también proporcionar fuente.

(PD: Sé que es un tema un poco viejo, pero de todos modos, la publicación podría ser útil para alguien)

7
dim_tz

D. DeMenthon ideó un algoritmo para calcular la pose de un objeto (su posición y orientación en el espacio) a partir de puntos característicos en una imagen 2D al conocer el modelo del objeto - este es su problema exacto :

Describimos un método para encontrar la pose de un objeto a partir de una sola imagen. Suponemos que podemos detectar y hacer coincidir en la imagen cuatro o más puntos de características no coplanares del objeto, y que conocemos su geometría relativa en el objeto.

El algoritmo se conoce como Posit y se describe en su artículo clásico "Pose de objeto basado en modelo en 25 líneas de código" (disponible en su sitio web , sección 4).

Enlace directo al artículo: http://www.cfar.umd.edu/~daniel/daniel_papersfordownload/Pose25Lines.pdf Implementación de OpenCV: http://opencv.willowgarage.com/ wiki/Posit

La idea es aproximar repetidamente la proyección en perspectiva mediante una proyección ortográfica a escala hasta converger en una pose precisa.

5
Julien-L

Desde el espacio 2-D habrá 2 rectángulos válidos que se pueden construir. Sin conocer la proyección matricial original, no sabrá cuál es la correcta. Es lo mismo que el problema de la "caja": ves dos cuadrados, uno dentro del otro, con los 4 vértices internos conectados a los 4 vértices externos respectivos. ¿Estás mirando una caja de arriba hacia abajo o de abajo hacia arriba?

Dicho esto, está buscando una matriz de transformación T donde ...

{{x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3}, {x4, y4, z4}} x T = {{x1, y1}, {x2, y2}, { x3, y3}, {x4, y4}}

(4 x 3) x T = (4 x 2)

Entonces T debe ser una matriz (3 x 2). Entonces tenemos 6 incógnitas.

Ahora construya un sistema de restricciones en T y resuelva con Simplex. Para construir las restricciones, usted sabe que una línea que pasa por los primeros dos puntos debe ser paralela a la línea que pasa a los segundos dos puntos. Usted sabe que una línea que pasa por los puntos 1 y 3 debe ser paralela a las líneas que pasan por los puntos 2 y 4. Usted sabe que una línea que pasa por 1 y 2 debe ser ortogonal a una línea que pasa por los puntos 2 y 3. Usted sabe que la longitud de la línea de 1 y 2 debe ser igual a la longitud de la línea de 3 y 4. Usted sabe que la longitud de la línea de 1 y 3 debe ser igual a la longitud de la línea de 2 y 4.

Para hacer esto aún más fácil, conoce el rectángulo, de modo que conoce la longitud de todos los lados.

Eso debería darle muchas restricciones para resolver este problema.

Por supuesto, para volver, puedes encontrar T-inverso.

@Rob: Sí, hay un número infinito de proyecciones, pero no un número infinito de proyectos donde los puntos deben satisfacer los requisitos de un rectángulo.

@nlucaroni: Sí, esto solo se puede resolver si tiene cuatro puntos en la proyección. Si el rectángulo se proyecta a solo 2 puntos (es decir, el plano del rectángulo es ortogonal a la superficie de proyección), entonces esto no se puede resolver.

Hmmm ... debería irme a casa y escribir esta pequeña joya. Esto suena divertido.

Actualizaciones:

  1. Hay un número infinito de proyecciones a menos que arregles uno de los puntos. Si fija los puntos del rectángulo original, entonces hay dos posibles rectángulos originales.
4
Jarrett Meyer

Para seguir el enfoque de Rons: puede encontrar sus valores z si sabe cómo ha rotado su rectángulo.

El truco es encontrar la matriz proyectiva que hizo la proyección. Afortunadamente esto es posible e incluso barato de hacer. La matemática relevante se puede encontrar en el documento "Mapeos proyectivos para deformación de imágenes" de Paul Heckbert.

http://pages.cs.wisc.edu/~dyer/cs766/readings/heckbert-proj.pdf

De esta manera, puede recuperar la parte homogénea de cada vértice que se perdió durante la proyección.

Ahora todavía te quedan cuatro líneas en lugar de puntos (como explicó Ron). Sin embargo, dado que conoce el tamaño de su rectángulo original, no se pierde nada. Ahora puede conectar los datos del método de Ron y del enfoque 2D a un solucionador de ecuaciones lineales y resolver para z. Obtiene los valores z exactos de cada vértice de esa manera.

Nota: Esto solo funciona porque:

  1. La forma original era un rectángulo.
  2. Conoces el tamaño exacto del rectángulo en el espacio 3D.

Es un caso especial realmente.

Espero que ayude, Nils

2
Nils Pipenbrinck

Suponiendo que los puntos son realmente parte de un rectángulo, estoy dando una idea genérica:

Encuentre dos puntos con una distancia máxima: estos probablemente definan una diagonal (excepción: casos especiales donde el rectángulo es casi paralelo al plano YZ, dejado para el estudiante). Llámalos A, C. Calcula los ángulos MALO, BCD. Estos, en comparación con los ángulos rectos, le dan orientación en el espacio 3D. Para obtener información sobre la distancia z, debe correlacionar los lados proyectados con los lados conocidos, y luego, basándose en el método de proyección 3d (¿es 1/z?), Está en el camino correcto para conocer las distancias.

2
tzot

Sacaré mi libro de Álgebra lineal cuando llegue a casa si nadie responde. Pero @ D G, no todas las matrices son invertibles. Las matrices singulares no son invertibles (cuando determinante = 0). Esto realmente sucederá todo el tiempo, ya que una matriz de proyección debe tener valores propios de 0 y 1, y ser cuadrada (ya que es idempotente, entonces p ^ 2 = p).

Un ejemplo fácil es, [[0 1] [0 1]] ya que el determinante = 0, ¡y esa es una proyección en la línea x = y!

1
nlucaroni

La proyección que tiene sobre la superficie 2D tiene infinitos rectángulos 3D que se proyectarán a la misma forma 2D.

Piénselo de esta manera: tiene cuatro puntos 3D que forman el rectángulo 3D. Llámalos (x0, y0, z0), (x1, y1, z1), (x2, y2, z2) y (x3, y3, z3). Cuando proyecta estos puntos en el plano x-y, suelta las coordenadas z: (x0, y0), (x1, y1), (x2, y2), (x3, y3).

Ahora, si desea volver a proyectar en el espacio 3D, debe realizar una ingeniería inversa de lo que eran z0, ..., z3. Pero cualquier conjunto de coordenadas z que a) mantenga la misma distancia x-y entre los puntos, yb) mantenga la forma en que funcionará un rectángulo. Entonces, cualquier miembro de este conjunto (infinito) hará: {(z0 + i, z1 + i, z2 + i, z3 + i) | i <- R}.

Editar @Jarrett: Imagina que resolviste esto y terminaste con un rectángulo en el espacio 3D. Ahora, imagina deslizar ese rectángulo hacia arriba y hacia abajo en el eje z. Esas infinitas cantidades de rectángulos traducidos tienen la misma proyección x-y. ¿Cómo sabes que encontraste el "correcto"?

Editar # 2: Muy bien, esto es de un comentario que hice sobre esta pregunta, un enfoque más intuitivo para razonar sobre esto.

Imagínese sosteniendo un pedazo de papel sobre su escritorio. Imagina que cada esquina del papel tiene un puntero láser sin peso que apunta hacia el escritorio. El papel es el objeto 3D, y los puntos del puntero láser en el escritorio son la proyección 2D.

Ahora, ¿cómo puede saber qué tan alto del escritorio está el papel mirando solo los puntos del puntero láser?

No puedes Mueva el papel hacia arriba y hacia abajo. Los punteros láser seguirán brillando en los mismos puntos en el escritorio, independientemente de la altura del papel.

Encontrar las coordenadas z en la proyección inversa es como tratar de encontrar la altura del papel basándose solo en los puntos del puntero láser en el escritorio.

1
Rob Dickerson
1
Ray Tayek

Si sabe que la forma es un rectángulo en un plano, puede limitar aún más el problema. Ciertamente no puede averiguar "qué" plano, por lo que puede elegir que se encuentre en el plano donde z = 0 y una de las esquinas está en x = y = 0, y los bordes son paralelos al eje x/y.

Los puntos en 3d son, por lo tanto, {0,0,0}, {w, 0,0}, {w, h, 0} y {0, h, 0}. Estoy bastante seguro de que no se encontrará el tamaño absoluto, por lo que solo la relación w/h es relevante, por lo que esta es una incógnita.

En relación con este plano, la cámara debe estar en algún punto cx, cy, cz en el espacio, debe estar apuntando en una dirección nx, ny, nz (un vector de longitud uno, por lo que uno de estos es redundante) y tener una longitud focal/ancho_de_imagen factor de w. Estos números se convierten en una matriz de proyección 3x3.

Eso da un total de 7 incógnitas: w/h, cx, cy, cz, nx, ny y w.

Tienes un total de 8 conocimientos: los 4 pares x + y.

Entonces esto se puede resolver.

El siguiente paso es usar Matlab o Mathmatica.

1
spitzak

Sí, Monte Carlo funciona, pero encontré una mejor solución para este problema. Este código funciona perfectamente (y usa OpenCV):

Cv2.CalibrateCamera(new List<List<Point3f>>() { points3d }, new List<List<Point2f>>() { points2d }, new Size(height, width), cameraMatrix, distCoefs, out rvecs, out tvecs, CalibrationFlags.ZeroTangentDist | CalibrationFlags.FixK1 | CalibrationFlags.FixK2 | CalibrationFlags.FixK3);

Esta función toma los puntos 3d y 2d conocidos, el tamaño de la pantalla y devuelve la rotación (rvecs [0]), la traducción (tvecs [0]) y la matriz de valores intrínsecos de la cámara. Es todo lo que necesitas.

1
Inflight

Cuando proyecta de 3D a 2D, pierde información.

En el caso simple de un solo punto, la proyección inversa le daría un rayo infinito a través del espacio 3d.

La reconstrucción estereoscópica generalmente comenzará con dos imágenes 2d y se proyectará de nuevo a 3D. Luego busque una intersección de los dos rayos 3D producidos.

La proyección puede tomar diferentes formas. Ortogonal o perspectiva. ¿Supongo que estás asumiendo una proyección ortogonal?

En su caso, suponiendo que tuviera la matriz original, tendría 4 rayos en el espacio 3D. Entonces podrá restringir el problema mediante las dimensiones de su rectángulo 3d e intentar resolverlo.

La solución no será única, ya que una rotación alrededor de cualquier eje que sea paralela al plano de proyección 2D será ambigua en su dirección. En otras palabras, si la imagen 2d es perpendicular al eje z, al girar el rectángulo 3d en sentido horario o antihorario alrededor del eje x se produciría la misma imagen. Del mismo modo para el eje y.

En el caso de que el plano del rectángulo sea paralelo al eje z, tiene aún más soluciones.

Como no tiene la matriz de proyección original, se introduce una mayor ambigüedad por un factor de escala arbitrario que existe en cualquier proyección. No puede distinguir entre una escala en la proyección y una traslación en 3d en la dirección del eje z. Esto no es un problema si solo está interesado en las posiciones relativas de los 4 puntos en el espacio 3d cuando se relacionan entre sí y no con el plano de la proyección 2d.

En una perspectiva, las cosas se ponen más difíciles ...

1
morechilli

Gracias a @Vegard por una excelente respuesta. Limpié un poco el código:

import pandas as pd
import numpy as np

class Point2:
    def __init__(self,x,y):
        self.x = x
        self.y = y

class Point3:
    def __init__(self,x,y,z):
        self.x = x
        self.y = y
        self.z = z

# 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):
    #print 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 Point2(720 * (x / w + 1) / 2., 576 - 576 * (y / w + 1) / 2.)

# 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(mat, amount):
    from copy import deepcopy
    from random import randrange, uniform
    mat2 = deepcopy(mat)
    mat2[randrange(4)][randrange(4)] += uniform(-amount, amount)
    return mat2

def approximate(mat, amount, n=1000):
    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(1000):
    mat,est = approximate(mat, 1)
    print mat
    print est

La llamada aproximada con .1 no funcionó para mí, así que la saqué. Lo ejecuté por un tiempo también, y la última vez que lo revisé fue a las

[[0.7576315397559887, 0, 0.11439449272592839, -0.314856490473439], 
[0.06440497208710227, 1, -0.5607502645413118, 0.38338196981556827], 
[0, 0, 1, 0], 
[0.05421620936883742, 0, -0.5673977598434641, 2.693116299312736]]

con un error de alrededor de 0.02.

1
user423805