math - corrigiendo la distorsión fisheye programáticamente
graphics geometry (6)
(Cartel original, proporcionando una alternativa)
La siguiente función asigna coordenadas de destino (rectilíneas) a coordenadas fuente (distorsionadas por ojo de pez). (Agradecería ayuda para revertirlo)
Llegué a este punto a través de prueba y error: no comprendo fundamentalmente por qué este código funciona, ¡se agradecen las explicaciones y la precisión mejorada !
def dist(x,y):
return sqrt(x*x+y*y)
def correct_fisheye(src_size,dest_size,dx,dy,factor):
""" returns a tuple of source coordinates (sx,sy)
(note: values can be out of range)"""
# convert dx,dy to relative coordinates
rx, ry = dx-(dest_size[0]/2), dy-(dest_size[1]/2)
# calc theta
r = dist(rx,ry)/(dist(src_size[0],src_size[1])/factor)
if 0==r:
theta = 1.0
else:
theta = atan(r)/r
# back to absolute coordinates
sx, sy = (src_size[0]/2)+theta*rx, (src_size[1]/2)+theta*ry
# done
return (int(round(sx)),int(round(sy)))
Cuando se usa con un factor de 3.0, distorsiona con éxito las imágenes utilizadas como ejemplos (no hice ningún intento de interpolación de calidad):
Enlace muerto
(Y esto es de la publicación del blog, para comparar :)
ACTUALIZACIÓN DE ESTADO DE BOUNTY:
Descubrí cómo mapear una lente lineal , desde las coordenadas de destination hasta las coordenadas de source .
¿Cómo se calcula la distancia radial desde el centro para pasar de ojo de pez a rectilíneo?
1). De hecho, me cuesta revertirlo y asignar las coordenadas de origen a las coordenadas de destino. ¿Qué es lo inverso, en el código en el estilo de las funciones de conversión que publiqué?
2). También veo que mi distorsión es imperfecta en algunas lentes, presumiblemente aquellas que no son estrictamente lineales. ¿Cuál es el equivalente de las coordenadas origen-destino para esas lentes? De nuevo, más código que solo fórmulas matemáticas, por favor ...
Pregunta como se dijo originalmente:
Tengo algunos puntos que describen las posiciones en una imagen tomada con una lente ojo de pez.
Quiero convertir estos puntos en coordenadas rectilíneas. Quiero distorsionar la imagen.
Encontré esta descripción de cómo generar un efecto ojo de pez, pero no cómo revertirlo.
También hay una publicación en el blog que describe cómo usar herramientas para hacerlo; estas fotos son de eso:
(1) : SOURCE Enlace original de la foto
Entrada: imagen original con distorsión ojo de pez para corregir.
(2) : DESTINATION Enlace original de la foto
Salida: imagen corregida (técnicamente también con corrección de perspectiva, pero eso es un paso separado).
¿Cómo se calcula la distancia radial desde el centro para pasar de ojo de pez a rectilíneo?
El talón de mi función se ve así:
Point correct_fisheye(const Point& p,const Size& img) {
// to polar
const Point centre = {img.width/2,img.height/2};
const Point rel = {p.x-centre.x,p.y-centre.y};
const double theta = atan2(rel.y,rel.x);
double R = sqrt((rel.x*rel.x)+(rel.y*rel.y));
// fisheye undistortion in here please
//... change R ...
// back to rectangular
const Point ret = Point(centre.x+R*cos(theta),centre.y+R*sin(theta));
fprintf(stderr,"(%d,%d) in (%d,%d) = %f,%f = (%d,%d)/n",p.x,p.y,img.width,img.height,theta,R,ret.x,ret.y);
return ret;
}
Alternativamente, podría de alguna manera convertir la imagen de ojo de pez a rectilínea antes de encontrar los puntos, pero estoy completamente aturdido por la documentación de OpenCV . ¿Existe una manera directa de hacerlo en OpenCV, y funciona lo suficientemente bien como para hacerlo en un video en vivo?
Encontré este archivo pdf y probé que las matemáticas son correctas (a excepción de la línea vd = *xd**fv+v0 which should say vd = **yd**+fv+v0 ).
http://perception.inrialpes.fr/CAVA_Dataset/Site/files/Calibration_OpenCV.pdf
No utiliza todos los últimos coeficientes que OpenCV tiene disponibles, pero estoy seguro de que se podría adaptar con bastante facilidad.
double k1 = cameraIntrinsic.distortion[0];
double k2 = cameraIntrinsic.distortion[1];
double p1 = cameraIntrinsic.distortion[2];
double p2 = cameraIntrinsic.distortion[3];
double k3 = cameraIntrinsic.distortion[4];
double fu = cameraIntrinsic.focalLength[0];
double fv = cameraIntrinsic.focalLength[1];
double u0 = cameraIntrinsic.principalPoint[0];
double v0 = cameraIntrinsic.principalPoint[1];
double u, v;
u = thisPoint->x; // the undistorted point
v = thisPoint->y;
double x = ( u - u0 )/fu;
double y = ( v - v0 )/fv;
double r2 = (x*x) + (y*y);
double r4 = r2*r2;
double cDist = 1 + (k1*r2) + (k2*r4);
double xr = x*cDist;
double yr = y*cDist;
double a1 = 2*x*y;
double a2 = r2 + (2*(x*x));
double a3 = r2 + (2*(y*y));
double dx = (a1*p1) + (a2*p2);
double dy = (a3*p1) + (a1*p2);
double xd = xr + dx;
double yd = yr + dy;
double ud = (xd*fu) + u0;
double vd = (yd*fv) + v0;
thisPoint->x = ud; // the distorted point
thisPoint->y = vd;
Implementé ciegamente las fórmulas desde aquí , así que no puedo garantizar que haga lo que necesita.
Use auto_zoom para obtener el valor del parámetro de zoom .
def dist(x,y):
return sqrt(x*x+y*y)
def fisheye_to_rectilinear(src_size,dest_size,sx,sy,crop_factor,zoom):
""" returns a tuple of dest coordinates (dx,dy)
(note: values can be out of range)
crop_factor is ratio of sphere diameter to diagonal of the source image"""
# convert sx,sy to relative coordinates
rx, ry = sx-(src_size[0]/2), sy-(src_size[1]/2)
r = dist(rx,ry)
# focal distance = radius of the sphere
pi = 3.1415926535
f = dist(src_size[0],src_size[1])*factor/pi
# calc theta 1) linear mapping (older Nikon)
theta = r / f
# calc theta 2) nonlinear mapping
# theta = asin ( r / ( 2 * f ) ) * 2
# calc new radius
nr = tan(theta) * zoom
# back to absolute coordinates
dx, dy = (dest_size[0]/2)+rx/r*nr, (dest_size[1]/2)+ry/r*nr
# done
return (int(round(dx)),int(round(dy)))
def fisheye_auto_zoom(src_size,dest_size,crop_factor):
""" calculate zoom such that left edge of source image matches left edge of dest image """
# Try to see what happens with zoom=1
dx, dy = fisheye_to_rectilinear(src_size, dest_size, 0, src_size[1]/2, crop_factor, 1)
# Calculate zoom so the result is what we wanted
obtained_r = dest_size[0]/2 - dx
required_r = dest_size[0]/2
zoom = required_r / obtained_r
return zoom
La descripción que mencionas indica que la proyección mediante una cámara de agujero de alfiler (una que no introduce distorsión de la lente) está modelada por
R_u = f*tan(theta)
y la proyección con cámaras de ojo de ojo de pez comunes (es decir, distorsionadas) está modelada por
R_d = 2*f*sin(theta/2)
Ya conoce R_d y theta y si conoce la longitud focal de la cámara (representada por f), corregir la imagen equivaldría a calcular R_u en términos de R_d y theta. En otras palabras,
R_u = f*tan(2*asin(R_d/(2*f)))
es la fórmula que estás buscando. La estimación de la distancia focal f puede resolverse calibrando la cámara u otros medios, como permitir que el usuario proporcione comentarios sobre la corrección de la imagen o el conocimiento de la escena original.
Para resolver el mismo problema usando OpenCV, debería obtener los parámetros intrínsecos de la cámara y los coeficientes de distorsión de la lente. Véase, por ejemplo, el Capítulo 11 de Learning OpenCV (no olvide verificar la correction ). Luego puede usar un programa como este (escrito con los enlaces de Python para OpenCV) para invertir la distorsión de la lente:
#!/usr/bin/python
# ./undistort 0_0000.jpg 1367.451167 1367.451167 0 0 -0.246065 0.193617 -0.002004 -0.002056
import sys
import cv
def main(argv):
if len(argv) < 10:
print ''Usage: %s input-file fx fy cx cy k1 k2 p1 p2 output-file'' % argv[0]
sys.exit(-1)
src = argv[1]
fx, fy, cx, cy, k1, k2, p1, p2, output = argv[2:]
intrinsics = cv.CreateMat(3, 3, cv.CV_64FC1)
cv.Zero(intrinsics)
intrinsics[0, 0] = float(fx)
intrinsics[1, 1] = float(fy)
intrinsics[2, 2] = 1.0
intrinsics[0, 2] = float(cx)
intrinsics[1, 2] = float(cy)
dist_coeffs = cv.CreateMat(1, 4, cv.CV_64FC1)
cv.Zero(dist_coeffs)
dist_coeffs[0, 0] = float(k1)
dist_coeffs[0, 1] = float(k2)
dist_coeffs[0, 2] = float(p1)
dist_coeffs[0, 3] = float(p2)
src = cv.LoadImage(src)
dst = cv.CreateImage(cv.GetSize(src), src.depth, src.nChannels)
mapx = cv.CreateImage(cv.GetSize(src), cv.IPL_DEPTH_32F, 1)
mapy = cv.CreateImage(cv.GetSize(src), cv.IPL_DEPTH_32F, 1)
cv.InitUndistortMap(intrinsics, dist_coeffs, mapx, mapy)
cv.Remap(src, dst, mapx, mapy, cv.CV_INTER_LINEAR + cv.CV_WARP_FILL_OUTLIERS, cv.ScalarAll(0))
# cv.Undistort2(src, dst, intrinsics, dist_coeffs)
cv.SaveImage(output, dst)
if __name__ == ''__main__'':
main(sys.argv)
También tenga en cuenta que OpenCV usa un modelo de distorsión de lente muy diferente al de la página web a la que se vinculó.
Si crees que tus fórmulas son exactas, puedes computir una fórmula exacta con trigonometría, así:
Rin = 2 f sin(w/2) -> sin(w/2)= Rin/2f
Rout= f tan(w) -> tan(w)= Rout/f
(Rin/2f)^2 = [sin(w/2)]^2 = (1 - cos(w))/2 -> cos(w) = 1 - 2(Rin/2f)^2
(Rout/f)^2 = [tan(w)]^2 = 1/[cos(w)]^2 - 1
-> (Rout/f)^2 = 1/(1-2[Rin/2f]^2)^2 - 1
Sin embargo, como dice @jmbr, la distorsión real de la cámara dependerá de la lente y el zoom. En lugar de confiar en una fórmula fija, es posible que desee intentar una expansión polinómica:
Rout = Rin*(1 + A*Rin^2 + B*Rin^4 + ...)
Ajustando primero A, luego coeficientes de orden superior, puede calcular cualquier función local razonable (la forma de la expansión aprovecha la simetría del problema). En particular, debería ser posible calcular los coeficientes iniciales para aproximar la función teórica anterior.
Además, para obtener buenos resultados, deberá usar un filtro de interpolación para generar su imagen corregida. Siempre que la distorsión no sea demasiado grande, puede usar el tipo de filtro que usaría para reescalar la imagen linealmente sin mucho problema.
Editar: según su solicitud, el factor de escala equivalente para la fórmula anterior:
(Rout/f)^2 = 1/(1-2[Rin/2f]^2)^2 - 1
-> Rout/f = [Rin/f] * sqrt(1-[Rin/f]^2/4)/(1-[Rin/f]^2/2)
Si trazas la fórmula anterior junto con el bronceado (Rin / f), puedes ver que tienen una forma muy similar. Básicamente, la distorsión de la tangente se vuelve severa antes de que el pecado (w) sea muy diferente de w.
La fórmula inversa debería ser algo así como:
Rin/f = [Rout/f] / sqrt( sqrt(([Rout/f]^2+1) * (sqrt([Rout/f]^2+1) + 1) / 2 )
Tomé lo que hizo JMBR y, básicamente, lo revirtió. Tomó el radio de la imagen distorsionada (Rd, es decir, la distancia en píxeles desde el centro de la imagen) y encontró una fórmula para Ru, el radio de la imagen sin distorsión.
Quieres ir por el otro lado. Para cada píxel en la imagen sin distorsión (procesada), desea saber cuál es el píxel correspondiente en la imagen distorsionada. En otras palabras, dado (xu, yu) -> (xd, yd). A continuación, reemplaza cada píxel de la imagen sin distorsión con su píxel correspondiente de la imagen distorsionada.
Comenzando donde lo hizo JMBR, hago lo contrario, encontrando Rd como una función de Ru. Yo obtengo:
Rd = f * sqrt(2) * sqrt( 1 - 1/sqrt(r^2 +1))
donde f es la distancia focal en píxeles (lo explicaré más adelante), r = Ru/f .
La distancia focal para mi cámara fue de 2.5 mm. El tamaño de cada píxel en mi CCD era de 6 um cuadrado. f fue, por lo tanto, 2500/6 = 417 píxeles. Esto se puede encontrar por prueba y error.
Finding Rd le permite encontrar el píxel correspondiente en la imagen distorsionada usando coordenadas polares.
El ángulo de cada píxel desde el punto central es el mismo:
theta = arctan( (yu-yc)/(xu-xc) ) donde xc, yc son los puntos centrales.
Entonces,
xd = Rd * cos(theta) + xc
yd = Rd * sin(theta) + yc
Asegúrese de saber en qué cuadrante se encuentra.
Aquí está el código C # que utilicé
public class Analyzer
{
private ArrayList mFisheyeCorrect;
private int mFELimit = 1500;
private double mScaleFESize = 0.9;
public Analyzer()
{
//A lookup table so we don''t have to calculate Rdistorted over and over
//The values will be multiplied by focal length in pixels to
//get the Rdistorted
mFisheyeCorrect = new ArrayList(mFELimit);
//i corresponds to Rundist/focalLengthInPixels * 1000 (to get integers)
for (int i = 0; i < mFELimit; i++)
{
double result = Math.Sqrt(1 - 1 / Math.Sqrt(1.0 + (double)i * i / 1000000.0)) * 1.4142136;
mFisheyeCorrect.Add(result);
}
}
public Bitmap RemoveFisheye(ref Bitmap aImage, double aFocalLinPixels)
{
Bitmap correctedImage = new Bitmap(aImage.Width, aImage.Height);
//The center points of the image
double xc = aImage.Width / 2.0;
double yc = aImage.Height / 2.0;
Boolean xpos, ypos;
//Move through the pixels in the corrected image;
//set to corresponding pixels in distorted image
for (int i = 0; i < correctedImage.Width; i++)
{
for (int j = 0; j < correctedImage.Height; j++)
{
//which quadrant are we in?
xpos = i > xc;
ypos = j > yc;
//Find the distance from the center
double xdif = i-xc;
double ydif = j-yc;
//The distance squared
double Rusquare = xdif * xdif + ydif * ydif;
//the angle from the center
double theta = Math.Atan2(ydif, xdif);
//find index for lookup table
int index = (int)(Math.Sqrt(Rusquare) / aFocalLinPixels * 1000);
if (index >= mFELimit) index = mFELimit - 1;
//calculated Rdistorted
double Rd = aFocalLinPixels * (double)mFisheyeCorrect[index]
/mScaleFESize;
//calculate x and y distances
double xdelta = Math.Abs(Rd*Math.Cos(theta));
double ydelta = Math.Abs(Rd * Math.Sin(theta));
//convert to pixel coordinates
int xd = (int)(xc + (xpos ? xdelta : -xdelta));
int yd = (int)(yc + (ypos ? ydelta : -ydelta));
xd = Math.Max(0, Math.Min(xd, aImage.Width-1));
yd = Math.Max(0, Math.Min(yd, aImage.Height-1));
//set the corrected pixel value from the distorted image
correctedImage.SetPixel(i, j, aImage.GetPixel(xd, yd));
}
}
return correctedImage;
}
}