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.

using Pkg
Pkg.activate(".")

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

We include our routines.

include("CoxHomotopy.jl")

polymake version 4.2
Copyright (c) 1997-2020
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.

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][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.

(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

@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.

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