Программное исправление искажения «рыбий глаз»

В описании, которое вы упомянули, указано, что проекция с помощью камеры с точечным отверстием (без искажения объектива) моделируется

R_u = f*tan(theta)

а проекция обычных камер с объективом «рыбий глаз» (то есть искаженная) моделируется

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

Вы уже знаете R_d и тета, и если бы вы знали фокусное расстояние камеры (представленное f), то исправление изображения равносильно вычислению R_u с точки зрения R_d и тета. Другими словами,

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

формула, которую вы ищете. Оценка фокусного расстояния f может быть решена путем калибровки камеры или других средств, таких как предоставление пользователю возможности оставить отзыв о том, насколько хорошо исправлено изображение, или использование знаний из исходной сцены.

Чтобы решить ту же проблему с помощью OpenCV, вам нужно будет получить внутренние параметры камеры и коэффициенты искажения объектива. См., например, главу 11 Изучение OpenCV (не забудьте проверить исправление). Затем вы можете использовать программу, подобную этой (написанную с привязками 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 использует модель искажения объектива, которая сильно отличается от модели на веб-странице, на которую вы ссылаетесь.

(Исходный постер, представляющий альтернативу)

Следующая функция сопоставляет целевые (прямолинейные) координаты с исходными (искаженными «рыбий глаз») координатами. (Буду признателен за помощь в восстановлении)

Я пришел к этому методом проб и ошибок: я не совсем понимаю, почему этот код работает, благодарны пояснения и повышенная точность!

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 он успешно искажает изображения, используемые в качестве примеров (я не пытался интерполировать качество):

Неработающая ссылка

(А это из поста в блоге, для сравнения :)

Если вы считаете, что ваши формулы точны, вы можете вычислить точную формулу с помощью триггера, например:

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)

Если вы построите приведенную выше формулу рядом с тангенсом (Rin/f), вы увидите, что они очень похожи по форме. По сути, искажение касательной становится серьезным до того, как sin(w) станет сильно отличаться от w.

Обратная формула должна быть примерно такой:

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

Я вслепую реализовал формулы из здесь, поэтому не могу гарантировать, что они сделают то, что вам нужно.

Используйте 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

Я взял то, что сделал 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 мм. Размер каждого пикселя на моей ПЗС-матрице был равен 6 мкм в квадрате. Следовательно, f было 2500/6 = 417 пикселей. Это можно найти методом проб и ошибок.

Поиск Rd позволяет найти соответствующий пиксель в искаженном изображении с использованием полярных координат.

Угол каждого пикселя от центральной точки одинаков:

theta = arctan( (yu-yc)/(xu-xc) ) где xc, yc — центральные точки.

Потом,

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

Убедитесь, что вы знаете, в каком квадранте вы находитесь.

Вот код С#, который я использовал

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

Я нашел этот 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;