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

Question

ОБНОВЛЕНИЕ СТАТУСА КОРПУСА:

Я обнаружил, как сопоставить линейную линзу , от destination координат до source координат.

Как рассчитать радиальное расстояние от центра для перехода от рыбий глаз к прямолинейному?

• 1). Я действительно изо всех сил пытаюсь повернуть его вспять и сопоставить координаты источника с координатами пункта назначения. Что обратного в коде в стиле конвертируемых мной функций?
  

• 2). I also see that my undistortion is imperfect on some lenses - presumably those that are not strictly linear. What is the equivalent to-and-from source-and-destination coordinates for those lenses? Again, more code than just mathematical formulae please...

---

Question as originally stated:

I have some points that describe positions in a picture taken with a fisheye lens.

I want to convert these points to rectilinear coordinates. I want to undistort the image.

I've found это описание о том, как создать эффект «рыбий глаз», но не о том, как его обратить вспять.

Есть также сообщение в блоге , которое описывает, как использовать инструменты для этого; эти картинки из этого:

(1) : SOURCE Оригинальная ссылка на фотографию

Input : Original image with fish-eye distortion to fix.

(2) : DESTINATION Оригинальная ссылка на фото

Output : Corrected image (technically also with perspective correction, but that's a separate step).

How do you calculate the radial distance from the centre to go from fisheye to rectilinear?

My function stub looks like this:

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

Alternatively, I could somehow convert the image from fisheye to rectilinear before finding the points, but I'm completely befuddled by the Документация OpenCV . Есть ли простой способ сделать это в OpenCV, и достаточно ли он эффективен, чтобы делать это в прямом эфире?

Answer 1

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

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, вам необходимо получить внутренние параметры камеры и коэффициенты искажения объектива. См., Например, главу 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 использует модель искажения объектива, совершенно отличную от модели на веб-странице, на которую вы ссылаетесь.

Answer 2

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

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

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

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

Мертвая ссылка

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

Answer 3

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

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), вы увидите, что они очень похожи по форме. По сути, искажение касательной становится серьезным до того, как грех (w) сильно отличается от w.

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

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

Answer 4

Я взял то, что сделал 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

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

Вот код 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;
      }
}

Answer 5

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

Answer 6

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

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