Calibración de cámaras

1Qué es calibrar

Los apuntes de transformaciones 3D terminaron con un modelo de la cámara:

Lleva un punto M′ = [X′ Y′ Z′ 1]^T, dado en las coordenadas de algún objeto, a su píxel w = [u v 1]^T. Tiene once parámetros. Cinco son intrínsecos, α, β, γ, u₀ y v₀, en K y P. Seis son extrínsecos, los ángulos ω_X, ω_Y, ω_Z y las tres componentes de t′, en S′. Calibrar es estimarlos a partir de imágenes.

La idea es fotografiar un objeto cuya forma se conoce con exactitud, un objeto de calibración, y elegir N puntos en él. Sus coordenadas M′_i se conocen por el objeto, en milímetros, y sus píxeles w_i se miden en la imagen. Reúne los parámetros en un vector θ. Para cualquier valor supuesto de θ, el modelo predice dónde debería aparecer cada punto:

donde F construye K, P y S′ a partir de θ, las multiplica y divide por λ. Un buen θ deja cada ŵ_i cerca del w_i medido, así que calibrar es buscar el θ que hace mínima la distancia promedio entre ellos:

J se mide en píxeles, así que dice directamente cuánto se equivoca el modelo, en promedio.

Esto se llama calibración fotogramétrica, porque se apoya en un objeto de tamaño conocido. La autocalibración prescinde de él: solo relaciona los mismos puntos vistos en varias imágenes. Pero no puede saber el tamaño real de la escena, porque un objeto pequeño cerca de la cámara da las mismas imágenes que una copia más grande más lejos. Cuando el resultado debe estar en milímetros, como para un robot que toma piezas, se usa un objeto de calibración.

La mesa

El ejemplo de toda esta página es una mesa de 800 mm de fondo, 1250 mm de ancho y 900 mm de alto, fotografiada de frente en una imagen de 1842 × 1408 píxeles. Sus ocho puntos conocidos son los extremos superior e inferior de sus cuatro patas, numerados del 0 al 7. El sistema de coordenadas de la mesa tiene su origen en la base de la pata trasera izquierda, la esquina 3, con X′ apuntando desde las patas traseras hacia las delanteras y la cámara, Y′ hacia la derecha a lo largo de la mesa, y Z′ hacia arriba, un sistema de mano derecha. Así, la base de la pata delantera izquierda, la esquina 0, está en (800, 0, 0), y el extremo superior de la pata delantera derecha, la esquina 5, está en (800, 1250, 900). Los píxeles se marcaron a mano sobre la foto, así que cada uno se desvía unos pocos píxeles.

2Calibración lineal

Al multiplicarla, KPS′ es una sola matriz de 3 × 4. Llámala A, la matriz de proyección general (en código suele llamarse P):

En vez de los once parámetros, el método lineal estima las doce entradas de A. No todas son libres: A y kA dan los mismos píxeles para cualquier número k ≠ 0, porque k solo cambia λ. Así que una entrada se puede fijar, y lo usual es a₃₄ = 1. Quedan once incógnitas, tantas como parámetros tiene la cámara.

¿Por qué se puede fijar a₃₄? La tercera fila de A da λ. Para el origen del objeto, M′ = [0 0 0 1]^T, da λ = a₃₄, así que a₃₄ es la profundidad de ese origen delante de la cámara, multiplicada por la escala de A. Fijarla en 1 solo elige la escala. La única excepción es un a₃₄ verdadero igual a 0, cuando el origen del objeto está a la altura del centro óptico, en el plano que pasa por C paralelo a la imagen. Ninguna escala convierte 0 en 1. Un origen que esté a la vista, como en la mesa, evita este caso.

Dos ecuaciones por punto

La tercera fila da λ = a₃₁X′ + a₃₂Y′ + a₃₃Z′ + 1. Al reemplazarla en las dos primeras filas, λ desaparece:

Ambas ecuaciones son lineales en las incógnitas, porque u, v, X′, Y′ y Z′ son números conocidos. Reúne las incógnitas en a = [a₁₁ a₁₂ … a₃₃]^T, once entradas fila por fila, y pásalas a la izquierda:

Cada punto agrega estas dos filas, y N puntos dan un sistema lineal

con 2N ecuaciones, donde Q tiene once columnas y r lista los u y v medidos.

Cuántos puntos

Once incógnitas necesitan al menos once ecuaciones, así que al menos seis puntos: cinco dan solo diez. Y los puntos no pueden estar todos en un plano. Si lo estuvieran, por ejemplo todos en el suelo, donde Z′ = 0, entonces a₁₃, a₂₃ y a₃₃ multiplicarían a 0 en todas las filas, y nada los determinaría. Los puntos de un plano solo determinan la homografía entre ese plano y la imagen, como en los apuntes de transformaciones 2D, y una homografía tiene ocho números, no once. Las ocho esquinas de la mesa están en dos niveles, el suelo y la cubierta, y sí determinan A.

Mínimos cuadrados

Los píxeles marcados se desvían unos pocos píxeles, así que con más de seis puntos ningún a cumple todas las ecuaciones. Como con las homografías, se toma el a que hace mínima la suma de cuadrados ‖Qa − r‖²:

Las columnas de Q difieren en tamaño en millones: un 1 junto a −uY′, que aquí llega a dos millones. Formar Q^TQ eleva al cuadrado esa diferencia y puede perder dígitos, así que las bibliotecas numéricas, como np.linalg.lstsq, llegan al mismo â sin formarla.

¿Qué minimizan aquí los mínimos cuadrados? La primera ecuación del punto i deja el residuo (primera fila de A)M′_i − u_iλ_i. El primer término es λ_iû_i, así que el residuo es λ_i(û_i − u_i): el error en píxeles, multiplicado por λ_i. Con a₃₄ = 1, λ_i es la profundidad del punto dividida por la profundidad del origen del objeto, así que los puntos más lejanos pesan más. El resultado se acerca al θ con el menor J, pero no es el mismo, y esta es una de las razones para el método no lineal de la sección 4.

La elección a₃₄ = 1 también se puede evitar, como con las homografías. Deja las doce entradas como incógnitas, agrega −u y −v como duodécima columna de Q, y resuelve Qa = 0 con ‖a‖ = 1. La respuesta es el vector singular derecho de Q con el menor valor singular. Para la mesa, las dos A coinciden en tres cifras significativas, y sus errores promedio difieren en menos de 0.001 píxeles.

3Dónde está la cámara, y sus once parámetros

A contiene la cámara completa, pero sus entradas lo mezclan todo. Dos hechos del modelo de cámara la separan. Escribe A = [B b], con B su bloque izquierdo de 3 × 3 y b su última columna.

El centro óptico

El centro óptico C no tiene imagen: todos los rayos parten de él, así que A lo envía al vector cero. En las coordenadas del objeto, AC = BC + b = 0, y por lo tanto

Esto no necesita ningún parámetro: encuentra dónde estaba la cámara a partir de A sola.

Separar K y R′

El modelo dice que A es un múltiplo de KPS′. Como KP = [K 0], ese producto es K[R′ t′], y

para algún número k. Dos propiedades separan B. K es triangular superior, con última fila [0 0 1]. Y las filas de R′, llámalas r₁, r₂ y r₃, son los ejes X, Y y Z de la cámara escritos en el sistema del objeto: vectores unitarios, perpendiculares entre sí. Llama m₁, m₂ y m₃ a las filas de B. Fila por fila desde abajo, B = kKR′ dice

y las incógnitas salen de a una.

• r₃ tiene largo 1, así que k = ±‖m₃‖ y r₃ = m₃/k.
• Un producto punto con r₃ extrae la parte de m₂/k a lo largo de r₃, porque r₂ es perpendicular a él: v₀ = m₂·r₃/k. Lo que queda, m₂/k − v₀r₃, es βr₂, cuyo largo es β y cuya dirección es r₂.
• Del mismo modo, u₀ = m₁·r₃/k y γ = m₁·r₂/k, y lo que queda de m₁/k es αr₁.

Esto es el proceso de Gram–Schmidt aplicado a las filas de B, de abajo hacia arriba. El signo de k es el que deja el objeto delante de la cámara, con λ positivo para sus puntos. Luego t′ = K⁻¹b/k, y los ángulos salen de R = R′^T, cuya primera fila es

Esta da ω_Y y ω_Z, y ω_X sale de la última columna de R.

Los largos son positivos, así que α y β siempre salen positivos, y un espejo en los datos aparece en R′: si det R′ = −1, R′ convierte una mano derecha en una izquierda, y ninguna cámara real da esos píxeles.

Hacia dónde apuntan los seis números

En esta página, como en los apuntes de transformaciones 3D, el sistema de coordenadas de la cámara es el original y el del objeto es el nuevo. Los ángulos ω giran los ejes de la cámara hasta los del objeto, y t′ es donde la cámara ve el origen del objeto. Muchos programas describen la misma cámara al revés, con el sistema de la cámara como el nuevo: un punto M′ del objeto está en R′M′ + t′ en el sistema de la cámara. Su traslación es el mismo t′, pero su matriz de rotación es R′, la traspuesta de R, así que sus tres ángulos son números distintos para la misma cámara. Antes de usar ángulos de otra fuente, revisa cuál de las dos matrices construyen.

Para la mesa, esto importa de otra forma más. Su eje X′ apunta casi directamente hacia la cámara, así que ω_Y sale cerca de 90°, cerca del bloqueo del cardán de los apuntes de transformaciones 3D: con 90°, ω_X y ω_Z girarían en torno a la misma recta. A 79° los ángulos siguen bien definidos, pero algunos cambios pequeños se agrandan. Mover la esquina 5 diez píxeles hacia abajo gira la cámara medio grado, pero cambia ω_X en 2.5° y ω_Z en 2.3°, ambos en el mismo sentido y unas cinco veces más. Esta es una de las razones por las que los programas de calibración guardan una rotación como un eje y un ángulo.

4Calibración no lineal

El método lineal nunca mira los parámetros. Estima doce números y minimiza un error en las ecuaciones, no en la imagen. La calibración no lineal trabaja con los parámetros mismos,

y minimiza J(θ), la distancia promedio en píxeles. Para cada θ, F construye K, convierte los ángulos en R = R_XR_YR_Z, forma S′ = [R^T t′; 0^T 1], y proyecta cada punto con λ_iŵ_i = KPS′M′_i.

Los ángulos entran a través de senos y cosenos, y cada ŵ_i se divide por λ_i, que también depende de θ. Así que, a diferencia de los mínimos cuadrados, el mínimo de J no tiene fórmula. Lo encuentra una búsqueda iterativa: parte de un valor inicial θ₀ y sigue cambiando θ mientras J baje. Algunas búsquedas, como el símplex de Nelder–Mead, solo comparan valores de J; otras siguen su pendiente. En ambos casos, una búsqueda encuentra el punto más bajo cerca de donde parte, que puede no ser el más bajo de todos, así que el valor inicial importa. La sección 3 da uno bueno: los parámetros de la A lineal.

¿Para qué buscar, si la A lineal ya ajusta? Por tres razones.

• J mide lo que importa, distancias en la imagen, con todos los puntos contando lo mismo.
• El modelo puede crecer. La distorsión del lente, en la sección siguiente, o un manipulador que gira el objeto entre imágenes, agregan parámetros que ninguna matriz A puede contener.
• Se puede incorporar lo que se sabe de la cámara. Para una cámara que se sabe sin torcimiento, γ se puede fijar en 0, algo que las doce entradas libres de A no pueden expresar.

Un promedio de distancias, o de cuadrados

J promedia distancias, no distancias al cuadrado, y las dos llevan a cámaras distintas. Para la mesa, el menor J es 2.69 píxeles, y se alcanza con cinco de las ocho esquinas ajustadas exactamente y las otras tres hasta a 10.6 píxeles. Minimizar en cambio las distancias al cuadrado reparte el error entre las ocho esquinas, ninguna a más de 5 píxeles, con una distancia promedio de 3.40 píxeles. Un promedio de distancias se comporta como una mediana: tolera unos pocos errores grandes si el resto ajusta exactamente, así que un clic malo lo arrastra menos. Un cuadrado hace que un error grande cueste más que varios pequeños, y la mayoría de los programas de calibración minimizan cuadrados e informan su raíz cuadrática media.

El promedio de distancias tiene además pliegues agudos, donde una distancia llega a 0 y su pendiente salta. Pueden detener una búsqueda: para la mesa, el símplex de Nelder–Mead se detuvo en cualquier punto entre 2.9 y 3.3 píxeles, según dónde partiera. Resolver una sucesión de problemas de mínimos cuadrados ponderados, cada uno pesando el error al cuadrado de un punto por uno sobre su distancia actual, llegó a 2.69 desde cada uno de varios puntos de partida.

Una imagen especular que ajusta

Una cámara equivocada ajusta los píxeles exactamente igual de bien que la correcta. Gira la cámara media vuelta en torno a su propio eje Y, de modo que mire en sentido contrario al objeto, y cambia los signos de β y γ. En las coordenadas de la cámara, el giro cambia los signos de X y Z de cada punto. Dividir −X por −Z no cambia nada, pero ahora Y se divide por −Z, y los nuevos signos de β y γ, los dos parámetros que multiplican a Y, lo compensan. La matriz de la nueva cámara es exactamente −A, y −A da la misma imagen que A. Pero ahora cada λ_i es negativo: todos los puntos están detrás de la cámara. Una búsqueda que parte de un valor inicial así, con un β negativo, se queda con él. Por eso, después de calibrar, revisa que α y β sean positivos y que cada λ_i sea positivo.

5Distorsión del lente

El modelo hasta aquí supone que los rayos pasan en línea recta por el centro óptico. Un lente real los curva un poco, sobre todo cerca de los bordes de la imagen, donde su superficie está más inclinada respecto del eje óptico. Las rectas de la escena se ven entonces curvas en la imagen.

La distorsión se modela como un paso más entre la proyección y la conversión a píxeles. P da el punto de la imagen ideal (x, y), una función de distorsión lo mueve al punto de la imagen real (x′, y′), el que registra el sensor, y K lo convierte a píxeles:

El modelo completo tiene cuatro pasos: S′ lleva el punto a las coordenadas de la cámara, P lo proyecta sobre la imagen ideal, la distorsión lo mueve a la imagen real, y K da su píxel. Solo la imagen real se ve.

Radial y tangencial

La distorsión suele separarse en dos partes, en torno a un centro donde se anula, casi siempre el punto principal. En coordenadas polares (r, θ) en torno a ese centro, una parte radial δ_r(r) mueve un punto a lo largo de su rayo desde el centro y depende solo de su distancia r. Una parte tangencial δ_θ(θ) lo gira en torno al centro y depende solo de su ángulo θ. De vuelta en x e y, con el centro en el origen:

Los primeros términos mueven el punto en δ_r en la dirección (cos θ, sin θ); los segundos lo mueven en un arco r δ_θ en la dirección perpendicular (−sin θ, cos θ). Ambas partes suelen ser polinomios de orden mayor que uno. Un δ_r proporcional a r solo cambiaría la escala de toda la imagen, lo que α y β ya hacen, y un δ_θ constante solo la giraría, lo que los ángulos ya hacen.

Una distorsión radial que empuja los puntos hacia afuera cada vez más hacia los bordes curva hacia adentro las rectas de una cuadrícula, como un cojín. Una que los tira hacia adentro las curva hacia afuera, como un barril, algo común en lentes gran angular.

Un ejemplo

Un lente forma la imagen de nueve puntos en una cuadrícula, con x e y en {−1, 0, 1}. El centro y los cuatro puntos medios de los lados quedan donde estaban, pero las cuatro esquinas pasan de (±1, ±1) a (±1.2, ±1.2). Cada punto queda sobre su propio rayo desde el centro, así que no hay distorsión tangencial: δ_θ = 0. La distancia al centro cambia solo en las esquinas, de √2 a 1.2√2, así que δ_r(1) = 0 y δ_r(√2) = 0.2√2. Un polinomio que se anula en r = 0 y r = 1 es δ_r(r) = k(r³ − r), y en r = √2 da k√2, así que k = 0.2. Como cos θ = x/r y sin θ = y/r,

Comprueba una esquina: x′ = 1 + 0.2 · 1 · (1 + 1 − 1) = 1.2. El punto medio de un lado, (1, 0), tiene x² + y² − 1 = 0 y no se mueve. Otros polinomios también ajustan los nueve puntos. Este tiene solo potencias impares de r, lo que mantiene δ_x y δ_y como polinomios en x e y.

Por qué solo el método no lineal puede incluirla

Una distorsión curva las rectas, y ninguna matriz puede hacer eso: A envía cada recta a una recta. Así que una calibración lineal no puede modelar la distorsión. Una no lineal agrega sus coeficientes a θ y los estima junto con el resto.

6Valores atípicos y RANSAC

Los mínimos cuadrados confían en todos los puntos. Si una esquina se marca en el lugar equivocado, por ejemplo la esquina 6, 150 píxeles a la izquierda, el ajuste se inclina hacia ella y todas las esquinas lo sufren. Los errores crecen de 2–5 píxeles a 9–67 píxeles, J crece de 3.69 a 35.7 píxeles, y la cámara encontrada se mueve 22 cm. La esquina equivocada ni siquiera se destaca: su error, 56 píxeles, es menor que el de la esquina 5, que se marcó bien. Un punto equivocado en los datos se llama valor atípico (outlier).

RANSAC, por random sample consensus (consenso de muestras aleatorias), da vuelta la pregunta. En vez de ajustar todos los puntos a la vez, repite cuatro pasos:

1. Elige una muestra al azar con tan pocos puntos como necesite el modelo: seis esquinas para una cámara.
2. Ajusta A a la muestra.
3. Cuenta cuántos de todos los puntos reproyecta esta A dentro de una tolerancia de sus píxeles. Estos son sus puntos consistentes (inliers).
4. Quédate con la muestra con más puntos consistentes vista hasta ahora.

Tras suficientes intentos, es probable que alguna muestra no tenga valores atípicos, y su A concuerda con todos los puntos buenos, lo que le da la mayor cantidad de puntos consistentes. Un ajuste final por mínimos cuadrados sobre esos puntos da la respuesta. La muestra se mantiene lo más pequeña posible porque una muestra pequeña tiene más probabilidad de estar libre de valores atípicos, así que se necesitan menos intentos.

El mismo método estimó homografías en los apuntes de transformaciones 2D, con muestras de cuatro pares de puntos.

RANSAC con solo ocho esquinas

Ocho esquinas son pocas para un modelo con once parámetros. Seis esquinas dan doce ecuaciones para once incógnitas, así que A ajusta muy de cerca casi cualquier grupo de seis esquinas, y una muestra al azar dice poco sobre el resto. Con tan pocas esquinas, en cambio, no hace falta muestrear: la figura de la sección 3 ajusta A a cada conjunto de ocho, siete o seis esquinas, 37 conjuntos en total, y se queda con el conjunto más grande cuyas esquinas caen todas dentro de la tolerancia de su propio ajuste. Ese conjunto es el consenso que RANSAC busca con muestras al azar. Cambia la figura a RANSAC y arrastra una esquina lejos. Para la esquina 6 marcada 150 píxeles fuera de lugar, el conjunto más grande son las otras siete esquinas, su ajuste reproyecta cada una a menos de 3.4 píxeles, y la esquina 6 queda fuera como valor atípico, a 162 píxeles de donde la pone el ajuste.

La tolerancia es de 20 píxeles, porque una cámara ajustada a siete de las esquinas predice la octava con unos 7 a 24 píxeles de error aunque todos los clics sean correctos, y una tolerancia más estrecha declararía valores atípicos a esquinas buenas. Además, una esquina arrastrada solo un poco a veces se explica igual de bien dejando fuera una esquina vecina y desplazando el punto principal. Así que, cuando califican dos conjuntos del mismo tamaño, la figura prefiere la cámara con el punto principal más cerca del centro de la imagen y el menor torcimiento. Arrastra la esquina más lejos, y es ella la que queda fuera.

Marcar puntos a mano es lento. Los programas encuentran puntos correspondientes por sí solos, detectando las esquinas de un tablero de ajedrez o emparejando puntos clave entre imágenes, y algunas de esas parejas son erróneas. RANSAC es lo que hace eso práctico. Entre dos vistas de una misma escena, la restricción epipolar da una prueba rápida para una pareja propuesta.

7La calibración en la práctica

Los objetos de calibración suelen ser planos: un tablero de ajedrez impreso, cuyas esquinas interiores un programa encuentra por sí solo con precisión de una fracción de píxel. Por la sección 2, una vista de un plano no basta, porque solo determina una homografía. Pero cada vista del tablero, sostenido en una pose distinta, da otra homografía, y todas comparten la misma K. Unas pocas vistas en poses distintas determinan los cinco parámetros intrínsecos, y luego el R′ y el t′ de cada vista.

OpenCV hace todo esto en una sola llamada, cv.calibrateCamera. Recibe las esquinas del tablero en las coordenadas del propio tablero, con Z′ = 0, y sus píxeles en cada vista. Devuelve K, los coeficientes de su modelo de distorsión del lente y, para cada vista, R′ y t′. OpenCV toma el sistema de la cámara como el nuevo, la forma inversa descrita en la sección 3, y guarda R′ como un vector de rotación: un eje cuyo largo es el ángulo girado en torno a él. La matriz A de cada vista es entonces [K 0][R′ t′; 0^T 1].

Los sistemas de imágenes de rayos X se calibran de la misma forma, a partir de imágenes de un tablero de calibración. La fuente de rayos X hace el papel del centro óptico, como mostraron los apuntes de transformaciones 3D.

8Resumen

Las ideas clave junto a sus equivalentes en NumPy. En código, los puntos 3D son las columnas de M, sus píxeles las columnas de w, y A es P.