程序化纠正鱼眼畸变

Question

程序化纠正鱼眼畸变

mathgraphicsgeometryprojection

56

赏金状态更新：

我发现了如何将线性镜头从目标坐标映射到源坐标。，即从destination坐标到source坐标的映射。

如何计算从鱼眼到矩形的径向距离？

1）我实际上很难反过来，将源坐标映射到目标坐标。请问在我发布的转换函数风格中，逆变换是什么？
2）我还发现我的去畸变对某些镜头不完美-可能是那些不严格线性的镜头。对于这些镜头，相应的源坐标和目标坐标是什么？请再次提供更多代码而不仅仅是数学公式...

我有一些描述鱼眼镜头拍摄图像中位置的点。

我想将这些点转换为直角坐标系。我想去除图像失真。

我找到了如何产生鱼眼效果的此描述，但没有找到如何反转它。

还有一个博客文章描述了如何使用工具来完成它; 这些图片来自那里:

(1) : SOURCE 原始照片链接

Input : 带鱼眼失真的原始图像。

(2) : DESTINATION 原始照片链接

Output : 纠正后的图像（技术上还包括透视校正，但这是一个单独的步骤）。

如何计算从中心到边缘的径向距离以从鱼眼转换为直角？

我的函数存根看起来像这样：

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文档让我完全困惑了。在OpenCV中有一种简单的方法吗？它的性能是否足够好，可以用于实时视频流？

- Will

我不太明白您在寻找什么。鱼眼图是从球面到图片平面的映射，反之则是从图片平面回到球面，对吗？您要寻找哪个直角坐标系？ - mtrw

1

@mtrw 我的源图像是鱼眼畸变的，我想要对其进行去畸变处理。 - Will

1

Will，你最终得到了一个确定的答案吗？我非常想看看你最终写出的任何代码。 - Nestor

@Will，我知道这个问题已经问了很久，但我有一个类似的问题 - 你是否得到了包括透视校正在内的功能代码片段转换？ - jogall

@jogal 我最终使用了我在下面发布的解决方案；祝你好运！ - Will

显示剩余3条评论

7个回答

9

（原帖提供替代方案）

以下函数将目标（直角坐标系）坐标映射到源（鱼眼畸变）坐标。 （我希望能够反向操作，需要帮助）

我通过试错得出了这个结果：我并不完全理解这段代码为什么有效，欢迎解释和改进以提高准确性！

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的因子时，它可以成功地去畸变所示的图像（我没有尝试进行质量插值）： ~~链接已失效~~ （以下是博客文章中的内容，供参考：）使用Panotools

- Will

1

你的代码之所以能够工作，是因为你将（rx，ry）按比例因子theta进行了缩放（现在它是一个比率而不是角度）。如果原始镜头具有从视角到图像偏移的“线性映射”（正如维基百科文章所述），我相信atan（r）/ r是正确的。 - comingstorm

2

你的映射反转应该按比例缩放tan(r')/r'的因子，其中r'是未失真图像中心的半径。 - comingstorm

2

而这两种方法都行的原因是，如果你有向量v' = v*k(|v|)，且你想让|v'|=f(|v|)，你可以对第一个式子取绝对值：|v'|=|v|*k(|v|)=f(|v|) -- 所以k(|v|)=f(|v|)/|v|。 - comingstorm

@comingstorm，非线性映射的等效方法是什么？ - Will

@stkent，所以我正在使用R进行这个操作，它运行得很好，除了根据我选择的因子，我会看到在拉伸或缩放图像时叠加在上面的'波浪'。也就是说，我得到了所需的桶形畸变，但是之后还有一些类似于灰色波浪的伪像。想知道这是否是某种边界问题或相移类型的问题？示例：https://imgur.com/a/lBNXPYh - SqueakyBeak

4

如果您认为您的公式是精确的，您可以使用三角函数计算出一个精确的公式，如下所示：

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 )

- comingstorm

3

我盲目地实现了这里的公式，所以我不能保证它能满足您的需求。

使用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

- Roman Zenka

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个像素。这可以通过试错方法找到。

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

- Barry Vant-Hull

2

我发现了这个pdf文件，我已经证明了数学是正确的（除了一行代码vd = *xd**fv+v0应该改为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;

- Myke Smyth

0

这可以作为一个优化问题来解决。只需在图像中绘制应该是直线的曲线。存储每个曲线的轮廓点。现在我们可以将鱼眼矩阵作为最小化问题来解决。最小化点中的曲线，这将给我们一个鱼眼矩阵。它有效。

也可以通过调整鱼眼矩阵使用滑动条手动完成！这里有一个使用OpenCV进行手动校准的鱼眼GUI代码。

- Jerin Antony

404页面未找到。 - Andrea Nicolai

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jmbr · Accepted Answer

你提到的描述指出，通过针孔相机（不引入镜头畸变的相机）的投影可以建模为

R_u = f*tan(theta)

常见鱼眼镜头拍摄的投影（即，畸变）由以下模型描述：

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

您已经知道R_d和theta，如果您知道相机的焦距（用f表示），那么校正图像就相当于计算R_u与R_d和theta之间的关系。换句话说，

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

这是你要找的公式。可以通过校准相机或其他手段（例如让用户反馈图像纠正的情况或使用原始场景的知识）来估算焦距f。

为了使用OpenCV解决同样的问题，你需要获取相机的内部参数和镜头畸变系数。例如，参见学习OpenCV的第11章（不要忘记检查更正）。然后，你可以使用类似于这个程序（使用OpenCV的Python绑定编写）的程序来反转镜头畸变：

#!/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使用的镜头畸变模型与您链接的网页中使用的模型非常不同。