Как исправить дисторсию на видео

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

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

Основная идея в том, что Virtual Dub позволяет работать с последовательностью файлов изображений. Это касается и сохранения видео, и загрузки. Поэтому общий сценарий очевидный: разложить видео на последовательность картинок, у картинок исправить дисторсию, а затем собрать из них видео.

Сперва про открытие видео Vision+ в Virtual Dub. Есть разные способы, но я пользуюсь обычным Virtual Dub без расширений. Видео с Vision+ при съемке записывается в контейнере mp4. Для открытия в Virtual Dub перегоняю видео в контейнер avi при помощи программы MP4Cam2AVI. И еще, если правильно помню, нужно поставить кодек x264vfw для понимания H264 Virtual Dub’ом.

Когда видео-файл открыт в Virtual Dub, то создание последовательности файлов изображений выполняется через File -> Export -> Image Sequence…

И теперь задача резко упрощается для меня, поскольку исправлять дисторсию изображений можно несколькими способами по вкусу. Этот вопрос не рассматриваю. Лично я пользуюсь программами для работы с панорамами (только на выходе не панорама, а последовательность изображений, представляющих каждый слой отдельно). Важно в выходных картинках подобрать правильное разрешение и соотношение сторон.

После генерации последовательность исправленных изображений загружается в Virtual Dub, после этого следует обычный процесс обработки.

Да, этот метод далек от идеала, поскольку требует основательных накладных расходов (место и время для сохранения изображений, их обработки, загрузки). Альтернативным способом является использование внутри Virtual Dub соответствующих фильтров. Например, Barrel Distortion Correction или De-fishr Plugin. В более продвинутом ПО редактирования видео есть свои способы правки дисторсии.

О результате. Сравнение кадра до и после исправления дисторсии:

И тестовое видео (снято на Иркутском водохранилище):

До исправления:

После исправления:

Время на прочтение
5 мин

Количество просмотров 24K

В данной статье будет рассказываться о применении библиотеки машинного зрения (openCV) для удаления эффекта радиального искажения (дисторсии) с фото и видео. Данный эффект также известен как эффект рыбьего глаза (fisheye) или distortion. Решение написать данную статью было принято после нескольких дней поиска информации в интернете. Не смотря на то, что есть гайды на английском языке, они не объясняют как правильно установить openCV, чтобы все работало. В статье присутствует готовый код.

Сразу привожу фото итогового результата. Слева оригинальное фото, справа — обработанное:

Сборка и установка openCV

Первое, что нужно сделать, это грамотно установить библиотеку openCV. Для этого скачиваем из официального репозитория два проекта — openCV и opencv_contrib.

git clone https://github.com/opencv/opencv.git
git clone https://github.com/opencv/opencv_contrib.git

Пока загружается openCV, устанавливаем видеокодек ffmpeg:

sudo apt-get install ffmpeg

Заходим в папку openCV, создаем подпапку buid и заходим в нее. Вся работа по сборке и установке библиотеки openCV будет производиться из этой директории.

cd opencv
mkdir build
cd build/

Для сборки библиотеки выполняем следующие команды:

cmake .. -DOPENCV_EXTRA_MODULES_PATH=/путь к папке opencv_contrib/modules/ /путь к папке opencv/
make -j5
sudo make install

У меня сборка заняла около полутора часов, установка — несколько минут. Обратите внимание: если у вас возникла ошибка при сборке(выполнение команды cmake), для нового запуска необходимо удалить файл CMakeCache.txt. После установки можем проверить все ли правильно получилось. Для этого можно вызвать рабочую среду python и импортировать библиотеку openCV. Если никаких ошибок не произошло, то вы все сделали правильно. Вторая строчка покажет, какая версия у вас установлена. На момент написания статьи я пользовался 3 версией библиотеки.

import cv2
print ("OpenCV version : {0}".format(cv2.__version__))

Калибровка камеры

Для того, чтобы убрать искажения, нам необходимо определить калибровочные коэффициенты для нашей камеры. Для этого необходимо скачать картинку с шахматной доской, сделать 5-6 снимков на камеру, изображения с которой мы хотим обработать. Все изображения необходимо конвертировать в формат png. Далее выполняем следующий код:

from __future__ import print_function
import numpy as np
import cv2
from common import splitfn
import os

if __name__ == '__main__':
    import sys
    import getopt
    from glob import glob

    args, img_mask = getopt.getopt(sys.argv[1:], '', ['debug=', 'square_size='])
    args = dict(args)
    args.setdefault('--debug', '/рабочая директория/')
    args.setdefault('--square_size', 1.0)
    if not img_mask:
        img_mask = '/папка с изображениями/*.png'
    else:
        img_mask = img_mask[0]

    img_names = glob(img_mask)
    debug_dir = args.get('--debug')
    if not os.path.isdir(debug_dir):
        os.mkdir(debug_dir)
    square_size = float(args.get('--square_size'))

    pattern_size = (9, 6)
    pattern_points = np.zeros((np.prod(pattern_size), 3), np.float32)
    pattern_points[:, :2] = np.indices(pattern_size).T.reshape(-1, 2)
    pattern_points *= square_size

    obj_points = []
    img_points = []
    h, w = 0, 0
    img_names_undistort = []
    for fn in img_names:
        print('processing %s... ' % fn, end='')
        img = cv2.imread(fn, 0)
        if img is None:
            print("Failed to load", fn)
            continue

        h, w = img.shape[:2]
        found, corners = cv2.findChessboardCorners(img, pattern_size)
        if found:
            term = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_COUNT, 30, 0.1)
            cv2.cornerSubPix(img, corners, (5, 5), (-1, -1), term)

        if not found:
            print('chessboard not found')
            continue

        img_points.append(corners.reshape(-1, 2))
        obj_points.append(pattern_points)

        print('ok')

    rms, camera_matrix, dist_coefs, rvecs, tvecs = cv2.calibrateCamera(obj_points, img_points, (w, h), None, None)

    print("nRMS:", rms)
    print("camera matrix:n", camera_matrix)
    print("distortion coefficients: ", dist_coefs.ravel())

    cv2.destroyAllWindows()

В результате выполнения данного скрипта в консоли появится сообщение об обработанных фото и отобразятся два важных параметра — camera matrix и distortion coefficients. Это и есть те калибровочные коэффициенты, которые нам нужны.

Обработка фото и видео

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

from __future__ import print_function
import numpy as np
import cv2
import glob
from matplotlib import pyplot as plt
from common import splitfn
import os

img_names_undistort = [img for img in glob.glob("/путь до папки с фотографиями/*.png")]
new_path = "/путь для сохранения обработанных изображений/"

camera_matrix = np.array([[1.26125746e+03, 0.00000000e+00, 9.40592038e+02],
                          [0.00000000e+00, 1.21705719e+03, 5.96848905e+02],
                          [0.00000000e+00, 0.00000000e+00, 1.00000000e+00]]);
dist_coefs = np.array([-0.49181345,  0.25848255, -0.01067125, -0.00127517, -0.01900726]);

i = 0

#for img_found in img_names_undistort:
while i < len(img_names_undistort):
        img = cv2.imread(img_names_undistort[i])
    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

        h,  w = img.shape[:2]
        newcameramtx, roi = cv2.getOptimalNewCameraMatrix(camera_matrix, dist_coefs, (w, h), 1, (w, h))

        dst = cv2.undistort(img, camera_matrix, dist_coefs, None, newcameramtx)

    dst = cv2.cvtColor(dst, cv2.COLOR_BGR2RGB)

        # crop and save the image
        x, y, w, h = roi
        dst = dst[y:y+h-50, x+70:x+w-20]

    name = img_names_undistort[i].split("/")
    name = name[6].split(".")
    name = name[0]
    full_name = new_path + name + '.jpg'

        #outfile = img_names_undistort + '_undistorte.png'
        print('Undistorted image written to: %s' % full_name)
        cv2.imwrite(full_name, dst)
    i = i + 1

from __future__ import print_function
import numpy as np
import cv2
import glob
from matplotlib import pyplot as plt
from common import splitfn
import os

FILENAME_IN = "videoin.mp4"
FILENAME_OUT = "videoout.mp4"
CODEC = 'mp4v' 

camera_matrix = np.array([[1.26125746e+03, 0.00000000e+00, 9.40592038e+02],
                          [0.00000000e+00, 1.21705719e+03, 5.96848905e+02],
                          [0.00000000e+00, 0.00000000e+00, 1.00000000e+00]]);
dist_coefs = np.array([-3.18345478e+01, 7.26874187e+02, -1.20480816e-01, 9.43789095e-02, 5.28916586e-01]);

print ("OpenCV version : {0}".format(cv2.__version__))
print((cv2.__version__).split('.'))
# Load video
video = cv2.VideoCapture(FILENAME_IN)

fourcc = cv2.VideoWriter_fourcc(*list(CODEC))

fps = video.get(cv2.CAP_PROP_FPS)

frame_count = video.get(cv2.CAP_PROP_FRAME_COUNT)

size = (int(video.get(cv2.CAP_PROP_FRAME_WIDTH)),
        int(video.get(cv2.CAP_PROP_FRAME_HEIGHT)))
sizew = (1676, 846)
writer = cv2.VideoWriter(FILENAME_OUT, fourcc, 25, sizew)

newcameramtx, roi = cv2.getOptimalNewCameraMatrix(camera_matrix, dist_coefs, (size[0], size[1]), 1, (size[0], size[1]))
x, y, w, h = roi
M = cv2.getRotationMatrix2D((size[0]/2,size[1]/2),5,1)

while video.grab() is True:
    print("On frame %i of %i."%(video.get(cv2.CAP_PROP_POS_FRAMES), frame_count))

    frame = video.retrieve()[1]
    frame = cv2.undistort(frame, camera_matrix, dist_coefs, None, newcameramtx)
    frame = frame[y:y+h-50, x+70:x+w-20]

    writer.write(frame)
video.release()
writer.release()

() translation by (you can also view the original English article)

Искажение широкоугольного объектива

Оптические искажения объектива на видео с дрона — это обычная вещь, поскольку, как правильно, система камеры на дроне имеет широкое поле обозрения (камера GoPro — классический пример такого эффекта). Это здорово иметь широкое поле обозрения для демонстрации ландшафта. И все же, если вы немного опустите камеру ниже или выше линии горизонта, вы немедленно заметите появление приличного оптического искажения, которое часто придает футажам вид «рыбьего глаза» по причине оптических изгибов. К счастью, это можно легко исправить в Adobe After Effects.

Optics Compensation

Использование эффекта Optics Compensation в After Effects быстро избавит ваше видео от любого широкоугольного искажения. Прежде всего, убедитесь, что параметр Reverse Lens Distortion у вас активный. Далее укажите значение Field of View: 75.

Опция Optimal Pixels

Если вы хотите пойти дальше в использовании этого эффекта и получить более широкую картинку без искажения, вы можете активировать опцию Optimal Pixels. Для достижения этого After Effects все же согнет футаж так, что внизу и вверху мы получим черные области. Самый простой способ это осуществить — установить некоторую пропорцию, формат кадра (например, 2,35). Я добавил изображение с пропорцией 2,35, которое вы можете использовать в вашем проекте над своим футажом.

Бонусная подсказка

Если вы хотите, чтобы ваше видео с дрона ещё более выделялось среди других последующим монтажом, обязательно посмотрите набор The Handy Seamless Transitions Pack от Videohive. Он включает более 490 переходов для After Effects, имеет хороший интерфейс, понятную навигацию по всем опциям. Подобные стили переходов приобрели невероятную популярность за последний год и, опираясь на видео превью, можно с легкостью сказать, что они выглядят даже лучше, чем на футаже воздушной съемки!

В этой статье упоминались:

• Широкоугольные объективы
• GoPro
• Handy Seamless Transitions
• Видео воздушной съемки с дрона от Ozark Drones

Что получилось,  Пробовал устранять дисторсию вручную с помощью 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/

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