Перейти к содержимому
Compvision.ru
ProgRoman

Устранение радиальной дисторсии

Recommended Posts

Добрый день, стоит такая задача есть камера типа "рыбий глаз" с него идёт видео и необходимо видео преобразовывать так что бы избавляться от этого эффекта рыбьего глаза или(радиальной дисторсии как я понял) тема для меня достаточно новая, как я понял так как видео идёт с одной камеры и она неподвижна, то откалибровать нормально невозможно.. поэтому пользоваться функцией 

undistort(frame,dst,cameraMatrix,distCoeffs);

так же невозможно, можно ли как-то корректно преобразовать не зная ни cameraMatrix ни distCoeffs...

---

В архиве ролик с камерыsmile.png

3.rar

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

1. Автоматически откалибровать возможно, почему нет? И undistort применить потом тоже.

2. Сохрани из видео один кадр и попробуй в том же Гимпе исправить искажения так, как тебе видится верным. Там есть фильтр типа "Исправить искажение линз". Получилось? Тогда открой исходники Гимпа и посмотри на функцию, которая это делает. Там всего три параметра (или четыре?), которые устанавливаются ползунками. Функция простая, можно переписать на OpenCV за полчаса.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Что получилось,smile.png  Пробовал устранять дисторсию вручную с помощью gimp, но как то особо не получилось возможно из за большого ракурса съёмки не знаю или может просто я не смог правильно параметры подобрать)  Посмотрел в интернете оказалось что снимки сделанные камерой рыбий глаз переводят в панорамные снимки

Mat convert(const Mat& src)
{
	int midy = src.size().height / 2;
	int midx = src.size().width / 2;
	int maxmag = (midy > midx ? midy : midx);
	int circum = 2 * CV_PI * maxmag;     // c = 2*pi*r

	Mat dst(Size(circum,maxmag),src.type());

	Vec3b black(0,0,0);
	for (int y = 0; y < maxmag; y++) 
	{
		for (int x = 0; x < circum; x++) 
		{
			double theta = -1.0 * x / maxmag;       
			double mag = maxmag - y;                
			
			int targety = int(midy + mag * cos(theta) + 0.5);
			int targetx = int(midx + mag * sin(theta) + 0.5);
			int ix =-(x * 2.0 * CV_PI / src.size().width);
			int iy = y * 1.0 * maxmag / src.size().height;

			if (targety < 0 || targety >= src.size().height || targetx < 0 || targetx >= src.size().width) 
			{
				dst.at<Vec3b>(Point(x,y)) = black;
			} 
			else 
			{
				dst.at<Vec3b>(Point(x,y)) = src.at<Vec3b>(targety,targetx );				
			}
		}
	}

	return dst;
}

Панораму я конечно получаю, но всё равно она сильно искажена, думаю её надо как-то делить на несколько частей и повторно пытаться убрать искажения прилагаю картинки(res1, res2,res3).

так же нашёл код который переводит изображение из fisheye в прямоугольник прилагаю класс

class Analyzer
 {
 private:
      vector<double> mFisheyeCorrect;
      int mFELimit;
      double mScaleFESize;

 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
		  mFELimit = 1500;
		  mScaleFESize = 0.9;
            //i corresponds to Rundist/focalLengthInPixels * 1000 (to get integers)
          for (int i = 0; i < mFELimit; i++)
          {
              double result = sqrt(1 - 1 / sqrt(1.0 + (double)i * i / 1000000.0)) * 1.4142136;
			  mFisheyeCorrect.push_back(result);
          }
      }

      Mat RemoveFisheye(Mat& aImage, double aFocalLinPixels)
      {
		  Mat correctedImage(aImage.rows, aImage.cols,aImage.type());
		  //The center points of the image
          double xc = aImage.size().width / 2.0;
          double yc = aImage.size().height / 2.0;

          bool xpos, ypos;
		  //Move through the pixels in the corrected image; 
		  //set to corresponding pixels in distorted image
          for (int i = 0; i < correctedImage.size().width; i++)
          {
			  for (int j = 0; j < correctedImage.size().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 = atan2(ydif, xdif);
                     //find index for lookup table
                  int index = (int)(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 = abs(Rd*cos(theta));
                  double ydelta = abs(Rd * sin(theta));
                     //convert to pixel coordinates
                  int xd = (int)(xc + (xpos ? xdelta : -xdelta));
                  int yd = (int)(yc + (ypos ? ydelta : -ydelta));
				  
                  xd = std::max(0, min(xd, aImage.size().width-1));
                  yd = std::max(0, min(yd, aImage.size().height-1));
                     //set the corrected pixel value from the distorted image
                  correctedImage.at<Vec3b>(Point(i, j)) = aImage.at<Vec3b>(Point(xd, yd));
              }
          }
          return correctedImage;
      }
};

/*
так же функция для вычисления фокусного расстояния
*/

double dist(double x, double y)
{
	return sqrt(x*x + y*y);
}

double calculate_focal_distance(int src_size_x, int src_size_y, double crop_factor)
{
	double f = dist(src_size_x, src_size_y)*crop_factor/CV_PI;

	return f;
}

Результаты работы этого класса  это изображения (rect1, rect2, rect3)

Сейчас пытаюсь делать преобразования из сферических координат, но что-то явно делаю не так буду разбираться код ниже

Mat convert6(const Mat& src)
{
	Point pfish;
	float theta, phi,r;
	Point3f psph;
	Mat dst(src.rows,src.cols,src.type());

	float FOV = CV_PI;// FOV of the fisheye, eg: 180 degrees
	float width = src.size().width;
	float heigh = src.size().height;
	Vec3b balck(0,0,0);
	for(int y = 0; y<heigh; y++)
	{
		for(int x = 0; x<width; x++)
		{
			//polar angles			
			theta = 2.0 * CV_PI * (x / width - 0.5); // -pi to pi
			phi = CV_PI * (y / heigh - 0.5);	// -pi/2 to pi/2

			//vector in 3D space
			psph.x = cos(phi)*sin(theta);
			psph.y = cos(phi)*cos(theta);
			psph.z = sin(phi);

			// Calculate fisheye angle and radius			
			theta = atan2(psph.z,psph.x);
			phi = atan2(sqrt(psph.x*psph.x+psph.z*psph.z),psph.y);
			r = width * phi / FOV; 

			// Pixel in fisheye space
			pfish.x = 0.5 * width + r * cos(theta);
			pfish.y = 0.5 * width + r * sin(theta);

			if(((pfish.x>0)&&(pfish.x<width))&&((pfish.y>0)&&(pfish.y<heigh)))
			{
				dst.at<Vec3b>(pfish) = src.at<Vec3b>(Point(x,y));
			}			
		}
	}

	return dst;
}

Результаты этого кода это изображения(sphr1,sphr2,sphr3)

Есть ещё код в матлабе, но не знаю может в силу того что я его плохо знаю как у меня так и не получилось адаптировать его для изображений привожу код

clear;

M=fspecial('gaussian',256,32); % generate fake image

 X0=size(M,1)/2; Y0=size(M,2)/2;
 [Y X z]=find(M);
 X=X-X0; Y=Y-Y0;
 theta = atan2(Y,X);
 rho = sqrt(X.^2+Y.^2);

 % Determine the minimum and the maximum x and y values:
 rmin = min(rho); tmin = min(theta);
 rmax = max(rho); tmax = max(theta);

 % Define the resolution of the grid:
 rres=128; % # of grid points for R coordinate. (change to needed binning)
 tres=128; % # of grid points for theta coordinate (change to needed binning)      
 
 F = TriScatteredInterp(rho,theta,z,'natural');

 %Evaluate the interpolant at the locations (rhoi, thetai).
 %The corresponding value at these locations is Zinterp:

 [rhoi,thetai] = meshgrid(linspace(rmin,rmax,rres),linspace(tmin,tmax,tres));
 Zinterp = F(rhoi,thetai);

 subplot(1,2,1); imagesc(M) ; axis square
 subplot(1,2,2); imagesc(Zinterp) ; axis square

Добавлю ещё как мне кажется полезные ссылки 

http://paulbourke.net/dome/fish2/

http://paulbourke.net/dome/fishtilt/

http://paulbourke.net/dome/2fish/

http://www.kscottz.com/fish-eye-lens-dewarping-and-panorama-stiching/

http://stackoverflow.com/questions/12924598/examples-to-convert-image-to-polar-coordinates-do-it-explicitly-want-a-slick-m

http://stackoverflow.com/questions/2477774/correcting-fisheye-distortion-programmatically

http://mathworld.wolfram.com/CylindricalEquidistantProjection.html

http://mathworld.wolfram.com/OrthographicProjection.html

 

post-2515-0-89436100-1418027981_thumb.jp

post-2515-0-35481100-1418027999_thumb.jp

post-2515-0-28481700-1418028017_thumb.jp

post-2515-0-14639200-1418028031_thumb.jp

post-2515-0-50882000-1418028045_thumb.jp

post-2515-0-19522200-1418028059_thumb.jp

post-2515-0-24721900-1418028090_thumb.jp

post-2515-0-80573700-1418028109_thumb.jp

post-2515-0-90494300-1418028147_thumb.jp

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

 

в  opencv можно задать любое преобразование через remap чтобы не заниматься самому интерполяцией.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Ок спасибо большое!

по первой ссылке есть код SelfCalibration.cpp от Smorodov который должен устранять дисторсии как я понял, скачал попробовал на тестовом изображении с решётками и всё получилось, попробовал на других изображениях.. как-то особо не получилось) привожу картинку  

post-2515-0-20498200-1418110215_thumb.jp

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Дык, код работает с решетками (или решеткообразными изображениями).

Фоткаете решетку, и калибруете камеру, правда не уверен будет ли она работать для очень сильных искажений.

 

Затем используете вычисленные параметры для исправления изображений.

Модель искажений можете задать свою. 

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

вам сначала надо определиться с типом дисторсии, по идее вам надо сделать подушку(в противоввес бочке)

 

еще формулы есть тут

http://docs.opencv.org/doc/tutorials/calib3d/camera_calibration/camera_calibration.html

можно поиграть с коэффициентами undistortPoints

http://docs.opencv.org/modules/imgproc/doc/geometric_transformations.html

 

гимп кстати не помогает, выкрутил main настройку на максимум, не знаю правда какая у них там формула.

http://docs.gimp.org/en/plug-in-lens-distortion.html

fisheye.jpg

fisheye_undistorted_gimp.jpg

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Вот код набросал чтобы поиграть с встроенной функцией коррекции дисторсии OpenCV:
 

#include <iostream>
#include <opencv2/opencv.hpp>


int main(int argc, char* argv[])
{
	cv::namedWindow("result");
	cv::Mat img = cv::imread("E:\\ImagesForTest\\fisheye.jpg");
	imshow("src",img);

	// 147,226,67,56,0,0,147,131


	int ifx =  128;
	int ify = 	128;

	int icx = 128;
	int icy = 128;

	int ik1 = 128;
	int ik2 = 128;

	int ip1 = 128;
	int ip2 = 128;

	cv::namedWindow("Control");
	// Create trackbars in "Control" window
	cv::createTrackbar("fx", "Control", &ifx, 256);
	cv::createTrackbar("fy","Control", &ify, 256);
	cv::createTrackbar("cx", "Control", &icx, 256);
	cv::createTrackbar("cy","Control", &icy, 256);
	cv::createTrackbar("k1", "Control", &ik1, 256);
	cv::createTrackbar("k2","Control", &ik2, 256);
	cv::createTrackbar("p1", "Control", &ip1, 256);
	cv::createTrackbar("p2","Control", &ip2, 256);
	int k=0;
	float fx, fy, cx, cy, k1, k2, p1, p2;

	cv::Mat intrinsics=cv::Mat::zeros(3,3,CV_64FC1);
	cv::Mat dist_coeffs=cv::Mat::zeros(1,4,CV_64FC1);
	cv::Mat mapx(img.size(),CV_32FC1);
	cv::Mat mapy(img.size(),CV_32FC1);
	cv::Mat dst(img.size(),img.type());

	while(k!=27)
	{
		fx=(ifx-128)/128.0;
		fy=(ify-128)/128.0;
		cx=(icx-128)/128.0;
		cy=(icy-128)/128.0;
		k1=(ik1-128)/128.0;
		k2=(ik2-128)/128.0;
		p1=(ip1-128)/128.0;
		p2=(ip2-128)/128.0;

		fx*=1000;
		fx+=500;

		fy*=1000;
		fy+=500;

		cx*=1000;
		cx+=img.cols;

		cy*=1000;
		cy+=img.rows;

		k1*=5;
		k2*=5;


		intrinsics.at<double>(0,0)=fx;
		intrinsics.at<double>(1,1)=fy;
		intrinsics.at<double>(2,2)=1.0;
		intrinsics.at<double>(0,2)=cx;
		intrinsics.at<double>(1,2)=cy;

		dist_coeffs.at<double>(0,0)=k1;
		dist_coeffs.at<double>(0,1)=k2;
		dist_coeffs.at<double>(0,2)=p1;
		dist_coeffs.at<double>(0,3)=p2;

		cv::initUndistortRectifyMap(intrinsics,dist_coeffs,cv::Mat(),intrinsics,img.size(),CV_32FC1,mapx,mapy);

		cv::remap(img,dst,mapx,mapy,0);

		imshow("dst",dst);
		k=cv::waitKey(10);
	}
	return 0;
}

Вот что у меня получилось после нескольких минут подстройки параметров:

post-1-0-84534100-1418132483_thumb.png

 

Результат не самый хороший, но углы почти прямоугольные.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Попробовал еще 1 метод из GIMP, стены выпрямились, а стол нет.

gimp_curve.PNG

 

Что если создать такую тулзу: отмечаем на изображениее кривые точками которые после преобразования должны перейти в прямые, вопрос только как потом получить "поле" пикселей для ремапа?

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

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

 

Попробовал еще 1 метод из GIMP, стены выпрямились, а стол нет.

А стол, скорее всего, такой и есть, с криволинейной кромкой.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Интересно, а математически как формулируется задача исправления такой дисторсии?

По найденным координатам углов доски строится матрица параметров камеры. Далее координаты каждого пикселя умножаются на неё переходя в исправленные координаты. Промежутки аппроксимируются. Так?

Может, лучше сделать что-то типа того, что предлагает mrgloom? Вывести на изображение виртуальную сетку с большим числом узлов, расставить узлы, искривив при этом сетку. После сделать тот же опенсивишный remap, только map будет заполнен:

1. вершинами сетки;

2. промежуточные значения для map интерполированы (сплайнами теми же).

Как такая идея полуручного исправления?

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Теоретически делается так:

Дано:

1) координаты точек сетки на калибровочной доске (x,y)

2) координаты тех же точек на изображении (x',y')

3) модель искажений в виде уравнений вида:

 

    a1*f1(x)+a2*f2(x)+...+an*fn(x)+ak*fk(y)+...+am*fm(y)=x'

    b1*f1(x)+b2*f2(x)+...+bn*fn(x)+bk*fk(y)+...+bm*fm(y)=y'

 

f1,f2,fn,fk,fm - некоторые функции, например x^2, x^3 и им подобные, но не обязательно, зависит от модели искажений.

 

Ищем:

Коэффициенты a1...am и b1...bm

 

Метод:

 

Решение переопределенной СЛАУ относительно коэффициентов.

 

    a1*f1(x1)+a2*f2(x1)+...+an*fn(x1)+ak*fk(y1)+...+am*fm(y1)=x1'

    b1*f1(x1)+b2*f2(x1)+...+bn*fn(x1)+bk*fk(y1)+...+bm*fm(y1)=y1' 

................

    a1*f1(xL)+a2*f2(xL)+...+an*fn(xL)+ak*fk(yL)+...+am*fn(yL)=xL'

    b1*f1(xL)+b2*f2(xL)+...+bn*fn(xL)+bk*fk(yL)+...+bm*fn(yL)=yL'

 

x1,xL,y1,yL координаты точек.

 

Система, в общем случае переопределенная, поэтому решается через SVD, что дает приближение по методу наименьших квадратов.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

4d7339dc0a90d158c7d43613c5321569a406193b

37d5f1e28b4dbaec6f828add97e33890c80b3344

462b15f3a7f036f979a52f32329304596ce7fc2f

 

по ссылке выше там 5 коэффициентов

For the distortion OpenCV takes into account the radial and tangential factors.

 

Мы тут вроде как пытаемся исправить радиальную дисторсию.

 

в вики более общая формула

http://en.wikipedia.org/wiki/Distortion_%28optics%29

161fb9bdab56a77849a8025f110ef562.png

81f0f81d0f2a2f830f71d5d4580d4613.png

 

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

 

 

а вообще вот тема была

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Подправил код со сферическими координатами теперь он работает гораздо лучше кидаю исправленный код и полученные картинки

Mat convert(const Mat& src)
{
	Point pfish;
	float theta, phi,r;
	Point3f psph;
	Mat dst(2*src.rows,src.cols,src.type());

	float FOV = CV_PI;// FOV of the fisheye, eg: 180 degrees
	float width = dst.size().width;
	float heigh = dst.size().height;
	Vec3b balck(0,0,0);
	for(int y = 0; y<heigh; y++)
	{
		for(int x = 0; x<width; x++)
		{
			//polar angles			
			theta = 2.0 * CV_PI * (x / width - 0.5); // -pi to pi
			phi = CV_PI * (y / heigh - 0.5);	// -pi/2 to pi/2			

			//vector in 3D space
			psph.x = cos(phi)*sin(theta);
			psph.y = cos(phi)*cos(theta);
			psph.z = sin(phi);

			// Calculate fisheye angle and radius			
			theta = atan2(psph.z,psph.x);
			phi = atan2(sqrt(psph.x*psph.x+psph.z*psph.z),psph.y);
			r = width * phi / FOV; 

			// Pixel in fisheye space
			pfish.x = 0.5 * width + r * cos(theta);
			pfish.y = 0.5 * width + r * sin(theta);

			if(((pfish.y>0)&&(pfish.y<int(heigh/2)))&&((pfish.x>0)&&(pfish.x<(int)(width))))
			{
				//dst.at<Vec3b>(pfish) = src.at<Vec3b>(Point(x,y));
				dst.at<Vec3b>(Point(x,y)) = src.at<Vec3b>(pfish);
			}
		}
	}

	return dst;
}

post-2515-0-91870800-1418370108_thumb.jp

post-2515-0-75163500-1418370126_thumb.jp

  • Like 1

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Попробовал метод из Image Deformation Using Moving Least Squares

http://www.mathworks.com/matlabcentral/fileexchange/12249-moving-least-squares

 

Проблема в том, чтобы получить приемлимый трансформ придётся задавать точки с опредеденной густотой по всему изображению.

мои эксперименты

https://dl.dropboxusercontent.com/u/8841028/fisheye/MLS_0.PNG

https://dl.dropboxusercontent.com/u/8841028/fisheye/MLS_1.PNG

https://dl.dropboxusercontent.com/u/8841028/fisheye/MLS_2.PNG

https://dl.dropboxusercontent.com/u/8841028/fisheye/MLS_3.PNG

  • Like 1

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Создайте учётную запись или войдите для комментирования

Вы должны быть пользователем, чтобы оставить комментарий

Создать учётную запись

Зарегистрируйтесь для создания учётной записи. Это просто!

Зарегистрировать учётную запись

Войти

Уже зарегистрированы? Войдите здесь.

Войти сейчас


  • Сейчас на странице   0 пользователей

    Нет пользователей, просматривающих эту страницу

×