Пропорции перспективно деформированного прямоугольника
Учитывая 2d изображение прямоугольника, искаженного в перспективе:
Я знаю, что форма изначально была прямоугольником, но я не знаю ее первоначальный размер.
Если я знаю пиксельные координаты углов на этом рисунке, как я могу рассчитать исходные пропорции, то есть частное (ширина / высота) прямоугольника?
(фон: цель состоит в том, чтобы автоматически не искажать фотографии прямоугольных документов, обнаружение края, вероятно, будет выполнено с помощью грубого преобразования)
ОБНОВИТЬ:
Было некоторое обсуждение того, возможно ли вообще определить соотношение ширина: высота с предоставленной информацией. Моя наивная мысль заключалась в том, что это должно быть возможно, так как я не могу придумать, как спроецировать, например, прямоугольник 1:4 на четырехугольник, изображенный выше. Соотношение явно близко к 1:1, поэтому должен быть способ определить его математически. У меня, однако, нет никаких доказательств этого, кроме моего интуитивного предположения.
Я еще не полностью понял аргументы, представленные ниже, но я думаю, что должно быть какое-то неявное предположение, что мы здесь упускаем, и это интерпретируется по-другому.
Однако после нескольких часов поисков я наконец нашел несколько работ, имеющих отношение к проблеме. Я изо всех сил пытаюсь понять математику, используемую там, пока безуспешно. В частности, первая статья, кажется, обсуждает именно то, что я хотел сделать, к сожалению, без примеров кода и очень плотной математики.
Чжэнъю Чжан, Ли-Вэй Хэ, "Сканирование белой доски и улучшение изображения" http://research.microsoft.com/en-us/um/people/zhang/papers/tr03-39.pdf стр.11
"Из-за искажения перспективы изображение прямоугольника выглядит как четырехугольник. Однако, поскольку мы знаем, что это прямоугольник в пространстве, мы можем оценить как фокусное расстояние камеры, так и соотношение сторон прямоугольника".
РОБЕРТ М. ХАРАЛИК. "Определение параметров камеры по перспективной проекции прямоугольника" http://portal.acm.org/citation.cfm?id=87146
"Мы показываем, как использовать двухмерную перспективную проекцию прямоугольника неизвестного размера и положения в трехмерном пространстве для определения параметров угла обзора камеры относительно планов прямоугольника".
11 ответов
Вот моя попытка ответить на мой вопрос после прочтения газеты
- Чжэнъю Чжан, Ли-Вэй Хэ, "Сканирование белой доски и улучшение изображения" http://research.microsoft.com/en-us/um/people/zhang/papers/tr03-39.pdf
Некоторое время я манипулировал уравнениями в SAGE и придумал этот псевдокод в стиле c:
// in case it matters: licensed under GPLv2 or later
// legend:
// sqr(x) = x*x
// sqrt(x) = square root of x
// let m1x,m1y ... m4x,m4y be the (x,y) pixel coordinates
// of the 4 corners of the detected quadrangle
// i.e. (m1x, m1y) are the cordinates of the first corner,
// (m2x, m2y) of the second corner and so on.
// let u0, v0 be the pixel coordinates of the principal point of the image
// for a normal camera this will be the center of the image,
// i.e. u0=IMAGEWIDTH/2; v0 =IMAGEHEIGHT/2
// This assumption does not hold if the image has been cropped asymmetrically
// first, transform the image so the principal point is at (0,0)
// this makes the following equations much easier
m1x = m1x - u0;
m1y = m1y - v0;
m2x = m2x - u0;
m2y = m2y - v0;
m3x = m3x - u0;
m3y = m3y - v0;
m4x = m4x - u0;
m4y = m4y - v0;
// temporary variables k2, k3
double k2 = ((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y - m1y*m4x) /
((m2y - m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x) ;
double k3 = ((m1y - m4y)*m2x - (m1x - m4x)*m2y + m1x*m4y - m1y*m4x) /
((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y - m3y*m4x) ;
// f_squared is the focal length of the camera, squared
// if k2==1 OR k3==1 then this equation is not solvable
// if the focal length is known, then this equation is not needed
// in that case assign f_squared= sqr(focal_length)
double f_squared =
-((k3*m3y - m1y)*(k2*m2y - m1y) + (k3*m3x - m1x)*(k2*m2x - m1x)) /
((k3 - 1)*(k2 - 1)) ;
//The width/height ratio of the original rectangle
double whRatio = sqrt(
(sqr(k2 - 1) + sqr(k2*m2y - m1y)/f_squared + sqr(k2*m2x - m1x)/f_squared) /
(sqr(k3 - 1) + sqr(k3*m3y - m1y)/f_squared + sqr(k3*m3x - m1x)/f_squared)
) ;
// if k2==1 AND k3==1, then the focal length equation is not solvable
// but the focal length is not needed to calculate the ratio.
// I am still trying to figure out under which circumstances k2 and k3 become 1
// but it seems to be when the rectangle is not distorted by perspective,
// i.e. viewed straight on. Then the equation is obvious:
if (k2==1 && k3==1) whRatio = sqrt(
(sqr(m2y-m1y) + sqr(m2x-m1x)) /
(sqr(m3y-m1y) + sqr(m3x-m1x))
// After testing, I found that the above equations
// actually give the height/width ratio of the rectangle,
// not the width/height ratio.
// If someone can find the error that caused this,
// I would be most grateful.
// until then:
whRatio = 1/whRatio;
Обновление: вот как эти уравнения были определены:
Ниже приведен код в SAGE. Доступ к нему можно получить через Интернет по адресу http://www.sagenb.org/home/pub/704/. (Sage действительно полезен в решении уравнений, и его можно использовать в любом браузере, проверьте это)
# CALCULATING THE ASPECT RATIO OF A RECTANGLE DISTORTED BY PERSPECTIVE
#
# BIBLIOGRAPHY:
# [zhang-single]: "Single-View Geometry of A Rectangle
# With Application to Whiteboard Image Rectification"
# by Zhenggyou Zhang
# http://research.microsoft.com/users/zhang/Papers/WhiteboardRectification.pdf
# pixel coordinates of the 4 corners of the quadrangle (m1, m2, m3, m4)
# see [zhang-single] figure 1
m1x = var('m1x')
m1y = var('m1y')
m2x = var('m2x')
m2y = var('m2y')
m3x = var('m3x')
m3y = var('m3y')
m4x = var('m4x')
m4y = var('m4y')
# pixel coordinates of the principal point of the image
# for a normal camera this will be the center of the image,
# i.e. u0=IMAGEWIDTH/2; v0 =IMAGEHEIGHT/2
# This assumption does not hold if the image has been cropped asymmetrically
u0 = var('u0')
v0 = var('v0')
# pixel aspect ratio; for a normal camera pixels are square, so s=1
s = var('s')
# homogenous coordinates of the quadrangle
m1 = vector ([m1x,m1y,1])
m2 = vector ([m2x,m2y,1])
m3 = vector ([m3x,m3y,1])
m4 = vector ([m4x,m4y,1])
# the following equations are later used in calculating the the focal length
# and the rectangle's aspect ratio.
# temporary variables: k2, k3, n2, n3
# see [zhang-single] Equation 11, 12
k2_ = m1.cross_product(m4).dot_product(m3) / m2.cross_product(m4).dot_product(m3)
k3_ = m1.cross_product(m4).dot_product(m2) / m3.cross_product(m4).dot_product(m2)
k2 = var('k2')
k3 = var('k3')
# see [zhang-single] Equation 14,16
n2 = k2 * m2 - m1
n3 = k3 * m3 - m1
# the focal length of the camera.
f = var('f')
# see [zhang-single] Equation 21
f_ = sqrt(
-1 / (
n2[2]*n3[2]*s^2
) * (
(
n2[0]*n3[0] - (n2[0]*n3[2]+n2[2]*n3[0])*u0 + n2[2]*n3[2]*u0^2
)*s^2 + (
n2[1]*n3[1] - (n2[1]*n3[2]+n2[2]*n3[1])*v0 + n2[2]*n3[2]*v0^2
)
)
)
# standard pinhole camera matrix
# see [zhang-single] Equation 1
A = matrix([[f,0,u0],[0,s*f,v0],[0,0,1]])
#the width/height ratio of the original rectangle
# see [zhang-single] Equation 20
whRatio = sqrt (
(n2*A.transpose()^(-1) * A^(-1)*n2.transpose()) /
(n3*A.transpose()^(-1) * A^(-1)*n3.transpose())
)
Упрощенные уравнения в с-коде, где определяется
print "simplified equations, assuming u0=0, v0=0, s=1"
print "k2 := ", k2_
print "k3 := ", k3_
print "f := ", f_(u0=0,v0=0,s=1)
print "whRatio := ", whRatio(u0=0,v0=0,s=1)
simplified equations, assuming u0=0, v0=0, s=1
k2 := ((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y - m1y*m4x)/((m2y
- m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x)
k3 := ((m1y - m4y)*m2x - (m1x - m4x)*m2y + m1x*m4y - m1y*m4x)/((m3y
- m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y - m3y*m4x)
f := sqrt(-((k3*m3y - m1y)*(k2*m2y - m1y) + (k3*m3x - m1x)*(k2*m2x
- m1x))/((k3 - 1)*(k2 - 1)))
whRatio := sqrt(((k2 - 1)^2 + (k2*m2y - m1y)^2/f^2 + (k2*m2x -
m1x)^2/f^2)/((k3 - 1)^2 + (k3*m3y - m1y)^2/f^2 + (k3*m3x -
m1x)^2/f^2))
print "Everything in one equation:"
print "whRatio := ", whRatio(f=f_)(k2=k2_,k3=k3_)(u0=0,v0=0,s=1)
Everything in one equation:
whRatio := sqrt(((((m1y - m4y)*m2x - (m1x - m4x)*m2y + m1x*m4y -
m1y*m4x)/((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y - m3y*m4x) -
1)*(((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y - m1y*m4x)/((m2y -
m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x) - 1)*(((m1y -
m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y - m1y*m4x)*m2y/((m2y - m4y)*m3x
- (m2x - m4x)*m3y + m2x*m4y - m2y*m4x) - m1y)^2/((((m1y - m4y)*m2x -
(m1x - m4x)*m2y + m1x*m4y - m1y*m4x)*m3y/((m3y - m4y)*m2x - (m3x -
m4x)*m2y + m3x*m4y - m3y*m4x) - m1y)*(((m1y - m4y)*m3x - (m1x -
m4x)*m3y + m1x*m4y - m1y*m4x)*m2y/((m2y - m4y)*m3x - (m2x - m4x)*m3y
+ m2x*m4y - m2y*m4x) - m1y) + (((m1y - m4y)*m2x - (m1x - m4x)*m2y +
m1x*m4y - m1y*m4x)*m3x/((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y
- m3y*m4x) - m1x)*(((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y -
m1y*m4x)*m2x/((m2y - m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x)
- m1x)) + (((m1y - m4y)*m2x - (m1x - m4x)*m2y + m1x*m4y -
m1y*m4x)/((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y - m3y*m4x) -
1)*(((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y - m1y*m4x)/((m2y -
m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x) - 1)*(((m1y -
m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y - m1y*m4x)*m2x/((m2y - m4y)*m3x
- (m2x - m4x)*m3y + m2x*m4y - m2y*m4x) - m1x)^2/((((m1y - m4y)*m2x -
(m1x - m4x)*m2y + m1x*m4y - m1y*m4x)*m3y/((m3y - m4y)*m2x - (m3x -
m4x)*m2y + m3x*m4y - m3y*m4x) - m1y)*(((m1y - m4y)*m3x - (m1x -
m4x)*m3y + m1x*m4y - m1y*m4x)*m2y/((m2y - m4y)*m3x - (m2x - m4x)*m3y
+ m2x*m4y - m2y*m4x) - m1y) + (((m1y - m4y)*m2x - (m1x - m4x)*m2y +
m1x*m4y - m1y*m4x)*m3x/((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y
- m3y*m4x) - m1x)*(((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y -
m1y*m4x)*m2x/((m2y - m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x)
- m1x)) - (((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y -
m1y*m4x)/((m2y - m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x) -
1)^2)/((((m1y - m4y)*m2x - (m1x - m4x)*m2y + m1x*m4y -
m1y*m4x)/((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y - m3y*m4x) -
1)*(((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y - m1y*m4x)/((m2y -
m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x) - 1)*(((m1y -
m4y)*m2x - (m1x - m4x)*m2y + m1x*m4y - m1y*m4x)*m3y/((m3y - m4y)*m2x
- (m3x - m4x)*m2y + m3x*m4y - m3y*m4x) - m1y)^2/((((m1y - m4y)*m2x -
(m1x - m4x)*m2y + m1x*m4y - m1y*m4x)*m3y/((m3y - m4y)*m2x - (m3x -
m4x)*m2y + m3x*m4y - m3y*m4x) - m1y)*(((m1y - m4y)*m3x - (m1x -
m4x)*m3y + m1x*m4y - m1y*m4x)*m2y/((m2y - m4y)*m3x - (m2x - m4x)*m3y
+ m2x*m4y - m2y*m4x) - m1y) + (((m1y - m4y)*m2x - (m1x - m4x)*m2y +
m1x*m4y - m1y*m4x)*m3x/((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y
- m3y*m4x) - m1x)*(((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y -
m1y*m4x)*m2x/((m2y - m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x)
- m1x)) + (((m1y - m4y)*m2x - (m1x - m4x)*m2y + m1x*m4y -
m1y*m4x)/((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y - m3y*m4x) -
1)*(((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y - m1y*m4x)/((m2y -
m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x) - 1)*(((m1y -
m4y)*m2x - (m1x - m4x)*m2y + m1x*m4y - m1y*m4x)*m3x/((m3y - m4y)*m2x
- (m3x - m4x)*m2y + m3x*m4y - m3y*m4x) - m1x)^2/((((m1y - m4y)*m2x -
(m1x - m4x)*m2y + m1x*m4y - m1y*m4x)*m3y/((m3y - m4y)*m2x - (m3x -
m4x)*m2y + m3x*m4y - m3y*m4x) - m1y)*(((m1y - m4y)*m3x - (m1x -
m4x)*m3y + m1x*m4y - m1y*m4x)*m2y/((m2y - m4y)*m3x - (m2x - m4x)*m3y
+ m2x*m4y - m2y*m4x) - m1y) + (((m1y - m4y)*m2x - (m1x - m4x)*m2y +
m1x*m4y - m1y*m4x)*m3x/((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y
- m3y*m4x) - m1x)*(((m1y - m4y)*m3x - (m1x - m4x)*m3y + m1x*m4y -
m1y*m4x)*m2x/((m2y - m4y)*m3x - (m2x - m4x)*m3y + m2x*m4y - m2y*m4x)
- m1x)) - (((m1y - m4y)*m2x - (m1x - m4x)*m2y + m1x*m4y -
m1y*m4x)/((m3y - m4y)*m2x - (m3x - m4x)*m2y + m3x*m4y - m3y*m4x) -
1)^2))
# some testing:
# - choose a random rectangle,
# - project it onto a random plane,
# - insert the corners in the above equations,
# - check if the aspect ratio is correct.
from sage.plot.plot3d.transform import rotate_arbitrary
#redundandly random rotation matrix
rand_rotMatrix = \
rotate_arbitrary((uniform(-5,5),uniform(-5,5),uniform(-5,5)),uniform(-5,5)) *\
rotate_arbitrary((uniform(-5,5),uniform(-5,5),uniform(-5,5)),uniform(-5,5)) *\
rotate_arbitrary((uniform(-5,5),uniform(-5,5),uniform(-5,5)),uniform(-5,5))
#random translation vector
rand_transVector = vector((uniform(-10,10),uniform(-10,10),uniform(-10,10))).transpose()
#random rectangle parameters
rand_width =uniform(0.1,10)
rand_height=uniform(0.1,10)
rand_left =uniform(-10,10)
rand_top =uniform(-10,10)
#random focal length and principal point
rand_f = uniform(0.1,100)
rand_u0 = uniform(-100,100)
rand_v0 = uniform(-100,100)
# homogenous standard pinhole projection, see [zhang-single] Equation 1
hom_projection = A * rand_rotMatrix.augment(rand_transVector)
# construct a random rectangle in the plane z=0, then project it randomly
rand_m1hom = hom_projection*vector((rand_left ,rand_top ,0,1)).transpose()
rand_m2hom = hom_projection*vector((rand_left ,rand_top+rand_height,0,1)).transpose()
rand_m3hom = hom_projection*vector((rand_left+rand_width,rand_top ,0,1)).transpose()
rand_m4hom = hom_projection*vector((rand_left+rand_width,rand_top+rand_height,0,1)).transpose()
#change type from 1x3 matrix to vector
rand_m1hom = rand_m1hom.column(0)
rand_m2hom = rand_m2hom.column(0)
rand_m3hom = rand_m3hom.column(0)
rand_m4hom = rand_m4hom.column(0)
#normalize
rand_m1hom = rand_m1hom/rand_m1hom[2]
rand_m2hom = rand_m2hom/rand_m2hom[2]
rand_m3hom = rand_m3hom/rand_m3hom[2]
rand_m4hom = rand_m4hom/rand_m4hom[2]
#substitute random values for f, u0, v0
rand_m1hom = rand_m1hom(f=rand_f,s=1,u0=rand_u0,v0=rand_v0)
rand_m2hom = rand_m2hom(f=rand_f,s=1,u0=rand_u0,v0=rand_v0)
rand_m3hom = rand_m3hom(f=rand_f,s=1,u0=rand_u0,v0=rand_v0)
rand_m4hom = rand_m4hom(f=rand_f,s=1,u0=rand_u0,v0=rand_v0)
# printing the randomly choosen values
print "ground truth: f=", rand_f, "; ratio=", rand_width/rand_height
# substitute all the variables in the equations:
print "calculated: f= ",\
f_(k2=k2_,k3=k3_)(s=1,u0=rand_u0,v0=rand_v0)(
m1x=rand_m1hom[0],m1y=rand_m1hom[1],
m2x=rand_m2hom[0],m2y=rand_m2hom[1],
m3x=rand_m3hom[0],m3y=rand_m3hom[1],
m4x=rand_m4hom[0],m4y=rand_m4hom[1],
),"; 1/ratio=", \
1/whRatio(f=f_)(k2=k2_,k3=k3_)(s=1,u0=rand_u0,v0=rand_v0)(
m1x=rand_m1hom[0],m1y=rand_m1hom[1],
m2x=rand_m2hom[0],m2y=rand_m2hom[1],
m3x=rand_m3hom[0],m3y=rand_m3hom[1],
m4x=rand_m4hom[0],m4y=rand_m4hom[1],
)
print "k2 = ", k2_(
m1x=rand_m1hom[0],m1y=rand_m1hom[1],
m2x=rand_m2hom[0],m2y=rand_m2hom[1],
m3x=rand_m3hom[0],m3y=rand_m3hom[1],
m4x=rand_m4hom[0],m4y=rand_m4hom[1],
), "; k3 = ", k3_(
m1x=rand_m1hom[0],m1y=rand_m1hom[1],
m2x=rand_m2hom[0],m2y=rand_m2hom[1],
m3x=rand_m3hom[0],m3y=rand_m3hom[1],
m4x=rand_m4hom[0],m4y=rand_m4hom[1],
)
# ATTENTION: testing revealed, that the whRatio
# is actually the height/width ratio,
# not the width/height ratio
# This contradicts [zhang-single]
# if anyone can find the error that caused this, I'd be grateful
ground truth: f= 72.1045134124554 ; ratio= 3.46538779959142
calculated: f= 72.1045134125 ; 1/ratio= 3.46538779959
k2 = 0.99114614987 ; k3 = 1.57376280159
Обновить
После прочтения вашего обновления и просмотра первой справки (сканирование белой доски и улучшение изображений) я вижу, где находится пропущенная точка.
Входные данные задачи - это четверка (A,B,C,D) и центр O проецируемого изображения. В статье это соответствует предположению u0=v0=0. Добавляя эту точку, проблема становится достаточно ограниченной, чтобы получить соотношение сторон прямоугольника.
Затем проблема переформулируется следующим образом: для заданной четверки (A,B,C,D) в плоскости Z=0 найдите положение глаза E(0,0,h), h>0 и трехмерную плоскость P, такую что проекция (A,B,C,D) на P является прямоугольником.
Обратите внимание, что P определяется E: чтобы получить параллелограмм, P должен содержать параллели с (EU) и (EV), где U=(AB)x(CD) и V=(AD)x(BC).
Экспериментально кажется, что эта проблема имеет в общем одно уникальное решение, соответствующее уникальному значению отношения w/h прямоугольника.
Предыдущий пост
Нет, вы не можете определить соотношение прямоугольников из проекции.
В общем случае четверка (A,B,C,D) четырех неколлинеарных точек плоскости Z=0 является проекцией бесконечно большого числа прямоугольников с бесконечно большим числом соотношений ширина / высота.
Рассмотрим две точки схода U, пересечение (AB) и (CD) и V, пересечение (AD) и (BC) и точку I, пересечение двух диагоналей (AC) и (BD). Чтобы проецировать как ABCD, параллелограмм центра I должен лежать на плоскости, содержащей линию, параллельную (UV) через точку I. На одной такой плоскости вы можете найти много прямоугольников, выступающих на ABCD, все с различным отношением w/h.
Посмотрите эти два изображения, сделанные с помощью Cabri 3D. В обоих случаях ABCD не изменяется (на серой плоскости Z=0), и синяя плоскость, содержащая прямоугольник, также не изменяется. Зеленая линия, частично скрытая, является (УФ) линией, а видимая зеленая линия параллельна ей и содержит I.
По этому ответу вы можете определить ширину / высоту. Вычислить трехмерную координату прямоугольника с координатой его тени?, Предположим, что ваш прямоугольник вращается на диагональной точке пересечения, рассчитайте его ширину и высоту. Но когда вы меняете расстояние между предполагаемой теневой плоскостью и реальной теневой плоскостью, пропорциональной прямоугольнику, то же самое с вычисленной шириной / высотой!
Размер на самом деле не нужен, и не пропорции. И знать, какая сторона вверх, является неуместным, учитывая, что он использует фотографии / сканы документов. Я сомневаюсь, что он собирается сканировать их задние стороны.
"Угловое пересечение" - это метод исправления перспективы. Это может быть полезно:
На вопрос о том, почему результаты дают ч / ш, а не ш / ч: мне интересно, правильно ли приведенное выше уравнение 20. Опубликовано это:
whRatio = sqrt (
(n2*A.transpose()^(-1) * A^(-1)*n2.transpose()) /
(n3*A.transpose()^(-1) * A^(-1)*n3.transpose())
)
Когда я пытаюсь выполнить это с OpenCV, я получаю исключение. Но все работает правильно, когда я использую следующее уравнение, которое для меня больше похоже на уравнение 20: Но на основе уравнения 20 это выглядит так, как должно быть:
whRatio = sqrt (
(n2.transpose()*A.transpose()^(-1) * A^(-1)*n2) /
(n3.transpose()*A.transpose()^(-1) * A^(-1)*n3)
)
Вот реализация Python https://www.microsoft.com/en-us/research/publication/2016/11/Digital-Signal-Processing.pdf .
import numpy as np
def find_rectangle_proportion(K, m1, m2, m3, m4):
# Detailed computations
# https://www.microsoft.com/en-us/research/publication/2016/11/Digital-Signal-Processing.pdf
# m1 = top_left, m2 = top_right, m3 = bottom_left, m4 = bottom_right
# Make homeneous coordinates
m1 = np.array([m1[0], m1[1], 1])
m2 = np.array([m2[0], m2[1], 1])
m3 = np.array([m3[0], m3[1], 1])
m4 = np.array([m4[0], m4[1], 1])
# (11)
k2 = np.dot(np.cross(m1, m4), m3) / np.dot(np.cross(m2, m4), m3)
# (12)
k3 = np.dot(np.cross(m1, m4), m2) / np.dot(np.cross(m3, m4), m2)
# (14)
n2 = k2 * m2 - m1
# (16)
n3 = k3 * m3 - m1
inv_K = np.linalg.inv(K)
KK = np.dot(inv_K.T, inv_K)
# ratio width/height (20)
ratio = np.sqrt(np.dot(n2, np.dot(KK, n2)) / np.dot(n3, np.dot(KK, n3)))
return ratio
if __name__ == "__main__":
top_left = [269, 25]
top_right = [800, 19]
bottom_right = [805, 748]
bottom_left = [273, 750]
# Load intrinsic matrices from calib.py
K = np.load("calibration_pixel_6a/intrinsic.npy")
ratio = find_rectangle_proportion(K, top_left, top_right, bottom_left, bottom_right)
print(f"ratio width/height = {ratio}")
Кажется, все еще есть некоторая путаница в этой интересной проблеме. Я хочу дать простое объяснение того, когда проблема может и не может быть решена.
Ограничения и степени свободы
Обычно, когда мы сталкиваемся с такой проблемой, первое, что нужно сделать, это оценить количество неизвестных степеней свободы (DoFs) N и количество независимых уравнений M, которые у нас есть для ограничения неизвестных степеней свободы. Невозможно решить проблему, если N превышает M (это означает, что ограничений меньше, чем неизвестных). Мы можем исключить все проблемы, которые в этом случае являются неразрешимыми. Если N не превышает M, то можно решить проблему с помощью уникального решения, но это не гарантировано (см. Пример со второго по последний абзац).
Давайте использовать p 1, p 2, p 3 и p 4, чтобы обозначить положения 4 углов плоской поверхности в мировых координатах. Давайте использовать R и t, чтобы быть трехмерным вращением и переводом, который преобразует их в координаты камеры. Давайте использовать K для обозначения внутренней матрицы камеры 3x3. Мы пока проигнорируем искажение линзы. 2D-положение i- го угла в изображении с камеры определяется как q i = f (K (Rp i + t)), где f - функция проекции f(x,y,z)=(x/z,y/z). Используя это уравнение, мы знаем, что каждый угол изображения дает нам два уравнения (т.е. два ограничения) для наших неизвестных: одно из компонента x в q i и одно из компонента y. Таким образом, у нас есть 8 ограничений для работы. Официальное название этих ограничений - ограничения репроекции.
Итак, каковы наши неизвестные DoFs? Конечно, R и t неизвестны, потому что мы не знаем положение камеры в мировых координатах. Поэтому у нас уже есть 6 неизвестных степеней свободы: 3 для R (например, рыскание, тангаж и крен) и 3 для t. Поэтому в оставшихся слагаемых может быть максимум двух неизвестных (K, p 1, p 2, p 3, p 4).
Разные проблемы
Мы можем построить различные задачи в зависимости от того, какие два члена в (K, p 1, p 2, p 3, p 4) мы будем считать неизвестными. На этом этапе выпишем K в обычном виде: K = (fx, 0, cx; 0, fy, cy; 0,0,1) где fx и fy - члены фокусного расстояния (fx/fy обычно называют соотношение сторон изображения) и (cx,cy) является главной точкой (центр проекции на изображении).
Мы могли бы получить одну проблему, используя fx и fy в качестве двух наших неизвестных, и предположить, что (cx, cy, p 1, p 2, p 3, p 4) все известны. Действительно, именно эта проблема используется и решается в рамках метода калибровки камеры OpenCV, используя изображения плоской цели шахматной доски. Это используется, чтобы получить начальную оценку для fx и fy, предполагая, что главная точка находится в центре изображения (что является очень разумным допущением для большинства камер).
В качестве альтернативы мы можем создать другую проблему, предполагая fx=fy, что опять-таки вполне разумно для многих камер, и предположим, что это фокусное расстояние (обозначаемое как f) является единственным неизвестным в K. Поэтому у нас все еще есть одно неизвестное, с которым можно играть (напомним, у нас может быть максимум два неизвестных). Итак, давайте использовать это, предположив, что мы знаем форму плоскости: как прямоугольник (что было первоначальным предположением в вопросе). Поэтому мы можем определить углы следующим образом: p 1 = (0,0,0), p 2 = (0, w, 0), p 3 = (h, 0,0) и p 4 = (h, w, 0), где h и w обозначают высоту и ширину прямоугольника. Теперь, поскольку у нас осталось только 1 неизвестное, давайте установим это как соотношение сторон плоскости: x=w/h. Теперь вопрос в том, можем ли мы одновременно восстановить x, f, R и t из 8 ограничений репроекции? Ответ оказывается да! И решение дано в статье Чжана, процитированной в вопросе.
Неоднозначность шкалы
Можно задаться вопросом, может ли быть решена другая проблема: если мы предположим, что K известно, а 2 неизвестных - это h и w. Могут ли они быть решены из уравнений репроекции? Ответ - нет, и потому, что существует неопределенность между размером плоскости и глубиной плоскости до камеры. В частности, если мы масштабируем углы p i по s и масштабируем t по s, то s отменяется в уравнениях перепроецирования. Поэтому абсолютный масштаб самолета не подлежит восстановлению.
Могут быть другие проблемы с различными комбинациями для неизвестных DoF, например, имеющих R, t, один из компонентов главной точки и ширину плоскости в качестве неизвестных. Однако необходимо подумать о том, какие случаи имеют практическое применение. Тем не менее, я еще не видел систематического набора решений для всех полезных комбинаций!
Больше очков
Мы можем подумать, что если бы мы добавили дополнительные точечные соответствия между плоскостью и изображением или использовали края плоскости, мы могли бы восстановить более 8 неизвестных степеней свободы. К сожалению, ответ - нет. Это потому, что они не добавляют никаких дополнительных независимых ограничений. Причина в том, что 4 угла полностью описывают преобразование из плоскости в изображение. Это можно увидеть, подгоняя матрицу гомографии, используя четыре угла, которые затем могут определять положения всех других точек на плоскости на изображении.
Невозможно узнать ширину этого прямоугольника, не зная расстояние "камеры".
маленький прямоугольник, видимый с расстояния в 5 сантиметров, выглядит так же, как огромный прямоугольник, если смотреть с метров
Нарисуйте правильный равнобедренный треугольник с этими двумя точками схода и третьей точкой под горизонтом (то есть на той же стороне горизонта, что и прямоугольник). Эта третья точка будет нашим источником, а две линии к точкам схода будут нашими осями. Назовите расстояние от начала координат до точки схода pi/2. Теперь продлите стороны прямоугольника от точек схода до осей и отметьте, где они пересекают оси. Выберите ось, измерьте расстояния от двух меток до начала координат, преобразуйте эти расстояния: x->tan(x), и разница будет "истинной" длиной этой стороны. Сделайте то же самое для другой оси. Возьмите соотношение этих двух длин, и все готово.
Dropbox имеет обширную статью в своем техническом блоге, где они описывают, как они решили проблему для своего приложения сканера.
https://blogs.dropbox.com/tech/2016/08/fast-document-rectification-and-enhancement/
Исправление документа
Мы предполагаем, что входной документ является прямоугольным в физическом мире, но если он не совсем направлен на камеру, получающиеся углы на изображении будут обычным выпуклым четырехугольником. Таким образом, чтобы удовлетворить нашу первую цель, мы должны отменить геометрическое преобразование, примененное процессом захвата. Это преобразование зависит от точки обзора камеры относительно документа (это так называемые внешние параметры), в дополнение к таким вещам, как фокусное расстояние камеры (внутренние параметры). Вот схема сценария захвата:
Чтобы отменить геометрическое преобразование, мы должны сначала определить упомянутые параметры. Если мы предположим, что камера хорошо симметрична (без астигматизма, без перекосов и т. Д.), Неизвестными в этой модели являются:
- 3D-расположение камеры относительно документа (3 степени свободы),
- 3D-ориентация камеры относительно документа (3 степени свободы),
- размеры документа (2 степени свободы) и
- фокусное расстояние камеры (1 степень свободы).
С другой стороны, x- и y-координаты четырех обнаруженных углов документа дают нам фактически восемь ограничений. Хотя на первый взгляд неизвестных (9) больше, чем ограничений (8), неизвестные не являются полностью свободными переменными - можно представить масштабирование документа физически и размещение его дальше от камеры, чтобы получить идентичную фотографию. Это отношение накладывает дополнительное ограничение, поэтому мы должны решить полностью ограниченную систему. (Реальная система уравнений, которую мы решаем, включает в себя несколько других соображений; соответствующая статья в Википедии дает хорошее резюме: https://en.wikipedia.org/wiki/Camera_resectioning)
Как только параметры были восстановлены, мы можем отменить геометрическое преобразование, примененное процессом захвата, чтобы получить хорошее прямоугольное изображение. Однако это потенциально трудоемкий процесс: для каждого выходного пикселя можно было бы посмотреть значение соответствующего входного пикселя в исходном изображении. Конечно, графические процессоры специально предназначены для таких задач: рендеринг текстуры в виртуальном пространстве. Существует преобразование вида - которое, как оказалось, является обратным преобразованию камеры, для которого мы только что решили! - с помощью которого можно отобразить полное входное изображение и получить исправленный документ. (Простой способ убедиться в этом - заметить, что после того, как на экране вашего телефона появится полное изображение ввода, вы можете наклонить и перевести телефон так, чтобы проекция области документа на экране была для вас прямолинейной.)
Наконец, напомним, что в отношении масштаба существовала неоднозначность: мы не можем сказать, был ли документ, например, размером с бумагу (8,5 х 11 дюймов) или плакатом (17 х 22 дюйма). Какими должны быть размеры выходного изображения? Чтобы устранить эту неоднозначность, мы подсчитываем количество пикселей в четырехугольнике во входном изображении и устанавливаем выходное разрешение в соответствии с этим количеством пикселей. Идея состоит в том, что мы не хотим слишком сильно увеличивать или уменьшать изображение.
Вам нужно больше информации, что преобразованная фигура может быть получена из любого параллелограмма с произвольной точки зрения.
Итак, я думаю, вам нужно сначала сделать какую-то калибровку.
Редактировать: для тех, кто сказал, что я был неправ, здесь приведено математическое доказательство того, что существуют бесконечные комбинации прямоугольников / камер, которые дают одну и ту же проекцию:
Чтобы упростить задачу (поскольку нам нужно только соотношение сторон), давайте предположим, что наш прямоугольник определяется следующими точками: R=[(0,0),(1,0),(1,r),(0,r)] (это упрощение аналогично преобразованию любой задачи в эквивалентную в аффинном пространстве).
Преобразованный многоугольник определяется как: T=[(tx0,ty0),(tx1,ty1),(tx2,ty2),(tx3,ty3)]
Существует матрица преобразования M = [[m00,m01,m02],[m10,m11,m12],[m20,m21,m22]] это удовлетворяет (Rxi,Ryi,1)*M=wi(txi,tyi,1)'
если мы расширим уравнение выше для точек,
за R_0 мы получаем: m02-tx0*w0 = m12-ty0*w0 = m22-w0 = 0
за R_1 мы получаем: m00-tx1*w1 = m10-ty1*w1 = m20+m22-w1 = 0
за R_2 мы получаем: m00+r*m01-tx2*w2 = m10+r*m11-ty2*w2 = m20+r*m21+m22-w2 = 0
и для R_3 мы получаем: m00+r*m01-tx3*w3 = m10+r*m11-ty3*w3 = m20 + r*m21 + m22 -w3 = 0
Пока у нас есть 12 уравнений, 14 неизвестных переменных (9 из матрицы, 4 из wiи 1 для соотношения r) а остальные являются известными значениями (txi а также tyi дано).
Даже если система не была занижена, некоторые из неизвестных умножаются между собой (r а также mi0 продукты) сделать систему нелинейной (вы можете преобразовать ее в линейную систему, присваивая новое имя каждому продукту, но у вас все равно останется 13 неизвестных, и 3 из них будут расширены до бесконечных решений).
Если вы можете найти какой-либо недостаток в рассуждениях или математике, пожалуйста, дайте мне знать.