Píxel retroproyectado a rayos 3D en coordenadas mundiales usando el método pseudoinverso

En lugar de tratar de depurar su código y verificar todas esas asignaciones inversas, describiré una forma de verificar sus propios resultados de manera objetiva. Si no tiene una buena idea de cuáles deberían ser los resultados, entonces realmente no veo cómo puede saber si son o no "razonables".

Asumiendo que no hay sesgo en la cámara, la matriz$K$tiene la forma

$$K=\begin{bmatrix}s_x&0&c_x\\0&s_y&c_y\\0&0&1\end{bmatrix}.$$ Los valores a lo largo de la diagonal son$x$- y$y$- factores de escala, y$(c_x,c_y)$son las coordenadas de la imagen del eje de la cámara, que se supone normal al plano de la imagen ($z=1$por convención). Entonces, en este sistema de coordenadas, el vector de dirección para un punto$(x,y)$en la imagen esta$(x-c_x,y-c_y,1)$y para obtener el vector de dirección correspondiente en el sistema de coordenadas de la cámara (externa), divida por los factores de escala respectivos:$((x-c_x)/s_x,(y-c_y)/s_y,1)$. Esto es exactamente lo que obtienes al aplicar$K^{-1}$, que se encuentra fácilmente que es $$K^{-1}=\begin{bmatrix}1/s_x&0&-c_x/s_x\\0&1/s_y&-c_y/s_y\\0&0&1\end{bmatrix}$$ usando tu método favorito. Finalmente, para transformar este vector en coordenadas mundiales, aplique$R^{-1}$, que es solo$R$'s transponer ya que es una rotación. El rayo resultante, por supuesto, se origina en la posición de la cámara en coordenadas mundiales. Debería ser sencillo codificar esta cascada explícitamente, después de lo cual puede compararla con los resultados que obtiene con cualquier otro método con el que esté experimentando.

En este caso concreto,$R$es solo la matriz de identidad, por lo que no hay nada más que hacer una vez que tenga el vector de dirección en las coordenadas de la cámara. Tenemos

$$s_x=282.363047 \\ s_y=280.10715905 \\ c_x=166.21515189 \\ c_y=108.05494375$$ por lo que la transformación de interno a externo es aproximadamente $$\begin{align}x&\to x/282.363-0.589 \\ y&\to y/280.107-0.386.\end{align}$$ Aplicando esto al punto$(20,20)$de tu pregunta anterior da$(-0.518,-0.314,1)$, que concuerda con el vector de dirección calculado allí. Tomando$(10,10)$en cambio resulta en$(-0.553,-0.350,1)$, que luego puede verificar con lo que haya producido su código, y así sucesivamente.

Aparte de todo eso, hay un problema cuando se usa el método pseudoinverso descrito por Zisserman. Da la siguiente ecuación para el rayo mapeado hacia atrás:

$$\mathbf X(\lambda)=P^+\mathbf x+\lambda\mathbf C.$$ Tenga en cuenta que el parámetro es un coeficiente de$\mathbf C$, la posición de la cámara en coordenadas mundiales, no del resultado de mapear hacia atrás el punto de la imagen$\mathbf x$. Convertido en coordenadas cartesianas, hay un factor de$\lambda+k$(para alguna constante$k$) en el denominador, por lo que esta no es una parametrización lineal simple. Para extraer un vector de dirección de esto, deberá convertir$P^+\mathbf x$en coordenadas cartesianas y luego restar$\mathbf C$.

Para ilustrar, aplicar$P^+$a$(10,10,1)$produce$(-0.553,-0.175,1.0,-0.175)$, entonces el rayo es$(-0.553,-t-0.175,1.0,t-0.175)$. En coordenadas cartesianas, el punto retroasignado es$(3.161,1.0,-5.713)$y restando la posición de la cámara da$(3.161,2.0,-5.713)$. Para comparar esto con el resultado conocido anterior, divida por la tercera coordenada:$(-0.553,-0.350,1.0)$, que está de acuerdo.

Actualización 2018.07.31: Para cámaras finitas, que es con lo que estás tratando, Zisserman sugiere una retroproyección más conveniente en el siguiente párrafo de la ecuación (6.14). La idea subyacente es que se descompone la matriz de la cámara como$P = \left[M\mid\mathbf p_4\right]$para que la retroproyección de un punto de imagen$\mathbf x$corta al plano en el infinito en$\mathbf D = ((M^{-1}\mathbf x)^T,0)^T$. Esto le da el vector de dirección del rayo retroproyectado en coordenadas mundiales y, por supuesto, el centro de la cámara está en$\tilde{\mathbf C}=-M^{-1}\mathbf p_4$, es decir, el rayo retroproyectado es

$$\tilde{\mathbf X}(\mu) = -M^{-1}\mathbf p_4+\mu M^{-1}\mathbf x = M^{-1}(\mu\mathbf x-\mathbf p_4).$$ Esta parametrización del rayo no sufre de la no linealidad mencionada anteriormente.

El método pseudoinverso funciona, a continuación se muestra el ejemplo para el píxel 215,180 (la esquina superior izquierda de la imagen es (0,0)), el rayo de este píxel va hacia la parte inferior derecha desde el punto de vista de una persona que mira desde el centro de la cámara hacia el eje Y. Debido al modelo de cámara estenopeica / proyección en perspectiva, fueron necesarios algunos cambios en el eje (pude cambiar mientras trazaba, pero el código a continuación es parte de otro análisis que tuve que realizar en un espacio 3D familiar).

from PIL import Image

from mpl_toolkits.mplot3d import Axes3D
import scipy.linalg as lin

K = [[ 282.363047,      0.,          166.21515189],
     [   0.,          280.10715905,  108.05494375],
     [   0.,            0.,            1.        ]]
K = np.array(K)
R = np.eye(3)
t = np.array([[0],[1.],[0]])
P = K.dot(np.hstack((R,t)))
C = np.array([0., 0., 1.])
p1 = np.array([215, 180, 1.])

X = np.dot(lin.pinv(P),p1)
X = X / X[3]
XX  = np.copy(X)
XX[1] = X[2]; XX[2] = X[1]; XX[2] = -XX[2]
w = 10
f = plt.figure()
ax = f.gca(projection='3d')
xvec = C - XX[:3] 
xvec = -xvec
ax.quiver(C[0], C[1], C[2], xvec[0], xvec[1], xvec[2],color='red')
ax.set_xlim(0,10);ax.set_ylim(0,10);ax.set_zlim(0,10)
ax.quiver(0., 0., 1., 0, 5., 0.,color='blue')
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_xlim(-w,w);ax.set_ylim(-w,w);ax.set_zlim(-w,w)
ax.view_init(elev=5, azim=100)
plt.savefig('out1.png')
ax.view_init(elev=5, azim=50)
plt.savefig('out2.png')

Dado que este es el primer y único resultado de Google que realmente ayuda, y tomó bastante tiempo crearlo, aquí está en Numpy

import numpy as np


def convert_nx2_to_homo_3xn(points_nx2: np.array):
    return np.hstack([points_nx2, np.ones((points_nx2.shape[0], 1))]).T


def convert_nx3_to_homo_4xn(points_nx3: np.array):
    return np.hstack([points_nx3, np.ones((points_nx3.shape[0], 1))]).T


def cast_2d_points_as_3d_rays(sub_pixels_nx2: np.array, proj_3x4: np.array):
    """
    see Harley & Zisserman pg 162, section 6.2.2, figure 6.14
    see https://math.stackexchange.com/a/597489/541203

    cast rays from camera center through the sub_pixels_nx2. Return camera center and ray directions.
    """
    m_3x3 = proj_3x4[:, :3]
    p4_3x1 = proj_3x4[:, 3]
    m_inv_3x3 = np.linalg.inv(m_3x3)

    # projection matrix to camera center
    camera_center_3x1 = np.expand_dims(-m_inv_3x3 @ p4_3x1, 1)
    camera_center_homo_4x1 = np.vstack([camera_center_3x1, 1])

    # projection matrix + pixel locations to ray directions
    sub_pixels_homo_3xn = convert_nx2_to_homo_3xn(sub_pixels_nx2)
    ray_directions_3xn = m_inv_3x3 @ sub_pixels_homo_3xn
    ray_directions_homo_4xn = convert_nx3_to_homo_4xn(ray_directions_3xn.T)
    ray_directions_homo_4xn[3, :] = 0  # this is a direction

    return ray_directions_homo_4xn, camera_center_homo_4x1

```