# Experiment 1: 8-point radial distortion problem¶

The $$8$$-point problem for cameras with radial distortion consists of $$9$$ equations in $$9$$ unknowns. Eight equations are given by

$\begin{split}\begin{pmatrix} p_{i,1}' & p_{i,2}' & 1 + r_i' \lambda \end{pmatrix} \underbrace{ \begin{pmatrix} f_{1,1} & f_{1,2} & f_{1,3}\\ f_{2,1} & f_{2,2} & f_{2,3}\\ f_{3,1} & f_{3,2} & 1 \end{pmatrix} }_{\boldsymbol{F}} \begin{pmatrix} p_{i,1} \\ p_{i,2} \\ 1 + r_i \lambda \end{pmatrix} = 0,\end{split}$

where the parameters $$p_{i,j}, p_{i,j}'$$ and $$r_i, r_i'$$ are known and represent distorted image coordinates and distortion radii, respectively. The true image coordinates are known only after the radial distortion parameter $:nbsphinx-math:lambda $ is recovered. The matrix :math:boldsymbol{F} is called the fundamental matrix and satisfies the additional constraint

$\det \boldsymbol{F} = 0.$

We activate the project’s environment.

[1]:

using Pkg
Pkg.activate(".")

 Activating environment at ~/Documents/Projects/cox-homotopies/code_cox_homotopies_fix/Project.toml


We include our routines.

[2]:

include("CoxHomotopy.jl")

polymake version 4.2
Ewgenij Gawrilow, Michael Joswig, and the polymake team
Technische Universität Berlin, Germany
https://polymake.org

This is free software licensed under GPL; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


[2]:

trackPathRandom (generic function with 1 method)


Next, we generate the equations.

[9]:

@polyvar(f[1:3,1:3], λ)
F=[f[1,1] f[1,2] f[1,3];
f[2,1] f[2,2] f[2,3];
f[3,1] f[3,2] 1
]

function generateSystem(; degenerate=false, epsilon=1e-7)
r=[rand() for i=1:8]
p = [vcat(rand(2,1),[1+λ*r[i]]) for i=1:8]
q12 = [rand(2,1) for i=1:8]
if degenerate
s = [-q12[i][1]+epsilon  for i=1:8]
else
s = [rand() for i=1:8]
end
q = [vcat(q12[i],[1+λ*s[i]]) for i=1:8]
vcat(vcat([(p[i]'*F*q[i])[1] for i=1:8]),[det(F)])
end

eq = generateSystem(degenerate=true, epsilon=1e-12);


The following snippet computes the polytope information and the generic orbit degree for this example.

[10]:

(P,Ft,a)=computePolytope(eq)
deg, Pmtx = orbitDegree(Ft)
println("The orbit degree is " * string(numerator(deg)))

The orbit degree is 4583


Solve the system using the Cox homotopy

[11]:

@time homSols, f, toricSols, homTargetResult = coxHomotopy(eq, P, Ft, a, monodromy_initial_timeout=0.1);

Tracking 18 paths... 100%|██████████████████████████████| Time: 0:00:01
# paths tracked:                  18
# non-singular solutions (real):  0 (0)
# singular endpoints (real):      12 (0)
# total solutions (real):         12 (0)
22.682467 seconds (57.74 M allocations: 1.861 GiB, 1.47% gc time)


We display the residual of the computed solutions, both in their toric and in their homogeneous coordinates.

[12]:

display(sort(get_residual(f,homSols, variables(f))))
display(sort(get_residual(eq,toricSols, variables(eq))))

16-element Array{Float64,1}:
3.697661377234898e-17
6.035588450010074e-17
6.56030666247518e-17
2.94632073128461e-16
2.9701530353310713e-16
2.2484394287606995e-15
6.258412536601062e-15
9.31617754484122e-15
7.506515742460871e-14
4.973963198697082e-13
1.5177313135249e-12
4.531831332865626e-12
6.777031120905949e-12
8.82514180939293e-12
9.022590033975171e-12
3.659391017275076e-11

16-element Array{Float64,1}:
1.1138044336203681e-15
1.1676695804299181e-15
1.759358273457284e-15
2.1034234396723264e-15
2.7113410871299463e-15
3.861958742080553e-15
4.528733562633992e-15
6.2884095731014595e-15
2.757477504586345e-12
5.75185381442267e-12
5.808784240103956e-12
6.551398579258292e-12
9.967747869098144e-12
1.0827524069551687e-11
1.5627352795273583e-11
1.7991132390563422e-10

[ ]:



[ ]:



[ ]: