Описание нелинейного преобразования между двумя изображениями с использованием гомографии

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

Матричное преобразование является линейным, а линейные преобразования всегда отображают прямые линии в прямые (http://en.wikipedia.org/wiki/Linear_map). Из рисунка видно, что преобразование является нелинейным, поэтому его нельзя описать с помощью матричной операции.

Тем не менее, вы можете использовать модель искажения объектива, подобную той, что используется в OpenCV (http://docs.opencv.org/doc/tutorials/calib3d/camera_dication/camera_dication.html), и получение коэффициентов не должно быть очень сложным. Вот что вы можете сделать в Matlab:

Назовите (x, y) координаты исходной точки (верхнее изображение) и (xp, yp) координаты искаженной точки (нижнее изображение), обе сдвинутые к центру изображения и разделенные на коэффициент масштабирования (то же самое для x и y), поэтому они лежат более или менее в интервале [-1, 1]. Модель искажения:

x = ( xp*(1 + k1*r^2 + k2*r^4 + k3*r^6) + 2*p1*xp*yp + p2*(r^2 + 2*xp^2));
y = ( yp*(1 + k1*r^2 + k2*r^4 + k3*r^6) + 2*p2*xp*yp + p1*(r^2 + 2*yp^2));

Где

r = sqrt(x^2 + y^2);

У вас есть 5 параметров: k1, k2, k3, p1, p2 для радиального и тангенциального искажения и 200 пар точек, так что вы можете решить нелинейную систему.

Убедитесь, что массивы x, y, xp и yp существуют в рабочей области, и объявите их глобальными:

global x y xp yp

Напишите функцию для оценки среднеквадратической ошибки с заданным набором произвольных коэффициентов искажения, скажем, она называется «dist»:

function val = dist(var)

global x y xp yp

val = zeros(size(xp));

k1 = var(1);
k2 = var(2);
k3 = var(3);
p1 = var(4);
p2 = var(5);

r = sqrt(xp.*xp + yp.*yp);
temp1 = x - ( xp.*(1 + k1*r.^2 + k2*r.^4 + k3*r.^6) + 2*p1*xp.*yp + p2*(r.^2 + 2*xp.^2));
temp2 = y - ( yp.*(1 + k1*r.^2 + k2*r.^4 + k3*r.^6) + 2*p2*xp.*yp + p1*(r.^2 + 2*yp.^2));
val = sqrt(temp1.*temp1 + temp2.*temp2);

Решите систему с помощью «fsolve»:

[coef, fval] = fsolve(@dist, zeros(5,1));

Значения в «coef» — это коэффициенты искажения, которые вы ищете. Чтобы исправить искажение новых точек (xp, yp), отсутствующих в исходном наборе, используйте уравнения:

r = sqrt(xp.*xp + yp.*yp);
x_corr = xp.*(1 + k1*r.^2 + k2*r.^4 + k3*r.^6) + 2*p1*xp.*yp + p2*(r.^2 + 2*xp.^2);
y_corr = yp.*(1 + k1*r.^2 + k2*r.^4 + k3*r.^6) + 2*p2*xp.*yp + p1*(r.^2 + 2*yp.^2);

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

Примечания:

• Координаты должны быть смещены к центру изображения, так как искажение симметрично относительно него.
• Нет необходимости нормировать до интервала [-1, 1], но обычно это делается для того, чтобы полученные коэффициенты искажения были более или менее одного порядка (работая со степенями 2, 4 и 6 пикселей). координатам потребуются очень маленькие коэффициенты).
• Этот метод не требует, чтобы точки на изображении находились в единой сетке.