# 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 first include our routines.

:

include("CoxHomotopy.jl")

:

trackPathRandom (generic function with 1 method)


Next, we generate the equations.

:

@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]+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]) for i=1:8]),[det(F)])
end

eq = generateSystem(degenerate=true);


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

:

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

The orbit degree is 4583

polymake: used package cdd
cddlib
Implementation of the double description method of Motzkin et al.
Copyright by Komei Fukuda.
http://www-oldurls.inf.ethz.ch/personal/fukudak/cdd_home/

polymake: used package ppl
The Parma Polyhedra Library ([[wiki:external_software#PPL]]): A C++ library for convex polyhedra
and other numerical abstractions.
http://www.cs.unipr.it/ppl/



Solve the system using the Cox homotopy

:

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

Tracking 6 paths... 100%|███████████████████████████████| Time: 0:00:00
# paths tracked:                  6
# non-singular solutions (real):  0 (0)
# singular endpoints (real):      2 (0)
# total solutions (real):         2 (0)
28.075607 seconds (100.54 M allocations: 3.724 GiB, 2.24% gc time)


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

:

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

16-element Array{Float64,1}:
4.723489140580794e-16
5.066357327780187e-16
8.325922430436857e-16
1.4612328645599375e-15
1.5312097326723611e-15
2.3637863043829518e-15
2.756752844177039e-15
9.60680008036377e-15
1.3927259770030184e-13
1.517291858886625e-13
2.7292206307646026e-13
6.385300390329031e-13
1.420555478198795e-7
5.802694235494295e-7
8.524184689962646e-7
3.589463117197928e-6

16-element Array{Float64,1}:
1.0856036488279697e-15
1.1435508643592667e-15
1.3980653454598973e-15
9.052233990134039e-15
1.2267538036369622e-14
3.0090359712477807e-14
3.813726654170083e-14
5.5759518006466014e-14
3.890573567900651e-13
6.284524738551312e-13
1.4822211423114135e-12
4.582949432004653e-12
4.855444113833696e-7
5.82429004776745e-7
8.524184690382938e-7
3.589481273994033e-6