【问题标题】:correcting fisheye distortion programmatically以编程方式校正鱼眼失真
【发布时间】:2011-01-29 11:37:35
【问题描述】:

赏金状态更新:

I discovered how to map a linear lens,从destination 坐标到source 坐标。

如何计算从中心到从鱼眼到直线的径向距离?

  • 1). 我实际上很难将其反转,并将源坐标映射到目标坐标。在我发布的转换函数样式的代码中,相反的是什么?

  • 2). 我还发现我的不失真在某些镜头上并不完美——大概是那些不是严格线性的镜头。这些镜头的等效往返源和目的地坐标是多少?再次,请提供比数学公式更多的代码...


问题如前所述:

我有一些点可以描述用鱼眼镜头拍摄的照片中的位置。

我想将这些点转换为直线坐标。我想不扭曲图像。

我找到了this description 如何生成鱼眼效果,但没有找到如何反转它。

还有一个blog post 描述了如何使用工具来做到这一点;这些图片来自于:

(1) : SOURCE Original photo link

输入:需要修复鱼眼失真的原始图像。

(2) : DESTINATION Original photo link

输出:校正后的图像(技术上也带有透视校正,但这是一个单独的步骤)。

如何计算从中心到从鱼眼到直线的径向距离?

我的函数存根如下所示:

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

或者,我可以在找到点之前以某种方式将图像从鱼眼转换为直线,但我完全被 OpenCV documentation 弄糊涂了。在 OpenCV 中是否有一种直接的方法可以做到这一点,并且它的性能是否足以在实时视频源中做到这一点?

【问题讨论】:

  • 我不太明白你在找什么。鱼眼从球体映射到图片平面。反向映射将从图片回到球体,对吗?你在找什么直线坐标?
  • @mtrw 我的源图像是鱼眼扭曲的,我想不扭曲它
  • 那么photo.net/learn/fisheye上的图片是你要找的吗?
  • 是的,正确的图片例如通过 OpenCV,或用于校正图片中任何点的公式。
  • 威尔,你有没有得到一个确凿的答案?我很想看看你最终得到的任何代码。

标签: math graphics geometry projection


【解决方案1】:

(原始海报,提供替代方案)

以下函数将目标(直线)坐标映射到源(鱼眼扭曲)坐标。 (我希望能帮助您扭转它)

我通过反复试验得到了这一点:我从根本上不明白为什么这段代码可以工作,感谢解释和提高准确性

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

当与 3.0 的因子一起使用时,它成功地使用作示例的图像不失真(我没有尝试质量插值):

死链接

(这是来自博客文章,用于比较:)

【讨论】:

  • 您的代码工作的原因是您将 (rx,ry) 缩放了一个因子 theta(现在是一个比率,而不是一个角度)。如果原始镜头具有(如维基文章所述)从视角到图像偏移的“线性映射”,我相信 atan(r)/r 是正确的。
  • 映射的反向应该是按 tan(r')/r' 的因子进行缩放,其中 r' 是未失真图像中心的半径。
  • 这两项工作的原因是,如果你有向量 v' = v*k(|v|),并且你想要 |v'|=f(|v|),你取第一个方程的绝对值: |v'|=|v|*k(|v|)=f(|v|) -- 使得 k(|v|)=f(|v|)/| v|
  • @comingstorm 那么非线性映射的等价物是什么?
  • @stkent,所以我在 R 中执行此操作及其工作,除了根据我选择的因素,我看到“波浪”叠加在我正在拉伸或未拉伸的图像上。例如,我得到了所需的桶形失真,但是在它上面有某种看起来像浅灰色波的假象。想知道这是某种边界问题还是相移类型的事情?示例:imgur.com/a/lBNXPYh
【解决方案2】:

我找到了这个 pdf 文件,并且我已经证明数学是正确的(vd = *xd**fv+v0 which should say vd = **yd**+fv+v0 行除外)。

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

它没有使用 OpenCV 可用的所有最新系数,但我相信它可以相当容易地适应。

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】:

    我接受了 JMBR 的做法,并且基本上扭转了它。他取了畸变图像的半径(Rd,即到图像中心的像素距离),并找到了一个公式,即未畸变图像的半径Ru。

    你想走另一条路。对于未失真(处理后的图像)中的每个像素,您想知道失真图像中对应的像素是什么。 换句话说,给定 (xu, yu) --> (xd, yd)。然后,您将未失真图像中的每个像素替换为失真图像中对应的像素。

    从 JMBR 开始,我做相反的事情,找到 Rd 作为 Ru 的函数。我明白了:

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

    其中 f 是以像素为单位的焦距(我稍后会解释),r = Ru/f

    我的相机的焦距是 2.5 毫米。我的 CCD 上每个像素的大小是 6 平方毫米。因此 f 为 2500/6 = 417 像素。这可以通过反复试验找到。

    Finding Rd 允许您使用极坐标在失真图像中找到对应的像素。

    每个像素与中心点的角度相同:

    theta = arctan( (yu-yc)/(xu-xc) ) 其中 xc, yc 是中心点。

    那么,

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

    确保你知道你在哪个象限。

    这是我使用的 C# 代码

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

    【讨论】:

      【解决方案4】:

      我盲目地实现了here中的公式,所以我不能保证它会满足你的需要。

      使用auto_zoom 获取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
      

      【讨论】:

        【解决方案5】:

        如果您认为您的公式是精确的,您可以使用 trig 计算一个精确的公式,如下所示:

        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
        

        但是,正如@jmbr 所说,实际的相机失真将取决于镜头和变焦。与其依赖固定公式,不如尝试多项式展开:

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

        通过先调整 A,然后调整高阶系数,您可以计算任何合理的局部函数(展开的形式利用了问题的对称性)。特别是,应该可以计算初始系数来逼近上述理论函数。

        此外,为了获得良好的效果,您需要使用插值过滤器来生成校正后的图像。只要失真不是太大,您就可以使用您将用来线性重新缩放图像的那种过滤器,没有太大问题。

        编辑:根据您的要求,上述公式的等效比例因子:

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

        如果将上述公式与 tan(Rin/f) 一起绘制,您会发现它们的形状非常相似。基本上,在 sin(w) 与 w 有很大不同之前,切线的失真会变得很严重。

        逆公式应该是这样的:

        Rin/f = [Rout/f] / sqrt( sqrt(([Rout/f]^2+1) * (sqrt([Rout/f]^2+1) + 1) / 2 )
        

        【讨论】:

          【解决方案6】:

          description you mention 声明针孔相机(不会引入镜头失真的相机)的投影是由

          R_u = f*tan(theta)
          

          而普通鱼眼镜头相机的投影(即失真)是由

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

          您已经知道 R_d 和 theta,如果您知道相机的焦距(用 f 表示),那么校正图像就相当于根据 R_d 和 theta 计算 R_u。换句话说,

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

          是您正在寻找的公式。估计焦距 f 可以通过校准相机或其他方式来解决,例如让用户提供有关图像校正效果的反馈或使用来自原始场景的知识。

          为了使用 OpenCV 解决相同的问题,您必须获得相机的内在参数和镜头畸变系数。例如,参见Learning OpenCV 的第 11 章(不要忘记查看correction)。然后你可以使用像这个程序这样的程序(用 Python 绑定 OpenCV 编写)来扭转镜头失真:

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

          另请注意,OpenCV 使用的镜头畸变模型与您链接到的网页中的模型完全不同。

          【讨论】:

          • 这取决于是否可以访问相关相机。不知道,我只是来个视频。另外,我觉得相机是批量生产的,单个的变化肯定不会很大吗?原文章中链接的工具,不需要有人拿着棋盘站在镜头前!?
          • 相机参数在同一个相机之间变化,只是通过调整变焦。此外,您可以依靠自动校准技术而不是使用棋盘。无论如何,我已经编辑了我的答案,以解决您问题的第一部分,为您提供您正在寻找的公式。
          • 谢谢jmbr!世界开始变得有意义。但是,我实际上无法让 R_u = f*tan(2*asin(R_d/(2*f))) 给我除了 NaN 之外的任何东西。我对三角学一无所知,这让我现在退缩了:(
          • 这是自动接受的 - 我一直在坚持对我提出的问题做出更明确的回答。
          • @Will:由于赏金系统已经改进,您现在也可以接受另一个答案
          猜你喜欢
          • 1970-01-01
          • 2015-12-24
          • 1970-01-01
          • 2018-04-15
          • 2017-04-27
          • 2018-06-05
          • 2016-05-20
          • 1970-01-01
          相关资源
          最近更新 更多