[PDF] Stable radial distortion calibration by polynomial matrix inequalities

View - Download Stable radial distortion calibration by polynomial matrix inequalities

Previous PDF

Next PDF

Stable radial distortion calibration by polynomial matrix inequalities programming

Jan Heller

1, Didier Henrion2;3;1, Tomas Pajdla1

Abstract

Polynomial and rational functions are the number one choice when it comes to modeling of radial distortion of lenses. However, several extrapolation and numerical issues may arise while using these functions that have not been covered by the literature much so far. In this paper, we identify these problems and show how to deal with them by enforcing nonnegativity of certain polynomials. Further, we show how to model these nonnegativities using polynomial matrix inequalities (PMI) and how to estimate the radial distortion parameters subject to PMI constraints using semidenite programming (SDP). Finally, we suggest several approaches on how to incorporate the proposed method into the overall camera calibration procedure.

1 Introduction

Radial distortion modeling is the most important non-linear part of the camera calibra- tion process [9]. The rst works on the topic came from the photogrammetric commu- nity [4, 5, 14]. Since then, a plethora of models has been suggested in the literature [16]. Among the proposed models, the ones based on polynomial and rational functions are the most popular. This popularity undoubtedly stems from the fact that these function are easily manipulated and yet provide sucient tting power for wide range or distortions. Unfortunately, the extrapolation qualities of polynomials can be quite unpredictable in situations where little or no data is available. However, even if data points are missing, the overall shape of the distortion is knowna prioriin many calibration scenarios,e.g., the lens introduces barrel or pincushion distortions. Based on sucha prioriinformation, the shape of the polynomial and rational distortion functions can be controlled by enforcing nonnegativity of certain polynomials. For example, in the case of pincushion distortion we can accomplish the desired shape by enforcing nonnegativity of the rst and the second derivatives of the distortion function on the whole eld of view of the camera. In this paper, we propose a radial distortion calibration procedure where a polynomial cost function,e.g., reprojection error, is minimized subject to such shape constraints.1 Faculty of Electrical Engineering, Czech Technical University in Prague, CZ-16627 Praha 6, Tech- nicka 2, Czech Republic.

2CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France.

3Universite de Toulouse; F-31400 Toulouse; France.

1arXiv.1409.5753v1 [math.OC] 19 Sep 2014

2 CAMERA RADIAL DISTORTION2

This shape optimization procedure is designed to stabilize the shape of the distortion function. It is based on polynomial matrix inequalities (PMI) programming and can be easily incorporated into an existing camera calibration procedure. In Section 2, we formally introduce the radial distortion function and present several ex- trapolation issues arising while using polynomial and rational distortion models. Next, in Section 3 we provide a minimal theoretical background needed for our shape stabilization approach. In Section 4, we demonstrate the proposed method on three types of radial distortion shapes and models and show how to incorporate the method into an overall camera calibration procedure. Finally, in Section 5 we experimentally validate our ap- proach and show that the method guarantees the correct shape of a distortion function without compromising the quality of the overall camera calibration as measured by the reprojection error.

2 Camera Radial Distortion

Let us suppose that a set of scene pointsXi2R3,i= 1;:::;nis observed by a camera. IfR2SO(3),t2R3are the camera extrinsic parameters, a scene pointXigets projected into an image point (xi;yi;1)>: i(xi;yi;1)>=RXi+t; i2R: In reality, some amount of radial distortion is always present and the camera observes a point (^xi;^yi;1)>which does not coincide with the ideal (and unobservable) point (xi;yi;1)>. In pixel coordinates, the camera observes a pointK(^xi;^yi;1)>, whereK2R33 is the matrix of intrinsic camera parameters, the so-called calibration matrix.Radial dis- tortion functionL:R!Ris a function of radiusr=px

2i+y2ithat models the radial

displacement of the ideal image point position from the center of the radial distortion as ^xi ^yi =L(r)xi y i :(1) The functionL(r) is only dened forr >0 andL(0) = 1,L(r)>0. For the purposes of demonstration of the proposed shape optimization procedure, we will useL(r) dened as follows L(r) =f(r)g(r)=1 +k1r+k2r2+k3r31 +k4r+k5k2+k6r3;(2) wherek= (k1;k2;:::;k6) is the vector of model parameters. This denition accommo- dates several models already proposed in the literature [12]. However, we will see that the shape optimization procedure holds for any rational function.

2.1 Extrapolation issues of radial distortion calibration

Let us motivate the need for the radial distortion shape optimization by demonstrating two examples of extrapolation issues arising while using polynomial and rational distortion models.

2 CAMERA RADIAL DISTORTION30.1 0.3 0.5 0.7 0.90.60.70.80.911.1

Distance from the center of distortionr

L(r) ??(a)

0.2 0.4 0.6 0.8 1 1.2-0.20.20.611.4

Distance from the center of distortionr

L(r) (b) (d) Figure 1:Calibration issues. Examples of issues arising while using polynomial and rational function for radial distortion calibration. See text for details. First, let's suppose a calibration scenario where images of a calibration target were taken, but the image projections of the known 3D points lie close to the center of the images with no points covering the corners of the images. Figure 1(a) shows in black the graph of the amount of barrel distortion introduced by the camera lens as a function of the distance from the center of the radial distortion. When a polynomial distortion model L(r) =f(r) is used, see Equation 2, in combination with an unconstrained calibration method [21, 2] (in red), the real distortion is tted successfully near the center of the image on intervals where the data points are available (left of the diamond symbol). However, the recovered polynomial quickly drifts away elsewhere (red circles depict the distances of the projections of the image corners). In green, a polynomial recovered by the method proposed in this paper is shown. Here, the negativity of the rst and the second derivatives of the polynomial on the whole eld of view was enforced. This caused the model to t the original distortion much closer on the whole eld of view. Figure 1(b) shows a synthetic checkerboard image (the upper left corner) and the same image distorted by the original barrel distortion (the upper right corner). In the lower left corner, the image is undistorted back using the polynomial recovered by [21, 2]. In the lower right corner, the image successfully undistorted by the polynomial recovered using the proposed shape optimization method is shown. Let us consider a similar calibration scenario to the one from the previous paragraph, this

3 POLYNOMIALS AND PMI PROGRAMMING4

time with a lens causing a mustache type radial distortion, see Figure 1(c). If the radial distortion model is used,L(r) =f(r)g(r), the classical calibration approach [21, 2] is able to correctly recover the original shape. However, the polynomialsf(r) andg(r) share a common root (red dash-dot lines), which causes a numerical instability presented as a sharp spike inL(r) around the common root|an issue we will call thezero-crossing problem. When the nonnegativity ofg(r) is enforced using the proposed approach, not only is the correct shape recovered, but since there is now no root in the eld of view interval (green dash-dot lines), the spike inL(r) is also gone. Figure 1(d) shows a similar arrangement as Figure 1(b), now with only the upper left part of the checkerboard shown. The numerical instability ofL(r) is presented as a notable ringing in the upper left corner of the checkerboard. One can argue that the common root is a consequence of the fact that the degrees off(r);g(r) are higher that needed and that a model with fewer coecients should be used. This may be true in some cases, however, we observed just as many situations where the lower degree polynomials resolved the zero-crossing problem only at the cost of a considerably higher reprojection error.

3 Polynomials and PMI Programming

In this section, we present a minimal theoretical background needed for the proposed shape optimization procedure.

3.1 Polynomials and polynomial matrices

An univariate polynomialp(x)2Rn[x] of degreen2Nis a real function dened as p(x) =pnxn+pn1xn1++p1x+p0=p> n(x); wherep= (p0;p1;:::;pn)>2Rn+1is the vector of coecients with a nonvanishing coecientpnand n(x) = (1;x;x2;:::;xn)>is the canonical basis. Letq(x)2R2n[x]. A symmetric matrixQ2Rn0n0,Q= (qi;j), wheren0=n+ 1, is calledGram matrix associated withq(x) and the basis n(x) [6] if q(x) = > n(x)Q n(x):(3) Generally, there is more than one Gram matrix associated with a polynomialq(x) and we will denote the set of such matrices asG(q(x)).The polynomialq(x) can be expressed in the elements ofQby simply expanding the right hand side of Equation 3 and by comparing the coecients. Letx= (x1;x2;:::;xd)2Rdbe a real vector and= (1;1;:::;d)2Ndan integer vector. Amonomialof degreen=Piis dened asx=Qn i=1xii:A multivariate polynomialp(x)2Rn[x] of degreen2Nis a mapping fromRdtoRdened as a linear combination of monomials up to degreen, p(x) =X jjnp x=X jjnp x11x22xdd= (p)>jjn(x)jjd=p> n(x);

3 POLYNOMIALS AND PMI PROGRAMMING5

wherep2Rmis the vector of coecients and n(x) is the canonical basis ofm=d+n d monomials up to degreen. By a polynomial matrix we will understand a symmetric matrix whose elements are polynomials. In the next,Sn(R[x]) will denote the set ofnn symmetric polynomial matrices. The degree ofP= (pi;j(x))2Sn(R[x]) is the largest degree of all the polynomial elements ofP, degP= maxi;jdegpi;j(x). Besides parameterizing polynomials by the associated Gram matrices, we will also need to \linearize" them,i.e., to substitute every monomialxby a new variabley2R. To do this, we dene theRiesz functional`y:Rn[x]!R[y], a linear functional that for a d-variate polynomial of degreen,p(x) =P px, returns anm-variate polynomial of degree one,`y(p(x)) =P py,m=d+n d. With a slight abuse of notation, we will also use`yas a matrix operator acting onSn(R[x]): ifP2Sn(R[x]), thenP0=`y(P) if and only ifp0i;j(y) =`y(pi;j(x)).

3.2 Polynomials positive on nite intervals

The shape optimization procedure presented in this paper is based on enforcing nonneg- ativity of certain polynomials. Since most of the real cameras have limited elds of view, we only need to control the behavior ofL(r) for valuesr2[0;r], where ris the maximal distance between the center of the radial distortion and an (undistorted) image point. For this, we need to characterize the set of univariate polynomials nonnegative on nite intervals. In [13], based on Markov-Lukacs theorem, Nesterov showed how to characterize such a set using positive semidenite Gram matrices: Theorem 1Let < ,p(x)2R[x]anddegp(x) = 2n. Thenp(x)0for allx2[;] if and only if p(x) =s(x) + (x)(x)t(x); wheres(x) = > n(x)S n(x),t(x) = > n1(x)T n1(x), such thatS;T0(i.e.,S2 G(s(x)),T2 G(t(x))are positive semidenite Gram matrices of polynomialss(x)and t(x), respectively). Ifdegp(x) = 2n+ 1, thenp(x)0for allx2[;]if and only if p(x) = (x)s(x) + (x)t(x); wheres(x) = > n(x)S n(x),t(x) = > n(x)T n(x), such thatS;T0. Even though Theorem 1 is an equivalence, we will only use it as an implication: as long as we will have matricesS;Tthat are positive semidenitive, Theorem 1 guarantees that a polynomialp(x) constructed using these matrices will be nonnegative on a given interval.

3.3 Polynomial Matrix Inequalities

According to Theorem 1, a polynomial is nonnegative on an interval as long the matrices S;Tare positive semidenite. By combining these constraints with a polynomial cost function, we get a problem of polynomial matrix inequalities (PMI) programming. A

PMI program can be formally dened as follows:

3 POLYNOMIALS AND PMI PROGRAMMING6

Problem 1 (Polynomial matrix inequalities program) minimizep(x) subject toGi(x)0; i= 1;:::;m; wherep(x)2R[x];Gi2Sni(R[x]): In general, Problem 1 is a hard non-convex problem. Note however, that if the cost functionp(x) and the matricesGi(x),i= 1;:::;mhave degree one, then Problem 1 reduces to a linear matrix inequality (LMI) program and as such is a semidenite program (SDP) solvable by any available SDP solver. In fact, most of the time the shape optimization problems in this paper lead to such a program. Sometimes still,Gi(x) will not be linear. In such cases, we will use the relaxation approach suggested by Henrion and Lasserre [10]. In [10], the authors proposed a hierarchy of LMI programsP1;P2;:::that produces a monotonically non-decreasing sequence of lower boundsp(x1)p(x2):::on Problem 1 that converges to the global minimump(x). Practically, the series converges top(x) in nitely many steps,i.e., there existsj2N, such thatp(xj) =p(x). The authors also showed how this situation can be detected and how the value ofxcan be extracted from the solution of the relaxation by the tools of linear algebra. Let us show here how to constructP,i.e., the LMI relaxation of Problem 1 of order; see [10] for the technical justication of this procedure. LetG2Sn(R[x]),n=Pm i=1ni denote a block diagonal matrix with matricesGion it's diagonal. Since (8i:Gi(x)0), G(x)0, we can replace the PMI constraintsGi(x)0 with one PMI constraintG(x)0. Next, we construct the so-calledmoment matrixM(y) andlocalizing matrixM(G;y)of

G, dened as

M (y) =`y( (x) > (x)); M (G;y) =`y(( (x) > (x)) G); where denotes the Kronecker product [10]. Let = 1 if degG2, =ddegGe2 otherwise.

Now, we can formally write the relaxationPas

Problem 2 (LMI relaxationPof order)

minimize`y(p(x)) subject toM (G;y)0; M (y)0: As the Riesz functional`ywas used to \linearize" both the cost function and the con- straints, we can easily see that Problem 2 is an LMI program.

4 SHAPE OPTIMIZATION FOR RADIAL DISTORTION CALIBRATION7

4 Shape optimization for radial distortion calibration

In this section, we show how to combine the results presented in Section 3 into the radial distortion shape optimization procedure. Technically, the procedure consists of minimization of a polynomial cost function in the vector of radial distortion parametersk subject to PMI constraints enforcing nonnegativity of certain polynomials in the radiusr. Such a minimization problem is a PMI program that can be dealt with using the approach from Section 3.3. As mentioned in Section 3.2, we only need to control the shape ofL(r) on the interval [0;r]. Note, that ris the maximal distance between the center of the radial distortion and undistortedimage points,i.e., the value of ris not known prior to the actual calibration. The value of ris therefore a user supplied parameter. Fortunately, the proposed method is not very sensitive to the value of this parameter and even a gross overestimate yields minima identical to the ground truth value.

4.1 Unconstrained radial distortion calibration

There are several ways how to determine the vector of parameterskof the distortion functionL(r) [9, 17]. All we need for our shape optimization approach is a polynomial cost function. Here, we will dene and use one of such possible cost functions. Let us rewrite Equation 1 usingL(r) from Equation 2 as g(r)^xi ^yi f(r)xi y i = g(r) ^xif(r)xi g(r) ^yif(r)yi! =0: By factoring out the vector of parameterskand by denoting A i= rxir2xir3xi^xir^xir2^xir3 ryir2yir3yi^yir^yir2^yir3! ;bi=xi^xi y i^yi we get a linear systemAik=bi. Now, we can stackA= (A>1;A>2;:::;A>n)>,b= (b>1;b>2;:::b>n)>and estimate the radial distortion parametersk= (k1;k2;:::;k6) as a solution to an overdetermined systemAk=bin the least square sense,i.e., by min- imizingkAkbk2. Note that for polynomial model,i.e.,g(x) = 1, this corresponds to the minimization of the reprojection error. Let us now express the minimization ofkAkbk2as an LMI program. By expanding kAikbik2= (Aikbi)>(Aikbi) =k>A>iAik2b>iAik+b>ibi and by denotingM=Pn i=1A>iAi,m=2Pn i=1A>ibi,c=Pn i=1b>ibi, we can write the polynomial form of the cost function as kAkbk2=k>Mk+m>k+c:(4) As expected, Equation 4 is a quadratic polynomial inkand by constructionM0,i.e., Mis a positive semidenite matrix. Even though the cost function is quadratic, it can be

4 SHAPE OPTIMIZATION FOR RADIAL DISTORTION CALIBRATION8

converted into a linear function using the Schur complement trick [3]: F= I Lk k >L>m>kc+

0,k>L>Lk+m>k+c

0: By decomposingMasM=L>L,e.g., using the Cholesky or the spectral decomposition [7] (recall thatM0), we can rewrite the minimization of Equation 4 as the following LMI program: Problem 3 (Unconstrained radial distortion calibration) minimize subject toF=I Lk k >L>m>kc+ 0:

4.2 Barrel distortion and the polynomial model

As we can see from the example of barrel radial distortion in Figure 1(a), this type of distortion can be characterized by the negativity of the rst and the second derivatives:

8r2[0;r]:L0(r)0&L00(r)0;(5)

where [0;r] spans the eld of view of the camera. If we consider the polynomial model L(r) =f(r), the constraints above mean that we need to enforce nonnegativity of poly- nomials f0(r) =k12k2r3k3r2;f00(r) =2k26k3r on the interval [0;r]. According to Theorem 1,f0(r)0 for8r2[0;r] i f0(r) =k12k2r3k3r2= 1(r)>S1 1(r) +r(rr)T1;(6) where S

1=s11s12

s 12s13

0;T1= (t11)0:

By expanding the right hand side of Equation 6 and by comparing the polynomial coe- cients, we get a parameterization ofkin the elements ofS1andT1: k1=s11

2k2= 2s12+ rt11

3k3=s13t119

)k= (s11;s1212 rt11;13 (t11s13);0;0;0):(7) Let's apply Theorem 1 tof00(r) to get the following constraint: f00(r) =2k26k3r=rS2+ (rr)T2;S2= (s21)0;T2= (t21)0:(8)

4 SHAPE OPTIMIZATION FOR RADIAL DISTORTION CALIBRATION9

By combining Equations 8 and 7, we can express the entries ofS2andT2in the entries of S 1;T1:

2k2= rt21

6k3=s21t21)

)(s

21=1r(2s12+ 2rs13rt11)

t

21=2r(s12+12

rt11)(9) Now, we have four PMI constraints on the shape ofL(r). If we combine these constraints along with the parameterization ofkfrom Equation 7 with Problem 3, we get a radial distortion calibration problem that enforces a barrel type distortion shape of the resulting distortion model:

Problem 4 (Barrel distortion calibration)

minimize subject toF0;S10;T1= (t11)0; S

2=1r(2s12+ 2rs13rt11)0;

T

2=2r(s12+12

rt11)0:

Problem 4 is a PMI program in 5 variables

;s11;s12;s13;t11. Since both the cost function and the PMI constraints have degree one, Problem 4 is in fact an SDP problem. Once it is solved, the unknown distortion parameterskcan be easily recovered using Equation 7.

4.3 Pincushion distortion and the division model

Let us make an analogous analysis for the pincushion distortion shape and the division modelL(r) =1g(r). This type of distortion is characterized by the nonnegativity of the rst and the second derivatives ofL(r) on the eld of view of the camera [0;r]. From the rst derivative we get the following constraint on the polynomial denominatorg(r): L

0(r) =g0(r)g

2(r))L0(r)0, g0(r)0:

The second derivative yields a bit more complicated constraint: L

00(r) =g(r)h(r)g

4(r)=h(r)g

3(r))L00(r)0,((g(r)0&h(r)0)_

(g(r)0&h(r)0); whereh(r) = 2(g0(r))2g(r)g00(r). However, since we know thatL(r)>0 by denition, we only need to consider the constraintsg(r)0;h(r)0. Let us start with the constraintg(r)0. According to Theorem 1,g(r)0 for8r2[0;r] i g(r) = 1 +k4r+k5r2+k6r3= 1(r)>S1 1(r) + (rr) 1(r)>T1 1(r);(10)

4 SHAPE OPTIMIZATION FOR RADIAL DISTORTION CALIBRATION10

where S

1=s11s12

s 12s13

0;T1=t11t12

t 12t13 0: This leads to the following parameterization ofkas well as to a constraint on the variable t 11:

1 = rt11

k

4=s11t11+ 2rt12

k

5= 2s122t12+ rt13

k

6=s13t139

>;)8 :k= (0;0;0;s11t11+ 2rt12;

2s122t12+ rt13;s13t13)

t

11=1r(11)

By applying Theorem 1 to the constraintg0(r)0, we get g0(r) =k42k5r3k6r2= 1(r)>S2 1(r) +r(rr)T2;(12) where S

2=s21s22

s 22s23
quotesdbs_dbs14.pdfusesText_20

[PDF] opencv camera

[PDF] opencv camera calibration

[PDF] opencv camera calibration c

[PDF] opencv camera calibration distortion model

[PDF] opencv distortion model

[PDF] opencv radial and tangential distortion

[PDF] opencv python tutorials documentation pdf

[PDF] opening business account

[PDF] openldap 2.4 setup

[PDF] openldap administrator's guide

[PDF] openldap create database

[PDF] openldap lib

[PDF] openldap mdb

[PDF] openldap sdk

[PDF] operant conditioning