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.

[2]:
include("CoxHomotopy.jl")
[2]:
trackPathRandom (generic function with 1 method)

Next, we generate the equations.

[3]:
@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);

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

[4]:
(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

[10]:
@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.

[11]:
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