correction programmatique de la distorsion fisheye

MISE À JOUR DU STATUT DE BOUNTY:

j'ai découvert comment cartographier une lentille linéaire , des coordonnées destination aux coordonnées source .

comment calculer la distance radiale à partir du centre pour passer de fisheye à rectiligne?

  • 1). j'ai du mal à l'inverser, et à cartographier les coordonnées de la source aux coordonnées de destination. Quel est l'inverse, en code dans le style des fonctions de conversion que j'ai posté?

  • 2). je vois aussi que mon undistortion est imparfaite sur certaines lentilles - probablement celles qui ne sont pas strictement linéaires. Quel est l'équivalent des coordonnées de source et de destination pour ces lentilles? Encore une fois, plus de code que de simples formules mathématiques s'il vous plaît...


Question initialement déclaré:

j'ai quelques points qui décrivent les positions dans une image prise avec un objectif fisheye.

je veux convertir ces points en coordonnées rectilignes. Je veux undistort l'image.

j'ai trouvé cette description de la façon de générer un effet fisheye, mais pas comment inverser la tendance.

il y a aussi un post de blog qui décrit comment utiliser les outils pour le faire; ces photos sont à partir de cela:

(1) : SOURCE lien photo Original



Entrée: image originale avec distorsion fish-eye à corriger.

(2) : DESTINATION lien photo Original



Sortie: image corrigée (techniquement aussi avec correction de perspective, mais c'est une étape séparée).

comment calculer la distance radiale à partir du centre pour passer de fisheye à rectiligne?

mon talon de fonction ressemble à ceci:

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;
}

alternativement, je pourrais d'une façon ou d'une autre convertir l'image de fisheye en rectiligne avant de trouver les points, mais je suis complètement embrouillé par la documentation OpenCV . Est-il un moyen simple de le faire dans OpenCV, et le fait d'effectuer assez bien pour le faire à un flux vidéo en direct?

48
demandé sur Community 2010-03-19 16:50:13

6 réponses

la description que vous mentionnez indique que la projection par une caméra à trous (qui n'introduit pas de distorsion de la lentille) est modélisée par

R_u = f*tan(theta)

et la projection par des caméras fisheye ordinaires (c'est-à-dire déformées) est modélisée par

R_d = 2*f*sin(theta/2)

vous connaissez déjà r_d et theta et si vous saviez La Distance focale de la caméra (représentée par f) alors corriger l'image reviendrait à calculer R_u in termes de R_d et thêta. En d'autres termes,

R_u = f*tan(2*asin(R_d/(2*f)))

est la formule que vous cherchez. Estimer la distance focale f peut être résolu en calibrant la caméra ou d'autres moyens tels que laisser l'Utilisateur fournir des commentaires sur la façon dont l'image est corrigée ou en utilisant la connaissance de la scène originale.

pour résoudre le même problème avec OpenCV, il faut obtenir les paramètres intrinsèques de la caméra et les coefficients de distorsion de la lentille. Voir, par exemple, le Chapitre 11 de Learning OpenCV (n'oubliez pas de cocher la correction ). Vous pouvez alors utiliser un programme comme celui-ci (écrit avec les fixations Python pour OpenCV) afin d'inverser la distorsion de l'objectif:

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

notez également Qu'OpenCV utilise un modèle de distorsion optique très différent de celui de la page Web à laquelle vous avez fait référence.

30
répondu jmbr 2010-03-21 21:55:44

(affiche originale, offrant une alternative)

la fonction suivante établit une correspondance entre les coordonnées de destination (rectilignes) et les coordonnées de source (déformées par fisheye). (je vous en serais reconnaissant de l'aide dans l'inversion de cette tendance)

j'en suis arrivé à ce point par tâtonnements: Je ne comprends pas fondamentalement pourquoi ce code fonctionne, des explications et une meilleure précision sont appréciées !

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

lorsqu'il est utilisé avec un facteur de 3.0, il efface avec succès les images utilisées comme exemples (je n'ai fait aucune tentative d'interpolation de qualité):

lien Mort

(et ceci est de l'article de blog, pour comparaison:)

Using Panotools

7
répondu Will 2018-01-19 18:14:07

si vous pensez que vos formules sont exactes, Vous pouvez calculer une formule exacte avec trig, comme cela:

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

cependant, comme le dit @jmbr, la distorsion réelle de la caméra dépendra de l'objectif et du zoom. Plutôt que de s'appuyer sur une formule fixe, vous pourriez essayer une expansion polynomiale:

Rout = Rin*(1 + A*Rin^2 + B*Rin^4 + ...)

en ajustant d'abord A, puis les coefficients d'ordre supérieur, vous pouvez calculer n'importe quelle fonction locale raisonnable (la forme de l'expansion prend l'avantage de la symétrie du problème). En particulier, il devrait être possible de calculer des coefficients initiaux pour rapprocher la fonction théorique ci-dessus.

Aussi, pour obtenir de bons résultats, vous aurez besoin d'utiliser un filtre d'interpolation pour générer votre image corrigée. Tant que la distorsion n'est pas trop grande, vous pouvez utiliser le type de filtre que vous utiliseriez pour modifier l'image de façon linéaire sans trop de problèmes.

Edit: selon votre demande, la facteur d'échelle équivalent pour la formule ci-dessus:

(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 vous tracez la formule ci-dessus à côté de tan(Rin/f), vous pouvez voir qu'ils sont très similaires dans la forme. Fondamentalement, la distorsion de la tangente devient sévère avant que le péché(w) devienne très différent de W.

la formule inverse devrait être quelque chose comme:

Rin/f = [Rout/f] / sqrt( sqrt(([Rout/f]^2+1) * (sqrt([Rout/f]^2+1) + 1) / 2 )
3
répondu comingstorm 2010-03-26 08:02:24

j'ai pris ce que JMBR a fait et je l'ai inversé. Il a pris le rayon de l'image déformée (Rd, c'est-à-dire la distance en pixels du centre de l'image) et a trouvé une formule pour Ru, le rayon de l'image non déformée.

Vous voulez aller dans l'autre sens. Pour chaque pixel dans l'image non déformée (image traitée), vous voulez savoir quel est le pixel correspondant dans l'image déformée. En d'autres termes, donné (xu, yu) --> (xd, yd). Vous remplacez ensuite chaque pixel dans le une image pixel correspondant de l'image déformée.

à partir de là où JMBR a fait, je fais l'inverse, en trouvant Rd en fonction de Ru. Je comprends:

Rd = f * sqrt(2) * sqrt( 1 - 1/sqrt(r^2 +1))

où f est la distance focale en pixels (Je l'expliquerai plus tard), et r = Ru/f .

la distance focale de mon appareil était de 2,5 mm. La taille de chaque pixel sur mon CCD était de 6 um carré. f était donc de 2500/6 = 417 pixels. Cela peut être trouvé par trial et erreur.

Finding Rd vous permet de trouver le pixel correspondant dans l'image déformée en utilisant des coordonnées polaires.

l'angle de chaque pixel par rapport au point central est le même:

theta = arctan( (yu-yc)/(xu-xc) ) où xc, yc sont les points centraux.

Puis,

xd = Rd * cos(theta) + xc
yd = Rd * sin(theta) + yc

assurez-vous de savoir dans quel quadrant vous êtes.

voici le code C que j'ai utilisé

 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;
      }
}
3
répondu Barry Vant-Hull 2013-06-07 00:10:58

j'ai trouvé ce fichier pdf, et j'ai prouvé que les mathématiques sont corrects (sauf pour la ligne vd = *xd**fv+v0 which should say vd = **yd**+fv+v0 ).

http://perception.inrialpes.fr/CAVA_Dataset/Site/files/Calibration_OpenCV.pdf

Il n'utilise pas tous les derniers coefficients OpenCV est disponible, mais je suis sûr qu'il pourrait être adapté assez facilement.

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;
3
répondu Myke Smyth 2017-03-17 04:45:05

j'ai aveuglément mis en œuvre les formules de ici , donc je ne peux pas garantir qu'il ferait ce dont vous avez besoin.

utilisez auto_zoom pour obtenir la valeur du paramètre 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
2
répondu Roman Zenka 2010-03-26 14:47:18