'Convert Matlab script to Python
I am trying to convert a script from the OCamCalib Matlab toolbox to python. This scripts estimate a inverse polynomial function used for the reverse path of light ray on a fisheye lens.
My python scripts fails to compute the polynomial values at the end.
The Matlab script: findinvpoly.m
%FINDINVPOLY finds the inverse polynomial specified in the argument.
% [POL, ERR, N] = FINDINVPOLY(SS, RADIUS, N) finds an approximation of the inverse polynomial specified in OCAM_MODEL.SS.
% The returned polynomial POL is used in WORLD2CAM_FAST to compute the reprojected point very efficiently.
%
% SS is the polynomial which describe the mirrror/lens model.
% RADIUS is the radius (pixels) of the omnidirectional picture.
% ERR is the error (pixel) that you commit in using the returned
% polynomial instead of the inverse SS. N is searched so that
% that ERR is < 0.01 pixels.
%
% Copyright (C) 2008 DAVIDE SCARAMUZZA, ETH Zurich
% Author: Davide Scaramuzza - email: [email protected]
function [pol, err, N] = findinvpoly(ss, radius)
if nargin < 3
maxerr = inf;
N = 1;
while maxerr > 0.01 %Repeat until the reprojection error is smaller than 0.01 pixels
N = N + 1;
[pol, err] = findinvpoly2(ss, radius, N);
maxerr = max(err);
end
else
[pol, err, N] = findinvpoly2(ss, radius, N)
end
function [pol, err, N] = findinvpoly2(ss, radius, N)
theta = [-pi/2:0.01:1.20];
r = invFUN(ss, theta, radius);
ind = find(r~=inf);
theta = theta(ind);
r = r(ind);
pol = polyfit(theta,r,N);
err = abs( r - polyval(pol, theta)); %approximation error in pixels
function r=invFUN(ss, theta, radius)
m=tan(theta);
r=[];
poly_coef=ss(end:-1:1);
poly_coef_tmp=poly_coef;
for j=1:length(m)
poly_coef_tmp(end-1)=poly_coef(end-1)-m(j);
rhoTmp=roots(poly_coef_tmp);
res=rhoTmp(find(imag(rhoTmp)==0 & rhoTmp>0 & rhoTmp<radius ));
if isempty(res) | length(res)>1
r(j)=inf;
else
r(j)=res;
end
end
The script calling the matlab function : main.m
coeff = [-1.729819e+03, 0.000000e+00, 2.218967e-04, -2.655544e-08, 2.169698e-11 ];
h = 2160;
w = 4096;
radius = sqrt(h*h + w*w);
[pol, err, N] = findinvpoly(coeff, radius);
pol
The output I have is :
pol =
1.0e+03 *
0.0002 0.0020 0.0057 0.0002 -0.0163 -0.0015 0.0345 0.0068 -0.0082 0.0463 0.0192 0.0751 0.2371 -0.0350 1.2974 2.4535
My Python script :
import numpy as np
import math
def find_inverse_polynome(poly_coefficients, radius, N=1, threshold_error = 0.01):
max_error = np.inf
while max_error > threshold_error:
N = N+1
poly, error = compute_inverse_polynome(poly_coefficients, radius, N)
max_error = error.max()
return poly, error, N
def compute_inverse_polynome(poly_coefficients, radius, N):
theta = np.arange(-np.pi, 1.2, 0.01)
r = np.array(inverse_poly(poly_coefficients, theta, radius))
indices = np.where(r != np.inf)[0]
theta = theta[indices]
r = r[indices]
poly = np.polyfit(theta, r, N)
print(theta)
exit()
error = abs(r - np.polyval(poly, theta)) # approximation of error in pixels
return poly, error
def inverse_poly(poly_coefficients, theta, radius):
m = np.tan(theta)
result = []
reversed_poly_coefficients = poly_coefficients[::-1]
reversed_poly_coefficients_temp = reversed_poly_coefficients
for idx, j in enumerate(range(len(m))):
reversed_poly_coefficients_temp[-2] = reversed_poly_coefficients[-2] - m[j]
rho_temp = np.roots(reversed_poly_coefficients_temp)
res = rho_temp[np.where((np.imag(rho_temp) == 0) & (rho_temp < radius) & (rho_temp>0))]
if res.size == 0 or res.size > 1:
result.append(np.inf)
else:
result.append(res)
return result
poly_coeff = [-1.729819e+03, 0.000000e+00, 2.218967e-04, -2.655544e-08, 2.169698e-11 ]
poly_coeff = [1707.2, 0, -0.00023066, 2.4174E-08, -1.8933E-11]
h = 2160
w = 4096
radius = math.sqrt(h*h + w*w)
a = find_inverse_polynome(poly_coeff, radius)
The current error I am having is :
Traceback (most recent call last):
File "/media/alix/System/Users/aleroy/Documents/DPhil/algorithms/opticalflow/generate_invpoly.py", line 52, in <module>
a = find_inverse_polynome(poly_coeff, radius)
File "/media/alix/System/Users/aleroy/Documents/DPhil/algorithms/opticalflow/generate_invpoly.py", line 9, in find_inverse_polynome
poly, error = compute_inverse_polynome(poly_coefficients, radius, N)
File "/media/alix/System/Users/aleroy/Documents/DPhil/algorithms/opticalflow/generate_invpoly.py", line 22, in compute_inverse_polynome
poly = np.polyfit(theta, r, N)
File "<__array_function__ internals>", line 5, in polyfit
File "/home/alix/Documents/dphil/venv/p38/lib/python3.8/site-packages/numpy/lib/polynomial.py", line 629, in polyfit
c, resids, rank, s = lstsq(lhs, rhs, rcond)
File "<__array_function__ internals>", line 5, in lstsq
File "/home/alix/Documents/dphil/venv/p38/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 2306, in lstsq
x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
TypeError: No loop matching the specified signature and casting was found for ufunc lstsq_n
EDIT : I was able to modify the code and it now outputs values. However the system does not converge. One condition was missing in the np.where of the inverse_poly() function. Now, I have this warning message :
RankWarning: Polyfit may be poorly conditioned
poly, error = compute_inverse_polynome(poly_coefficients, radius, N)
As requested in the comments, this is the current value of r:
[[2.44779367e+03+0.j]
[2.43423568e+03+0.j]
[2.40712477e+03+0.j]
[2.36649200e+03+0.j]
[2.31243826e+03+0.j]
[2.24520354e+03+0.j]
[2.16526289e+03+0.j]
[2.07344374e+03+0.j]
[1.97104699e+03+0.j]
[1.85993558e+03+0.j]
[1.74253978e+03+0.j]
[1.62173327e+03+0.j]
[1.50057400e+03+0.j]
[1.38196726e+03+0.j]
[1.26835405e+03+0.j]
[1.16151748e+03+0.j]
[1.06254023e+03+0.j]
[9.71883778e+02+0.j]
[8.89530856e+02+0.j]
[8.15138353e+02+0.j]
[7.48169221e+02+0.j]
[6.87991946e+02+0.j]
[6.33948462e+02+0.j]
[5.85396629e+02+0.j]
[5.41734273e+02+0.j]
[5.02410876e+02+0.j]
[4.66931443e+02+0.j]
[4.34855756e+02+0.j]
[4.05795080e+02+0.j]
[3.79407661e+02+0.j]
[3.55393816e+02+0.j]
[3.33491061e+02+0.j]
[3.13469549e+02+0.j]
[2.95127906e+02+0.j]
[2.78289525e+02+0.j]
[2.62799306e+02+0.j]
[2.48520811e+02+0.j]
[2.35333796e+02+0.j]
[2.23132080e+02+0.j]
[2.11821707e+02+0.j]
[2.01319365e+02+0.j]
[1.91551022e+02+0.j]
[1.82450751e+02+0.j]
[1.73959723e+02+0.j]
[1.66025327e+02+0.j]
[1.58600421e+02+0.j]
[1.51642675e+02+0.j]
[1.45114005e+02+0.j]
[1.38980084e+02+0.j]
[1.33209908e+02+0.j]
[1.27775427e+02+0.j]
[1.22651214e+02+0.j]
[1.17814185e+02+0.j]
[1.13243344e+02+0.j]
[1.08919561e+02+0.j]
[1.04825384e+02+0.j]
[1.00944862e+02+0.j]
[9.72633960e+01+0.j]
[9.37676047e+01+0.j]
[9.04452052e+01+0.j]
[8.72849075e+01+0.j]
[8.42763201e+01+0.j]
[8.14098661e+01+0.j]
[7.86767082e+01+0.j]
[7.60686813e+01+0.j]
[7.35782325e+01+0.j]
[7.11983671e+01+0.j]
[6.89225998e+01+0.j]
[6.67449111e+01+0.j]
[6.46597079e+01+0.j]
[6.26617875e+01+0.j]
[6.07463058e+01+0.j]
[5.89087481e+01+0.j]
[5.71449021e+01+0.j]
[5.54508347e+01+0.j]
[5.38228693e+01+0.j]
[5.22575668e+01+0.j]
[5.07517067e+01+0.j]
[4.93022710e+01+0.j]
[4.79064291e+01+0.j]
[4.65615237e+01+0.j]
[4.52650585e+01+0.j]
[4.40146865e+01+0.j]
[4.28081993e+01+0.j]
[4.16435173e+01+0.j]
[4.05186810e+01+0.j]
[3.94318424e+01+0.j]
[3.83812575e+01+0.j]
[3.73652794e+01+0.j]
[3.63823519e+01+0.j]
[3.54310030e+01+0.j]
[3.45098401e+01+0.j]
[3.36175442e+01+0.j]
[3.27528658e+01+0.j]
[3.19146198e+01+0.j]
[3.11016821e+01+0.j]
[3.03129852e+01+0.j]
[2.95475152e+01+0.j]
[2.88043081e+01+0.j]
[2.80824470e+01+0.j]
[2.73810591e+01+0.j]
[2.66993132e+01+0.j]
[2.60364170e+01+0.j]
[2.53916149e+01+0.j]
[2.47641859e+01+0.j]
[2.41534412e+01+0.j]
[2.35587228e+01+0.j]
[2.29794012e+01+0.j]
[2.24148742e+01+0.j]
[2.18645646e+01+0.j]
[2.13279196e+01+0.j]
[2.08044085e+01+0.j]
[2.02935217e+01+0.j]
[1.97947697e+01+0.j]
[1.93076813e+01+0.j]
[1.88318027e+01+0.j]
[1.83666965e+01+0.j]
[1.79119404e+01+0.j]
[1.74671260e+01+0.j]
[1.70318581e+01+0.j]
[1.66057535e+01+0.j]
[1.61884399e+01+0.j]
[1.57795551e+01+0.j]
[1.53787458e+01+0.j]
[1.49856666e+01+0.j]
[1.45999792e+01+0.j]
[1.42213508e+01+0.j]
[1.38494533e+01+0.j]
[1.34839622e+01+0.j]
[1.31245549e+01+0.j]
[1.27709093e+01+0.j]
[1.24227025e+01+0.j]
[1.20796086e+01+0.j]
[1.17412968e+01+0.j]
[1.14074293e+01+0.j]
[1.10776577e+01+0.j]
[1.07516204e+01+0.j]
[1.04289382e+01+0.j]
[1.01092096e+01+0.j]
[9.79200415e+00+0.j]
[9.47685546e+00+0.j]
[9.16325079e+00+0.j]
[8.85061852e+00+0.j]
[8.53831107e+00+0.j]
[8.22558216e+00+0.j]
[7.91155518e+00+0.j]
[7.59517858e+00+0.j]
[7.27516091e+00+0.j]
[6.94987288e+00+0.j]
[6.61719378e+00+0.j]
[6.27425842e+00+0.j]
[5.91701393e+00+0.j]
[5.53937971e+00+0.j]
[5.13147807e+00+0.j]
[4.67530431e+00+0.j]
[4.13134483e+00+0.j]
[3.37490430e+00+0.j]
[9.69109951e-01+0.j]
[1.03280934e+00+0.j]
[1.06639997e+00+0.j]
[1.08970138e+00+0.j]
[1.10772769e+00+0.j]
[1.12251856e+00+0.j]
[1.13511044e+00+0.j]
[1.14610442e+00+0.j]
[1.15588125e+00+0.j]
[1.16469793e+00+0.j]
[1.17273640e+00+0.j]
[1.18013037e+00+0.j]
[1.18698108e+00+0.j]
[1.19336713e+00+0.j]
[1.19935080e+00+0.j]
[1.20498229e+00+0.j]
[1.21030271e+00+0.j]
[1.21534615e+00+0.j]
[1.22014122e+00+0.j]
[1.22471214e+00+0.j]
[1.22907965e+00+0.j]
[1.23326161e+00+0.j]
[1.23727355e+00+0.j]
[1.24112901e+00+0.j]
[1.24483994e+00+0.j]
[1.24841688e+00+0.j]
[1.25186920e+00+0.j]
[1.25520528e+00+0.j]
[1.25843261e+00+0.j]
[1.26155796e+00+0.j]
[1.26458744e+00+0.j]
[1.26752657e+00+0.j]
[1.27038039e+00+0.j]
[1.27315350e+00+0.j]
[1.27585010e+00+0.j]
[1.27847405e+00+0.j]
[1.28102889e+00+0.j]
[1.28351789e+00+0.j]
[1.28594406e+00+0.j]
[1.28831021e+00+0.j]
[1.29061893e+00+0.j]
[1.29287263e+00+0.j]
[1.29507354e+00+0.j]
[1.29722377e+00+0.j]
[1.29932526e+00+0.j]
[1.30137985e+00+0.j]
[1.30338924e+00+0.j]
[1.30535505e+00+0.j]
[1.30727878e+00+0.j]
[1.30916185e+00+0.j]
[1.31100560e+00+0.j]
[1.31281129e+00+0.j]
[1.31458012e+00+0.j]
[1.31631320e+00+0.j]
[1.31801160e+00+0.j]
[1.31967633e+00+0.j]
[1.32130834e+00+0.j]
[1.32290854e+00+0.j]
[1.32447779e+00+0.j]
[1.32601691e+00+0.j]
[1.32752666e+00+0.j]
[1.32900780e+00+0.j]
[1.33046103e+00+0.j]
[1.33188702e+00+0.j]
[1.33328641e+00+0.j]
[1.33465981e+00+0.j]
[1.33600780e+00+0.j]
[1.33733096e+00+0.j]
[1.33862980e+00+0.j]
[1.33990484e+00+0.j]
[1.34115656e+00+0.j]
[1.34238545e+00+0.j]
[1.34359194e+00+0.j]
[1.34477647e+00+0.j]
[1.34593945e+00+0.j]
[1.34708128e+00+0.j]
[1.34820234e+00+0.j]
[1.34930298e+00+0.j]
[1.35038358e+00+0.j]
[1.35144446e+00+0.j]
[1.35248594e+00+0.j]
[1.35350835e+00+0.j]
[1.35451198e+00+0.j]
[1.35549712e+00+0.j]
[1.35646406e+00+0.j]
[1.35741305e+00+0.j]
[1.35834437e+00+0.j]
[1.35925825e+00+0.j]
[1.36015495e+00+0.j]
[1.36103469e+00+0.j]
[1.36189770e+00+0.j]
[1.36274419e+00+0.j]
[1.36357437e+00+0.j]
[1.36438845e+00+0.j]
[1.36518661e+00+0.j]
[1.36596906e+00+0.j]
[1.36673596e+00+0.j]
[1.36748750e+00+0.j]
[1.36822384e+00+0.j]
[1.36894514e+00+0.j]
[1.36965157e+00+0.j]
[1.37034328e+00+0.j]
[1.37102042e+00+0.j]
[1.37168312e+00+0.j]
[1.37233153e+00+0.j]
[1.37296577e+00+0.j]
[1.37358599e+00+0.j]
[1.37419230e+00+0.j]
[1.37478482e+00+0.j]
[1.37536368e+00+0.j]
[1.37592897e+00+0.j]
[1.37648081e+00+0.j]
[1.37701932e+00+0.j]
[1.37754457e+00+0.j]
[1.37805668e+00+0.j]
[1.37855575e+00+0.j]
[1.37904185e+00+0.j]
[1.37951508e+00+0.j]
[1.37997552e+00+0.j]
[1.38042326e+00+0.j]
[1.38085838e+00+0.j]
[1.38128094e+00+0.j]
[1.38169103e+00+0.j]
[1.38208872e+00+0.j]
[1.38247406e+00+0.j]
[1.38284714e+00+0.j]
[1.38320801e+00+0.j]
[1.38355672e+00+0.j]
[1.38389335e+00+0.j]
[1.38421795e+00+0.j]
[1.38453056e+00+0.j]
[1.38483124e+00+0.j]
[1.38512004e+00+0.j]
[1.38539701e+00+0.j]
[1.38566219e+00+0.j]
[1.38591561e+00+0.j]
[1.38615734e+00+0.j]
[1.38638739e+00+0.j]
[1.38660581e+00+0.j]
[1.38681263e+00+0.j]
[1.38700788e+00+0.j]
[1.38719160e+00+0.j]
[1.38736382e+00+0.j]
[1.38752455e+00+0.j]
[1.38767383e+00+0.j]
[1.38781167e+00+0.j]
[1.38793810e+00+0.j]
[1.38805314e+00+0.j]
[1.38815681e+00+0.j]
[1.38824912e+00+0.j]
[1.38833008e+00+0.j]
[1.38839971e+00+0.j]
[1.38845802e+00+0.j]
[1.38850502e+00+0.j]
[1.38854071e+00+0.j]
[1.38856510e+00+0.j]
[1.38857819e+00+0.j]
[1.38857999e+00+0.j]
[1.38857049e+00+0.j]
[1.38854970e+00+0.j]
[1.38851761e+00+0.j]
[1.38847422e+00+0.j]
[1.38841951e+00+0.j]
[1.38835349e+00+0.j]
[1.38827614e+00+0.j]
[1.38818745e+00+0.j]
[1.38808740e+00+0.j]
[1.38797598e+00+0.j]
[1.38785318e+00+0.j]
[1.38771897e+00+0.j]
[1.38757334e+00+0.j]
[1.38741626e+00+0.j]
[1.38724771e+00+0.j]
[1.38706766e+00+0.j]
[1.38687608e+00+0.j]
[1.38667294e+00+0.j]
[1.38645822e+00+0.j]
[1.38623188e+00+0.j]
[1.38599388e+00+0.j]
[1.38574418e+00+0.j]
[1.38548275e+00+0.j]
[1.38520955e+00+0.j]
[1.38492452e+00+0.j]
[1.38462763e+00+0.j]
[1.38431882e+00+0.j]
[1.38399805e+00+0.j]
[1.38366526e+00+0.j]
[1.38332040e+00+0.j]
[1.38296341e+00+0.j]
[1.38259423e+00+0.j]
[1.38221280e+00+0.j]
[1.38181905e+00+0.j]
[1.38141292e+00+0.j]
[1.38099434e+00+0.j]
[1.38056323e+00+0.j]
[1.38011952e+00+0.j]
[1.37966313e+00+0.j]
[1.37919398e+00+0.j]
[1.37871199e+00+0.j]
[1.37821706e+00+0.j]
[1.37770912e+00+0.j]
[1.37718806e+00+0.j]
[1.37665379e+00+0.j]
[1.37610620e+00+0.j]
[1.37554521e+00+0.j]
[1.37497068e+00+0.j]
[1.37438253e+00+0.j]
[1.37378062e+00+0.j]
[1.37316485e+00+0.j]
[1.37253508e+00+0.j]
[1.37189120e+00+0.j]
[1.37123307e+00+0.j]
[1.37056054e+00+0.j]
[1.36987349e+00+0.j]
[1.36917177e+00+0.j]
[1.36845522e+00+0.j]
[1.36772368e+00+0.j]
[1.36697701e+00+0.j]
[1.36621502e+00+0.j]
[1.36543754e+00+0.j]
[1.36464441e+00+0.j]
[1.36383542e+00+0.j]
[1.36301039e+00+0.j]
[1.36216911e+00+0.j]
[1.36131139e+00+0.j]
[1.36043700e+00+0.j]
[1.35954573e+00+0.j]
[1.35863735e+00+0.j]
[1.35771161e+00+0.j]
[1.35676828e+00+0.j]
[1.35580709e+00+0.j]
[1.35482777e+00+0.j]
[1.35383006e+00+0.j]
[1.35281367e+00+0.j]
[1.35177829e+00+0.j]
[1.35072362e+00+0.j]
[1.34964935e+00+0.j]
[1.34855512e+00+0.j]
[1.34744061e+00+0.j]
[1.34630544e+00+0.j]
[1.34514924e+00+0.j]
[1.34397162e+00+0.j]
[1.34277216e+00+0.j]
[1.34155046e+00+0.j]
[1.34030606e+00+0.j]
[1.33903849e+00+0.j]
[1.33774729e+00+0.j]
[1.33643194e+00+0.j]
[1.33509192e+00+0.j]
[1.33372667e+00+0.j]
[1.33233562e+00+0.j]
[1.33091818e+00+0.j]
[1.32947370e+00+0.j]
[1.32800152e+00+0.j]
[1.32650096e+00+0.j]
[1.32497128e+00+0.j]
[1.32341172e+00+0.j]
[1.32182147e+00+0.j]
[1.32019968e+00+0.j]
[1.31854548e+00+0.j]
[1.31685791e+00+0.j]
[1.31513599e+00+0.j]
[1.31337867e+00+0.j]
[1.31158484e+00+0.j]
[1.30975335e+00+0.j]
[1.30788294e+00+0.j]
[1.30597231e+00+0.j]
[1.30402007e+00+0.j]
[1.30202473e+00+0.j]
[1.29998472e+00+0.j]
[1.29789836e+00+0.j]
[1.29576386e+00+0.j]
[1.29357931e+00+0.j]
[1.29134267e+00+0.j]
[1.28905173e+00+0.j]
[1.28670414e+00+0.j]
[1.28429738e+00+0.j]
[1.28182871e+00+0.j]]
This is the current value of theta:
[-3.14159265e+00 -3.13159265e+00 -3.12159265e+00 -3.11159265e+00
-3.10159265e+00 -3.09159265e+00 -3.08159265e+00 -3.07159265e+00
-3.06159265e+00 -3.05159265e+00 -3.04159265e+00 -3.03159265e+00
-3.02159265e+00 -3.01159265e+00 -3.00159265e+00 -2.99159265e+00
-2.98159265e+00 -2.97159265e+00 -2.96159265e+00 -2.95159265e+00
-2.94159265e+00 -2.93159265e+00 -2.92159265e+00 -2.91159265e+00
-2.90159265e+00 -2.89159265e+00 -2.88159265e+00 -2.87159265e+00
-2.86159265e+00 -2.85159265e+00 -2.84159265e+00 -2.83159265e+00
-2.82159265e+00 -2.81159265e+00 -2.80159265e+00 -2.79159265e+00
-2.78159265e+00 -2.77159265e+00 -2.76159265e+00 -2.75159265e+00
-2.74159265e+00 -2.73159265e+00 -2.72159265e+00 -2.71159265e+00
-2.70159265e+00 -2.69159265e+00 -2.68159265e+00 -2.67159265e+00
-2.66159265e+00 -2.65159265e+00 -2.64159265e+00 -2.63159265e+00
-2.62159265e+00 -2.61159265e+00 -2.60159265e+00 -2.59159265e+00
-2.58159265e+00 -2.57159265e+00 -2.56159265e+00 -2.55159265e+00
-2.54159265e+00 -2.53159265e+00 -2.52159265e+00 -2.51159265e+00
-2.50159265e+00 -2.49159265e+00 -2.48159265e+00 -2.47159265e+00
-2.46159265e+00 -2.45159265e+00 -2.44159265e+00 -2.43159265e+00
-2.42159265e+00 -2.41159265e+00 -2.40159265e+00 -2.39159265e+00
-2.38159265e+00 -2.37159265e+00 -2.36159265e+00 -2.35159265e+00
-2.34159265e+00 -2.33159265e+00 -2.32159265e+00 -2.31159265e+00
-2.30159265e+00 -2.29159265e+00 -2.28159265e+00 -2.27159265e+00
-2.26159265e+00 -2.25159265e+00 -2.24159265e+00 -2.23159265e+00
-2.22159265e+00 -2.21159265e+00 -2.20159265e+00 -2.19159265e+00
-2.18159265e+00 -2.17159265e+00 -2.16159265e+00 -2.15159265e+00
-2.14159265e+00 -2.13159265e+00 -2.12159265e+00 -2.11159265e+00
-2.10159265e+00 -2.09159265e+00 -2.08159265e+00 -2.07159265e+00
-2.06159265e+00 -2.05159265e+00 -2.04159265e+00 -2.03159265e+00
-2.02159265e+00 -2.01159265e+00 -2.00159265e+00 -1.99159265e+00
-1.98159265e+00 -1.97159265e+00 -1.96159265e+00 -1.95159265e+00
-1.94159265e+00 -1.93159265e+00 -1.92159265e+00 -1.91159265e+00
-1.90159265e+00 -1.89159265e+00 -1.88159265e+00 -1.87159265e+00
-1.86159265e+00 -1.85159265e+00 -1.84159265e+00 -1.83159265e+00
-1.82159265e+00 -1.81159265e+00 -1.80159265e+00 -1.79159265e+00
-1.78159265e+00 -1.77159265e+00 -1.76159265e+00 -1.75159265e+00
-1.74159265e+00 -1.73159265e+00 -1.72159265e+00 -1.71159265e+00
-1.70159265e+00 -1.69159265e+00 -1.68159265e+00 -1.67159265e+00
-1.66159265e+00 -1.65159265e+00 -1.64159265e+00 -1.63159265e+00
-1.62159265e+00 -1.61159265e+00 -1.60159265e+00 -1.59159265e+00
-1.58159265e+00 -1.57159265e+00 -1.56159265e+00 -1.55159265e+00
-1.54159265e+00 -1.53159265e+00 -1.52159265e+00 -1.51159265e+00
-1.50159265e+00 -1.49159265e+00 -1.48159265e+00 -1.47159265e+00
-1.46159265e+00 -1.45159265e+00 -1.44159265e+00 -1.43159265e+00
-1.42159265e+00 -1.41159265e+00 -1.40159265e+00 -1.39159265e+00
-1.38159265e+00 -1.37159265e+00 -1.36159265e+00 -1.35159265e+00
-1.34159265e+00 -1.33159265e+00 -1.32159265e+00 -1.31159265e+00
-1.30159265e+00 -1.29159265e+00 -1.28159265e+00 -1.27159265e+00
-1.26159265e+00 -1.25159265e+00 -1.24159265e+00 -1.23159265e+00
-1.22159265e+00 -1.21159265e+00 -1.20159265e+00 -1.19159265e+00
-1.18159265e+00 -1.17159265e+00 -1.16159265e+00 -1.15159265e+00
-1.14159265e+00 -1.13159265e+00 -1.12159265e+00 -1.11159265e+00
-1.10159265e+00 -1.09159265e+00 -1.08159265e+00 -1.07159265e+00
-1.06159265e+00 -1.05159265e+00 -1.04159265e+00 -1.03159265e+00
-1.02159265e+00 -1.01159265e+00 -1.00159265e+00 -9.91592654e-01
-9.81592654e-01 -9.71592654e-01 -9.61592654e-01 -9.51592654e-01
-9.41592654e-01 -9.31592654e-01 -9.21592654e-01 -9.11592654e-01
-9.01592654e-01 -8.91592654e-01 -8.81592654e-01 -8.71592654e-01
-8.61592654e-01 -8.51592654e-01 -8.41592654e-01 -8.31592654e-01
-8.21592654e-01 -8.11592654e-01 -8.01592654e-01 -7.91592654e-01
-7.81592654e-01 -7.71592654e-01 -7.61592654e-01 -7.51592654e-01
-7.41592654e-01 -7.31592654e-01 -7.21592654e-01 -7.11592654e-01
-7.01592654e-01 -6.91592654e-01 -6.81592654e-01 -6.71592654e-01
-6.61592654e-01 -6.51592654e-01 -6.41592654e-01 -6.31592654e-01
-6.21592654e-01 -6.11592654e-01 -6.01592654e-01 -5.91592654e-01
-5.81592654e-01 -5.71592654e-01 -5.61592654e-01 -5.51592654e-01
-5.41592654e-01 -5.31592654e-01 -5.21592654e-01 -5.11592654e-01
-5.01592654e-01 -4.91592654e-01 -4.81592654e-01 -4.71592654e-01
-4.61592654e-01 -4.51592654e-01 -4.41592654e-01 -4.31592654e-01
-4.21592654e-01 -4.11592654e-01 -4.01592654e-01 -3.91592654e-01
-3.81592654e-01 -3.71592654e-01 -3.61592654e-01 -3.51592654e-01
-3.41592654e-01 -3.31592654e-01 -3.21592654e-01 -3.11592654e-01
-3.01592654e-01 -2.91592654e-01 -2.81592654e-01 -2.71592654e-01
-2.61592654e-01 -2.51592654e-01 -2.41592654e-01 -2.31592654e-01
-2.21592654e-01 -2.11592654e-01 -2.01592654e-01 -1.91592654e-01
-1.81592654e-01 -1.71592654e-01 -1.61592654e-01 -1.51592654e-01
-1.41592654e-01 -1.31592654e-01 -1.21592654e-01 -1.11592654e-01
-1.01592654e-01 -9.15926536e-02 -8.15926536e-02 -7.15926536e-02
-6.15926536e-02 -5.15926536e-02 -4.15926536e-02 -3.15926536e-02
-2.15926536e-02 -1.15926536e-02 -1.59265359e-03 8.40734641e-03
1.84073464e-02 2.84073464e-02 3.84073464e-02 4.84073464e-02
5.84073464e-02 6.84073464e-02 7.84073464e-02 8.84073464e-02
9.84073464e-02 1.08407346e-01 1.18407346e-01 1.28407346e-01
1.38407346e-01 1.48407346e-01 1.58407346e-01 1.68407346e-01
1.78407346e-01 1.88407346e-01 1.98407346e-01 2.08407346e-01
2.18407346e-01 2.28407346e-01 2.38407346e-01 2.48407346e-01
2.58407346e-01 2.68407346e-01 2.78407346e-01 2.88407346e-01
2.98407346e-01 3.08407346e-01 3.18407346e-01 3.28407346e-01
3.38407346e-01 3.48407346e-01 3.58407346e-01 3.68407346e-01
3.78407346e-01 3.88407346e-01 3.98407346e-01 4.08407346e-01
4.18407346e-01 4.28407346e-01 4.38407346e-01 4.48407346e-01
4.58407346e-01 4.68407346e-01 4.78407346e-01 4.88407346e-01
4.98407346e-01 5.08407346e-01 5.18407346e-01 5.28407346e-01
5.38407346e-01 5.48407346e-01 5.58407346e-01 5.68407346e-01
5.78407346e-01 5.88407346e-01 5.98407346e-01 6.08407346e-01
6.18407346e-01 6.28407346e-01 6.38407346e-01 6.48407346e-01
6.58407346e-01 6.68407346e-01 6.78407346e-01 6.88407346e-01
6.98407346e-01 7.08407346e-01 7.18407346e-01 7.28407346e-01
7.38407346e-01 7.48407346e-01 7.58407346e-01 7.68407346e-01
7.78407346e-01 7.88407346e-01 7.98407346e-01 8.08407346e-01
8.18407346e-01 8.28407346e-01 8.38407346e-01 8.48407346e-01
8.58407346e-01 8.68407346e-01 8.78407346e-01 8.88407346e-01
8.98407346e-01 9.08407346e-01 9.18407346e-01 9.28407346e-01
9.38407346e-01 9.48407346e-01 9.58407346e-01 9.68407346e-01
9.78407346e-01 9.88407346e-01 9.98407346e-01 1.00840735e+00
1.01840735e+00 1.02840735e+00 1.03840735e+00 1.04840735e+00
1.05840735e+00 1.06840735e+00 1.07840735e+00 1.08840735e+00
1.09840735e+00 1.10840735e+00 1.11840735e+00 1.12840735e+00
1.13840735e+00 1.14840735e+00 1.15840735e+00 1.16840735e+00
1.17840735e+00 1.18840735e+00 1.19840735e+00]
The Matlab script converges in a splut of a second, while my algorithm seems to run for a very long time. The output error is very large since the value of r is too small (see output above).
Solution 1:[1]
I actually happened to port most of the OCamCalib toolbox over to Python at some point, so here's what I ended up with (edited a bit to match your naming):
import cv2
import numpy as np
def find_inverse_polynome(poly_coefficients, radius, N=1, threshold_error=0.01):
max_error = np.inf
while max_error > threshold_error:
N += 1
poly, error = compute_inverse_polynome(poly_coefficients, radius, N)
max_error = error.max()
return poly, error, N
def compute_inverse_polynome(poly_coefficients, radius, N):
theta = np.arange(-np.pi / 2, 1.2, 0.01)
r = inverse_poly(poly_coefficients, theta)
valid = r < radius
theta = theta[valid]
r = r[valid]
poly = np.polyfit(theta, r, N)
error = np.abs(r - np.polyval(poly, theta)) # approximation of error in pixels
return poly, error
def inverse_poly(poly_coefficients, theta):
m = np.tan(theta)
r = np.zeros(len(m))
poly_coef = poly_coefficients[::-1]
poly_coef_tmp = poly_coef.copy()
for i in range(len(m)):
poly_coef_tmp[-2] = poly_coef[-2] - m[i]
rho_tmp = cv2.solvePoly(poly_coef_tmp[::-1], maxIters=50)[1][:, 0]
res = rho_tmp[(rho_tmp[:, 0] > 0) & (np.abs(rho_tmp[:, 1]) < 1e-10), 0]
r[i] = np.min(res) if len(res) > 0 else np.inf
return r
This makes use of cv2.solvePoly() instead of np.roots() due to the OpenCV's version being several times faster for some reason.
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|---|
| Solution 1 | Martin Valgur |