An Exact Formula for Calculating Inverse Radial Lens Distortions

This article presents a new approach to calculating the inverse of radial distortions. The method presented here provides a model of reverse radial distortion, currently modeled by a polynomial expression, that proposes another polynomial expression where the new coefficients are a function of the original ones. After describing the state of the art, the proposed method is developed. It is based on a formal calculus involving a power series used to deduce a recursive formula for the new coefficients. We present several implementations of this method and describe the experiments conducted to assess the validity of the new approach. Such an approach, non-iterative, using another polynomial expression, able to be deduced from the first one, can actually be interesting in terms of performance, reuse of existing software, or bridging between different existing software tools that do not consider distortion from the same point of view.


Introduction
Distortion is a physical phenomenon that in certain situations may greatly impact an image's geometry without impairing quality nor reducing the information present in the image. Applying the projective pinhole camera model is often not possible without taking into account the distortion caused by the camera lens. This phenomenon can be modelled by a radial distortion, the most prominent component, and a second, with a lesser effect, a decentering distortion which has both a radial and a tangential components. Radial distortion is caused by the spherical shape of the lens, whereas tangential distortion is caused by the decentering and non-orthogonality of the lens components with respect to the optical axis ( [1,2]). It is important to note that radial distortion is highly correlated with focal length [3] even if in literature it is not modelled within the intrinsic parameters of the camera [4]. This is due to the fact that the radial distortion model is not linear, in contrary to other intrinsic parameters. We can see in Figure 1 the displacement applied to a point caused by both radial and tangential distortion.
Decentering distortion was modelled by Conrady in 1919 [5] then remodelled by Brown in 1971 [6] and a radial distortion model was proposed by Brown in 1966 [7]. These distortion models have been adopted by the Photogrammetry as well as the Computer Vision communities for several decades. Most photogrammetric software such as PhotoModeler (EOS) uses these models (see Equations (1) and (2)) to correct observations visible on the images and provide ideal observations. Roughly, radial distortion can be classified in two families, barrel distortion and pincushion radial distortion. Regarding the k 1 coefficient in Formula (1), barrel distortion corresponds to a negative value of k 1 and pincushion distortion to a positive value of k 1 , for an application of the distortion and not a  In the center, a painting from Piet Mondrian [12] (which is now in the public domain since 1 January 2016); on the left, the painting with a barrel effect; and on the right, the same image with pincushion distortion.
Using these models to compensate the observations is now well known and many software dealing with images or panoramas propose plugins dedicated to distortion correction (mainly only radial distortion) [13] . However, although we have the equations to compensate the distortion, how to compute the inverse function in order to apply such a distortion is not obvious. For example, when an image of a known 3D point is computed using a calibrated camera, the 2D projected point can be easily computed, but we need then to apply the distortion to the image point in order to obtain an accurate projection of the original 3D point. This first application justifies the present work. How to determine the inverse of a closed form solution for distortion model equations? The second reason is the merging the work of the two communities involved, photogrammetry and computer vision. While having worked separately for years, the situation between the two communities has drastically changed for almost a decade [14]. This is also visible in the form of new commercial or open-source software dealing with photogrammetry or computer vision. For example, PhotoScan (from Agisoft) or MATLAB that comes with the camera calibration toolbox or, for the open-source side, OpenCV toolbox which provides a solution for multiview adjustment. These three software use the same Equations (1) and (2) to manage distortion, but the mathematical model here is used to apply distortion and not to compensate it. Determining an exact formula to calculate inverse lens distortion, which allows using the same software to apply and compensate distortion with two set of k n parameters can be very useful and, in fact, is the purpose of this work. This paper is organized as follows: in the next section, a demonstration on computing an exact formula for inverse radial distortion is presented. This approach gives a set of k 1 ...k n coefficients computed from the original polynomial distortion model. In Section 3, several applications of this formula are presented. First of all an experiment on this inversion formula is done using only the k 1 ...k n coefficient. Then, an application used to compute a formula for an inverse radial lens distortion is applied to an image coming from a metric camera (Wild P32) with a 3-micron distortion in the edge of the frame. The next experiment is done on a calibration grid built with black disk. Finally, a discussion on converting the distortion model between several photogrammetric software, PhotoModeler (EOS [15]), PhotoScan (Agisoft [16]) and OpenCV [17] is proposed.

Calibration Approach
Radial distortion is mainly considered in the camera calibration process. Since Duane Brown's first publication, a large quantity of work was performed in the field of camera calibration, opening the way for new methods. Several techniques were proposed using orthogonal planes, 2D objects with plannars patterns up to self-calibration with unknown 3D points. Interesting reviews were published on both the photogrammetry and computer vision sides by Fraser [18], Zhang [19] and more recently by Shortis [20].
When Brown [6] proposed a radial distortion model in 1971, he also proposed a way to calibrate cameras using a set of plumb lines. The idea of using a set of straight wires to compute a distortion model in a camera calibration process remains in use 45 years later in the fields of photogrammetry and computer vision; Hartley in 1993 [21,22] then Faugeras and Devernay [23] and recently Nomura [24], Clauss [25], Tardif [26], and Rosten [27].

Inverse Radial Distortion
After Conrady and Brown a lot of work was done to deal with removing distortion from images. As the problem is shared by photogrammetry community as well as computer vision we can refer to many books and papers on this topic. Including the famous Manual of photogrammetry [28] and Atkinson gives an overview of these problems for the two communities [29].
Nevertheless, the problem of reverse distortion is somewhat the poor relation of the problems of distortion. As mentioned by Heikkilä and Silvén [30], "only a few solutions to the back-projection problem can be found in literature, although the problem is evident in many applications." And in the same paper, "we can notice that there is no analytic solution to the inverse mapping".
In the particular case of high distortion as in wide-angle and fish-eye lenses, some non polynomial (and invertible) models have been proposed; for example Basu and Licardie introduced the Fish-Eye Transform (FET) in [31]. Also Faugeras and Devernay [23] propose another invertible model based on the Field-of-View. A complete description of these models can be found in the review written by Hugues [32] and also in [33].
Regarding the polynomial model, several solutions have been tested to perform inverse radial distortion and the solution can be classified in three main classes (even if other approaches can be found such as the use of a neural network [34]): • Approximation. Mallon [35], Heikkilä [30,36] and then Wei and Ma [37] proposed inverse approximations of a Taylor expansion including first order derivatives. According to Mallon and Welhan, "This is sometimes assumed to be the actual model and indeed suffices for small distortion levels." [35]. A global approach, inverse distortion plus image interpolation, is presented in a patent held by Adobe Systems Incorporated [38]. • Iterative. Starting from an initial guess, the distorted position is iteratively refined until a convenient convergence is determined [39][40][41]; • Look-up table. All the pixels are previously computed and a look-up table is generated to store the mapping (as for example in OpenCV).
All these methods involve restrictions and constraints on accuracy, time processing or memory usage.
Nevertheless, some very good results can be obtained. For example, implementing the iterative approach gives excellent results, however the processing time is drastically increased. Given in Peter Abeles' blog [40], the method is easy to implement. Results are shown in Figure 3. Iterative method applied to a Nikon D700 camera with a 14mm lens. Along the frame diagonal the points are first compensated by what is a normal process of the radial distortion using Equation (1) and then the distortion is applied with the iterative process [40] and the result compared to the original point. On Y-axis, the distance between the original point and the computed reverse point.
The iterative solution works by first estimating the radial distortion magnitude at the distorted point and then refining the estimate until it converges.
Algorithm 1 shows an implementation of this approach.

Algorithm 1 Iterative algorithm to compute the inverse distortion
Require: point P n P c = P n repeat r = ||P c || dr = 1 + k 2 r 2 + k 4 r 4 + ... P c = P n /dr until Convergence of P c return P c Only a few iterations are necessary.
The results presented in Figure 3 are in pixels. The used camera here is a Nikon D700 with a 14 mm lens. The calibration was done with PhotoModeler EOS and the results are k 1 = 1.532 × 10 −4 , k 2 = −9.656 × 10 −8 , k 3 = 7.245 × 10 −11 . The coefficients are expressed in millimeters. The center of autocollimation and the center of distortion are close to the image center. Only a few iterations are necessary to compute the inverse distortion. In this case, with a calibration made using PhotoModeler [15], the inverse of distortion represents the application of the distortion to a point projected from the 3D space onto the image.
This iterative approach is very interesting when the processing time is not an issue, as for example in the generation of a look-up table. Note, however, that a good initial value is needed.
According to these existing methods, we now want to obtain a formula for the inverse radial distortion when modelled by a polynomial form as described by Brown. The inverse polynomial form will be an expression of the original k 1 ...k 4 coefficients of the original distortion.

Lens Distortion Models
We consider the general model of distortion correction or distortion removal that can be written in the following form by separating radial and tangential/decentering components:

Radial Distorsion
In the following we will consider only radial distortion.

General Framework
Given a model of distortion or correction with parameters (k 1 , k 2 , k 3 , ...), our general objective is to find the inverse transformation. A natural assumption is to express the inverse transformation on the same form of the direct transformation, i.e., with parameters (k 1 , k 2 , k 3 , ...). Therefore we want to express each k i as a function of all the k j .

Radial distortion
Let us assume that there exists two transformations T 1 and T 2 : where r = x 2 + y 2 , r = x 2 + y 2 , and P and Q are power series: with a 0 = 1, a 1 = k 1 ,..., a n = k n (in order to use k as an index in the calculus of the Appendix).
In addition to starting at n = 0 we facilitate those calculi. We can scale r as r = αr in order to have the same domain of definition for P and Q. So Q reads: with b n = b n α 2n . In the following b n is used instead of b n but we keep in mind this change of variable. Given the definition of r and by using transformation T 2 in Equation (4) we obtain: and similarly with Equation (3): r = r|P(r)| P and Q are positive which allows removing the absolute value. Hence by injecting the last equation in the first we get: and at the end: It is possible to derive a very general relation between coefficients a n and b n but it is not exactly adapted to real situations where P is a polynomial of finite order. Therefore we can derive a slightly simpler relationship in the case where only a 1 , ..., a 4 are given. It is summarized in the following result: Proposition 1. Given the sequence a 1 , ..., a 4 it is possible to obtain the recursive relation: where we use the following intermediate coefficients: We will derive this expression in Appendix A and show how the coefficients b 1 , ..., b n can be computed both with symbolic and numeric algorithms in Appendix B.
Several remarks can be made about this result: The problem is symmetric in terms of P and Q, so the relations found for a n can of course be applied in the reverse order.

Remark 2
For any n the coefficient b n can be computed recursively. In Equation (9), the first summation is obtained thanks to a 0 , ..., a 4 and q(n − 1), ..., q(n − 4) that only depends on the sequence a n . Similarly the second summation involves b 0 , ..., b n−1 and values p(j, 2k) which both depend only the given sequence a n . Therefore the recursive formula for b n can be implemented at any order n. We provide the 4 first terms: All formula till b 9 are summarized in Appendix C.

Results and Experimental Section
In this section we propose three experiments to test this inverse formula for radial distortion in order to evaluate the relevance of such approach.
• First, we begin by testing the accuracy of the inverse formula by applying the forward/inverse formula recursively within a loop. Hence, the inverse of the inverse radial distortion is computed 10,000 times and compared to the original distortion coefficients. • Then, for a given calibrated camera, we compute the residual after applying and compensating the distortion along the camera frame. A residual curve shows the results of the inverse camera in all the frame. • The last experiment is the use of inverse distortion model on a image made with a metric camera built with a large eccentricity (a film-based Wild P32 camera), without distortion. We apply a strong distortion and then compensate it and finally compare it to the original image.

Inverse Distortion Loop
In this experiment, as the formula gives an inverse formula for radial distortion we do it twice and compare the final result with the original one. In a second step, we iterate this process 10,000 times and compare the final result with the original distortion.
The Table 1 shows the original radial distortion and the computed inverse parameters. The original distortion is obtained by using PhotoModeler to calibrate a Nikon D700 camera with a 14 mm lens from Sigma. The results for this step are presented in Table 2. In the columns 'Delta Loop 1' and 'Delta Loop 10000' we can see that k 1 and k 2 did not change and the delta on k 3 and k 4 are small with respect to the corresponding coefficients: the error is close to 1E-10 smaller than the corresponding coefficient. Note that k 4 was not present in the original distortion and as the inverse formula is in function of only k 1 ...k 4 , the loop is computed without the coefficients k 5 ...k 9 which influences the results, visible from k 4 . Table 2. Radial distortion inverse loop and residual between coefficients of the orignal distortion and after n inversions (n = 2 and n = 10,000).

Coefficient
Original This first experiments shows the inverse property of the formula and of course not the relevance of an inverse distortion model. But this experiment shows also the high stability of the inversion process. However, even if coefficients k 1 ...k 4 are sufficient in order to compensate distortion, the use of coefficients k 1 ...k 9 are important for the inversion stability.
The next two experiments show the relevance of this formula for the inverse radial distortion model.

Inverse Distortion Computation onto a Frame
This second experiment uses a Nikon D700 equipped with a 14 mm lens from Sigma. This camera is a full frame format, i.e., a 24 mm × 36 mm frame size. The camera was calibrated using PhotoModeler and the inverse distortion coefficients are presented in Table 1, where Column 1 gives the calibration result on the radial distortion, and Column 2 the computed inverse radial distortion.
Note that the distortion model provided by the calibration using PhotoModeler gives as a result a compensation of the radial distortion, in millimeters, limited to the frame.
The way to use this coefficient is to first express a 2D point on the image in the camera reference system, in millimeters, with the origin on the CoD (Center of Distortion), close to the center of the image. Then the polynomial model is applied from this point.
The inverse of this distortion is the application of such a radial distortion to a point theoretically projected onto the frame.
In all following experiments, the residuals are computed as follows: A 2D point p, is chosen inside the frame, its coordinate are previously computed in millimeters in the camera reference system with the origin on the CoD. Then p 1 is p compensated by the inverse of distortion. Finally p 2 is p 1 compensated by the original distortion.
The residual is the value dist(p, p 2 ) . The following results show the 2D distortion residual curve. For a set of points on the segment [O, max X/2] the residuals are computed and presented in Figure 4 as Y-axis. The X-axis represents the distance from the CoD. These data comes from the calibration process and are presented in Table 1.
The results shown in Figures 4 and 5 are given in pixels.
In Figure 4a, below, we present the residuals using only coefficients k 1 ..k 4 of the inverse distortion. The maximum residual is close to 4 pixels, but residuals are less than one pixel until close to the frame border. This can be used when using non configurable software where it is not possible to use more than 4 coefficients for radial distortion modeling.
In Figure 4, below, we present the residual computed from 0 to max X frame using coeficients k 1 ...k 9 for inverse distortion.
The results are very good, less than 0.07 pixel on the frame border along 0X axis and the performance is quite the same as for compensating the original distortion.
We can see that in almost all the images the residuals are close to the ones presented in Figure 5. Nevertheless, we can observe in Figure 4 higher residual in the corners, where the distance to the CoD is the greatest.  Here follows a brief analysis of the residuals: These two experiments show that the results are totally acceptable even if the residuals are higher in zone furthest from the CoD, i.e., in the diagonal of the frame. As shown in Table 3, only 2.7 % of the frames have residuals > 1 pixel.

Inverse Distortion Computation on an Image Done with a Metric Camera
This short experiment used an image taken with a Wild P32 metric camera in order to work on an image without distortion. The Wild P32 terrestrial camera is a photogrammetric camera designed for close-range photogrammetry, topography, architectural and other special photography and survey applications.
This camera was used as film based, the film is pressed onto a glass plate fixed to the camera body on which 5 fiducial marks are incised. The glass plate prevents any film deformation.
The film format is 65 mm × 80 mm and the focal length, fixed, is 64 mm. Designed for architectural survey the camera has a high eccentricity and the 5 fiducial marks were used in this paper to compute the CoD. In Figure 6a four fiducial marks are visible (the fifth is overexposed in the sky). The fiducial marks are organized as follows: one at the principal point (PP), three at 37.5 mm from the PP (left,right,top) and one at 17.5 mm (bottom).
This image was taken in 2000 in the remains of the Romanesque Aleyrac Priory, in northern Provence (France) [42]. Its semi-ruinous state gives a clear insight into the constructional details of its fine ashlar masonry as witnessed by this image taken using a Wild P32 during a photogrammetric survey.
As this image did not have any distortion, we used a polynomial distortion coming from another calibration and adapted it to the P32 file format (see Table 4). The initial values of the coefficients have been conserved and the distortion polynom expressed in millimeters is the compensation due at any point of the file format. The important eccentricity of the CoD is used in the image rectification: the COD is positioned on the central fiducial marks visible on the images in Figure 6a,b.
(a) (b) Figure 6. Distortion-less image taken using a small-format Wild P32 metric camera and application of an artificial distortion. (a) On the left, original image taken with P32 Wild metric camera; (b) On the right, pincushion distortion applied on this original image whithout interpolation. As images have not the same pixel size some vacant pixel are visible as black lines (see Hughues [32]). These lines surround the distortion center, here located on a fiducial mark, strongly shifted from the image center.
After scanning the image (the film was scanned by Kodak and the result file is a 4860 × 3575 pixel image), we first measure the five fiducial marks in pixels on the scanned image and then compute an affine transformation to pass from the scanned image in pixels to the camera reference system in millimeters where the central cross is located at (0.0, 0.0). This is done according to a camera calibration provided by the vendor, which gives the coordinates of each fiducial mark in millimeters in the camera reference system. Table 5 shows the coordinates of the fiducial marks and highlights the high eccentricity of the camera built for architectural survey. This operation is called internal orientation in photogrammetry and it is essential when using images coming from film-based camera that were scanned. The results of these measurements are shown in Table 5.
In Figure 6a we can see the original image taken in Aleyrac while in the Figure 6b we can see the result of the radial distortion inversion. Figure 7 shows the original image in grey and the image computed after a double inversion of the radial distortion model in green.
We can observe no visible difference in the image. This is correlated with the previous results in the second experiment, see Figure 4. Table 4. Radial distortion comensation and then application of the inverse used with the image taken with the P32 camera.

Conclusions and Discussion
The experiments presented in this article show the relevance of the proposed methodology and the reliability of the result. However, a significant difference exists depending on whether the set of coefficients k 1 ...k 4 or k 1 ...k 9 . For large distortion the number of parameters should be significant. See Figures 4 and 5 for the influence of the number of coefficients. We can note that since the formulation by Brown, the number of coefficients used to characterize the distortion has increased. In 2015 the Agisoft company added k 4 in their radial distortion model while at the same time many software still use only k 1 , k 2 and in 2016 they add p 3 and p 4 to the tangential distortion model.
Even when k 1 ...k 4 are sufficient for compensating the radial distortion, it is however necessary to increase the degree of the polynomial to correctly compute the inverse.

A Bridge between PhotoModeler and Agisoft for Radial Distortion
One of the applications for using such a formula to compute the inverse distortion coefficient in function of k 1 ...k 4 is to convert distortion models between two software programs that use the inverse distortion model, as for example PhotoModeler and PhotoScan from Agisoft. Indeed PhotoModeler uses the Brown distortion model to compensate for observations made on images and so to obtain a theoretical observation without distortion effect. In contrary, PhotoScan from Agisoft uses a similar model but it adds the distortion to a point projected onto the image. To convert a distortion model from PhotoModeler to PhotoScan, or vice versa, we need to compute the inverse distortion model. his implication in the iterative inverse distortion method and Jean-Philip Royer for his implementation of the inverse distortion in Python in PhotoScan Software in order to import and export camera distortion toward other photogrammetric software.
Author Contributions: Pierre Drap designed the research, implemented the inverse distortion method and analyzed the results. Julien Lefèvre proved the mathematical part. Both authors wrote the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In Equation (8)  For k fixed P(r) k can be rewritten as a product: P(r) k = 4 ∑ n 1 =0 a n 1 r 2n 1 ... 4 ∑ n k =0 a n k r 2n k which gives a more compact expression P(r) k = 4k ∑ m=0 r 2m ∑ n 1 +...+n k =m 0≤n i ≤4 a n 1 ...a n k = 4k ∑ m=0 r 2m p(m, k) Note that p(m, k) = 0 as soon as m > 4k. Then we obtain: Let us call t(k) the coefficient in the previous sum. Actually t(k) will turn out to be q(k), defined in Proposition 1. Last we can express the initial equality P(r)Q r(P(r) = 1: To obtain the equality: l ∑ k=0 a k t(l − k) = 0 for l ≥ 1 We decompose this sum: and since a 0 = 1 we conclude that t(l) satisfy the same recurrence relationship as q in Proposition 1. The two quantities are therefore equal since they have the same initial value. We have also an alternative expression for q given by the definition of t depending on b n and p(m, 2n) that gives: By remarking that p(0, 2k) = 1 we obtain Proposition 1.

Appendix B
There are two ways of implementing the result of Proposition 1. One can be interested in having a symbolic representation of coefficients b n with respect to a 1 , a 2 , a 3 , a 4 . But one can also simply obtain numeric values for b n given numeric values for a n .
The main ingredient in both cases is to compute efficiently the coefficients p(j, k). One can start by the trivial result that if n 1 + ... + n k = j then n 1 + ... + n k−1 = j − n k . But n k can take only 5 values between 0 and 4 (there are no other coefficient a n in our case). Therefore one can easily derive the recursive identity: p(j, k) = p(j, k − 1) + a 1 p(j − 1, k − 1) + a 2 p(j − 2, k − 1) + a 3 p(j − 3, k − 1) + a 4 p(j − 4, k − 1) To compute b n we need coefficients p(j, 2k) with the constraints j + k = n, 0 ≤ k and 1 ≤ j ≤ 8k. A table (n − 1) × 2(n − 1) is defined and coefficients are progressively filled by varying the coefficient k from 1 to 2(n − 1) thanks to dynamic programming. This step is summarized in Algorithm B1.
Algorithm B1 Computation of p(j, k) Require: coefficients a 1 , a 2 , a 3 , a 4 , integer N Define an array p of size N × 2N − 1 for k = 0 : 2(N − 1) do p(0, k) = 1 end for for k = 1 : 2(N − 1) do for j = 1 : 2(N − 1) do p(j, k) = p(j, k − 1) + a 1 p(j − 1, k − 1) + a 2 p(j − 2, k − 1) + a 3 p(j − 3, k − 1) + a 4 p(j − 4, k − 1) With p(j, k) = 0 as soon as j < 0 or k ≤ 0 end for end for return p The most tricky aspect is to manipulate formal terms with a 1 , a 2 , a 3 , a 4 . For that it is good to remark that coefficients b n are made of terms a n 1 1 a n 2 2 a n 3 3 a n 4 4 such as n 1 + 2n 2 + 3n 3 + 4n 4 = n. We can therefore have an a priori bound on each exponents, n, n/2, n/3, n/4 respectively. Given that, each multinomial term can be represented as a coefficient in a 4D array of bounded size. It is also very convenient to use a sparse representation for it due to many vanishing terms. Addition of terms are simply additions of 4D arrays of size bounded by n 4 . Multiplication requires shifting operations on dimensions of the array, basically multiplying by a n 1 1 corresponds to a translation by n 1 of the first dimension.