Numerical Solution for Radial Distortion Rectification in Optical Systems

Abstract

Accurate homography estimation is a crucial step for many computer vision applications. Nevertheless, nonlinear optical camera imaging effects can introduce radial distortion, making unfeasible the pinhole model for homography estimation. In this paper, an algorithm to rectify radially distorted images using the Maple software is proposed. First, the effect of radial distortion is modeled and analyzed. Next, an inverse distortion model is developed. The proposed algorithm allows us to estimate both the homography matrix and the distortion parameters by processing images of a calibration target using the Gauss-Newton approach. Successful estimation of homographies and distortion parameters to correct real-world images is reported.

This work was supported by the Consejo Nacional de Ciencia y Tecnología (CONACYT) by the projects Cátedras-880 and A1-S-28112. Authors thank the support of Instituto Politécnico Nacional by the project SIP-20210845.

Appendices

A Homogeneous coordinates

In this section we introduce the homogeneous coordinates and the direct linear transform (DLT) method for homography estimation. We use the homogeneous coordinates operator $\mathcal {H}[.]$ as in [12] to represent an n-dimensional point in homogeneous coordinates. To clarify, let a Cartesian’s coordinate point $\boldsymbol{u} \in \mathbb R^n$ specified as

$$\begin{aligned} \boldsymbol{u} = \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{bmatrix}, \end{aligned}$$

(18)

whose homogeneous representation is given by

$$\begin{aligned} \boldsymbol{v} = \mathcal {H}_s[\boldsymbol{u}] = \begin{bmatrix} v_1 \\ v_2 \\ \vdots \\ v_{n+1}\end{bmatrix} = \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ s\end{bmatrix}, \end{aligned}$$

(19)

such that $v \in \mathbb R^{n+1}$. Inversely to obtain the Cartesian’s coordinates from homogeneous coordinates representation we use the inverse homogeneous operator defined as

$$\begin{aligned} \boldsymbol{u} = \mathcal {H}_s^{-1}[\boldsymbol{v}] = \frac{s}{\mathcal {S}[\boldsymbol{v}]}\mathcal {H}_0^{-1}[\boldsymbol{v}], \end{aligned}$$

(20)

where $\mathcal {S}$ is the scale operator which returns the last element of $\boldsymbol{v}$ and the null inverse homogeneous operator specified as $\mathcal {H}_0^{-1}[\boldsymbol{v}] = \left[ v_1, v_2, \dots , v_n \right] ^T$.

1.1 A.1 Correspondence Between Planes

Homography estimation is critical in many vision-based applications such as camera calibration [16], perspective correction [12], among others [8]. Projective transformations are linear operations that relate points in homogeneous coordinate space. Homographies are a special kind of projective transformations that define direct and inverse mapping relationships between planes in homogeneous coordinates such as points, lines, and other objects. In this work, we use the point paradigm for correspondences between two planes $\left( \varPi _\mu , \varPi _\rho \right)$ which are denoted as $\mu _i \leftrightarrow \rho _i$ with $\mu _i \in \varPi _\mu$ and $\rho _i \in \varPi _\rho$. Equally important, consider direct and inverse mappings

$$\begin{aligned} \rho _i = \mathcal {H}^{-1}\left[ G\mathcal {H}[\mu _i] \right] , \end{aligned}$$

(21)

and

$$\begin{aligned} \mu _i = \mathcal {H}^{-1}\left[ G^{-1}\mathcal {H}[\rho _i] \right] , \end{aligned}$$

(22)

we may state the problem to estimate the homography that satisfies the equation for the direct mapping and inverse mapping given a set of correspondence points. To estimate a homography given a pair correspondences $\mu _i \leftrightarrow \rho _i$ we may minimize the algebraic distance

$$\begin{aligned} \begin{aligned} \mathcal {H}[\mu _i] \times \mathcal {H}[\rho _i]&= \boldsymbol{0}_3, \end{aligned} \end{aligned}$$

(23)

From Eq. (23), using the matrix form of cross product $\left[ \mathcal {H}[\mu _i] \right] _{\times } \mathcal {H}[\rho _i] = \boldsymbol{0}_3$ we may write

$$\begin{aligned} \begin{bmatrix} 0 &{} -1 &{} \mu _{y,i} \\ 1 &{} 0 &{} -\mu _{x,i} \\ -\mu _{y,i} &{} \mu _{x,i} &{} 0 \\ \end{bmatrix} \begin{bmatrix} \bar{\boldsymbol{g}_1}^T\mathcal {H}[\rho _i]\\ \bar{\boldsymbol{g}_2}^T\mathcal {H}[\rho _i]\\ \bar{\boldsymbol{g}_3}^T\mathcal {H}[\rho _i]\\ \end{bmatrix} = \boldsymbol{0}_3, \end{aligned}$$

(24)

where $\bar{\boldsymbol{g}_i}^T$ are the matrix rows of G, and rearranging terms from Eq. (24) we obtain the measurements matrix M and the homography vector $\boldsymbol{g}$ specified as follows

$$\begin{aligned} \underbrace{ \begin{bmatrix} \boldsymbol{0}^T &{} -\mathcal {H}[\rho _i]^T &{} \mu _{y,i}\mathcal {H}[\rho _i]^T\\ \mathcal {H}[\rho _i]^T &{} \boldsymbol{0}^T &{} -\mu _{x,i}\mathcal {H}[\rho _i]^T\\ -\mu _{y,i}\mathcal {H}[\rho _i]^T &{} \mu _{x,i}\mathcal {H}[\rho _i]^T &{} \boldsymbol{0}^T \end{bmatrix}}_{M} \underbrace{ \begin{bmatrix} \bar{\boldsymbol{g}_1}^T \\ \bar{\boldsymbol{g}_2}^T \\ \bar{\boldsymbol{g}_3}^T \end{bmatrix}}_{\boldsymbol{g}} = \boldsymbol{0}_3. \end{aligned}$$

(25)

Now from Eq. (25) the third row is a linear combination of the first two rows of the measurement matrix and promptly we may simplify it as

$$\begin{aligned} \underbrace{ \begin{bmatrix} \boldsymbol{0}^T &{} -\mathcal {H}[\rho _i]^T &{} \mu _{y,i}\mathcal {H}[\rho _i]^T\\ \mathcal {H}[\rho _i]^T &{} \boldsymbol{0}^T &{} -\mu _{x,i}\mathcal {H}[\rho _i]^T \end{bmatrix}}_{M} \underbrace{ \begin{bmatrix} \bar{\boldsymbol{g}_1}^T \\ \bar{\boldsymbol{g}_2}^T \\ \bar{\boldsymbol{g}_3}^T \end{bmatrix}}_{\boldsymbol{g}} = \boldsymbol{0}_3. \end{aligned}$$

(26)

The Eq. (26) is the principle of a well-known method for homography estimation called the Direct Linear Transformation (DLT). As we may observe, the DLT requires at least four correspondences points to estimate the homography matrix G, but we may supply more than four point correspondences and rewrite Eq. (26) as

$$\begin{aligned} \begin{bmatrix} \boldsymbol{0}^T &{} -\mathcal {H}[\rho _1]^T &{} \mu _{y,1}\mathcal {H}[\rho _1]^T\\ \boldsymbol{0}^T &{} -\mathcal {H}[\rho _2]^T &{} \mu _{y,2}\mathcal {H}[\rho _2]^T\\ \vdots &{} \vdots &{} \vdots \\ \boldsymbol{0}^T &{} -\mathcal {H}[\rho _i]^T &{} \mu _{y,i}\mathcal {H}[\rho _i]^T\\ \mathcal {H}[\rho _1]^T &{} \boldsymbol{0}^T &{} -\mu _{x,1}\mathcal {H}[\rho _1]^T\\ \mathcal {H}[\rho _2]^T &{} \boldsymbol{0}^T &{} -\mu _{x,2}\mathcal {H}[\rho _2]^T\\ \vdots &{} \vdots &{} \vdots \\ \mathcal {H}[\rho _i]^T &{} \boldsymbol{0}^T &{} -\mu _{x,i}\mathcal {H}[\rho _i]^T \end{bmatrix} \begin{bmatrix} \bar{\boldsymbol{g}_1}^T \\ \bar{\boldsymbol{g}_2}^T \\ \bar{\boldsymbol{g}_3}^T \end{bmatrix} = \boldsymbol{0}_{2n}, \end{aligned}$$

(27)

which we may be solved through the singular value decomposition method (SVD) [16]. A detailed Maple implementation for homogeneous coordinates transformation and homography estimation by the DLT method can be found in Appendix A.2.

1.2 A.2 Homogeneous Coordinates and DLT Maple Listings

For this Maple listings section, it is necessary to load the following libraries using

B Support Functions Maple Listings

C Numerical Solution for Radial Distortion Maple Listings