I've been trying to make Zhang's "Single-View Geometry of A Rectangle With Application to Whiteboard Image Rectification" work in Python.
Abstract and formula source from Microsoft
This came from another question here on the stack, but the code was written in Sage and no longer runs on the latest Sage version. And most of that thread was a debate about whether it was possible. It is, but making it work in Python has been a problem.
The Problem: The formula it uses to calculate the focal throws a math domain error because it tries to get a square root of a negative number. I'm not sure I've converted Zhang's formulas to Python correctly. Without the correct focal it isn't very accurate at determining the ratio of the length to width of the rectangle. With the correct focal it seems to work pretty well.
Edit: Here is the version that will run on the new Sage website. It doesn't throw an error for the frame, but it comes up with a focal and w/h ratio that are way off: https://cloud.sagemath.com/projects/a8af3fd5-1390-4e80-a12a-13f18f77fdc0/files/Sage Whiteboard rectification.sagews
Any help would be much appreciated. I've put many hours into this.
The following code is from Python 2.7
import numpy as np
from numpy import matrix
from numpy import *
from math import sqrt
from math import pow
def FindAspect():
# 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
# The shape in the image is a perfect square filmed with a focal of 960.
w = 1920
h = 1080
u = 960
v = 540
# These coords work, and give a ratio of 1.1499 using the focal formula
# If you input the focal as 960 then the second part correctly identifies
# the ratio very accurately. 1.0025
## m1x = 690.196
## m1y = 998.728
##
## m2x = 1112.58
## m2y = 798.383
##
## m3x = 690.201
## m3y = 81.2587
##
## m4x = 1112.6
## m4y = 281.54
#These coords throw an error due to focal formula sqrt of negative number.
m1x = 849.029
m1y = 792.363
m2x = 1170.147
m2y = 1019.212
m3x = 849.033
m3y = 287.516
m4x = 1170.154
m4y = 60.757
# 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 = u
v0 = v
# pixel aspect ratio; for a normal camera pixels are square, so s=1
s = 1
# homogenous coordinates of the quadrangle
m1 = array([m1x,m1y,1])
m2 = array([m2x,m2y,1])
m3 = array([m3x,m3y,1])
m4 = array([m4x,m4y,1])
# The following equations are later used in calculating the the focal length
# and the rectangle's aspect ratio.
# see [zhang-single] Equation 11, 12
k2 = np.dot(np.cross(m1, m4), m3) / np.dot(np.cross(m2, m4), m3)
k3 = np.dot(np.cross(m1, m4), m2) / np.dot(np.cross(m3, m4), m2)
# see [zhang-single] Equation 14,16
n2 = k2 * m2 - m1
n3 = k3 * m3 - m1
# the focal length of the camera.
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
)
)
)
print 'Focal: ', f
## f = 960
# standard pinhole camera matrix
# see [zhang-single] Equation 1
K = np.matrix([[f,0,u0],[0,f,v0],[0,0,1]])
#the width/height ratio of the original rectangle
# see [zhang-single] Equation 20
whRatio = sqrt(
(np.dot(np.asarray(n2 * (K.T**-1)).reshape(-1), np.asarray((K**-1).dot(n2)).reshape(-1))) /
(np.dot(np.asarray(n3 * (K.T**-1)).reshape(-1), np.asarray((K**-1).dot(n3)).reshape(-1)))
)
print whRatio
FindAspect()