I am trying to solve a system of linear equations related to camera calibration as provided in the paper https://downloads.hindawi.com/journals/mpe/2016/1392832.pdf. I am trying to solve Equation number 21 of the paper (described below) but with Polynomial model (also described below). With Polynomial model (Which is given in Equation 5 of the paper and I have described both Division model and Polynomial model in below equations), the equation 21 changes from (Division model) \begin{equation} s\cdot\begin{bmatrix}\frac{\hat{x_d}}{L(\hat{r_d},k)}\\\frac{\hat{y_d}}{L(\hat{r_d},k)}\\1\end{bmatrix} = \begin{bmatrix}\hat{h_1}^TP_w\\\hat{h_2}^TP_w\\\hat{h_3}^TP_w\end{bmatrix} \end{equation} to (Polynomial model) \begin{equation} s\cdot\begin{bmatrix}\hat{x_d}\cdot L(\hat{r_d},k)\\\hat{y_d}\cdot L(\hat{r_d},k)\\1\end{bmatrix} = \begin{bmatrix}\hat{h_1}^TP_w\\\hat{h_2}^TP_w\\\hat{h_3}^TP_w\end{bmatrix} \end{equation} where, $L(\hat{r_d},k) = (1+k_1r_d^2+k_2r_d^4+k_3r_d^6)$ and $k_1, k_2, k_3$ are numbers, $\hat{x_d}, \hat{y_d}$ are pixel coordinates (transformed pixel coordinates to be specific) and I have 1800 pixels, $\hat{r_d}^2 = \hat{x_d}^2 + \hat{y_d}^2$ and $ \hat{h_1^T}, \hat{h_2^T}, \hat{h_3^T} $ are $1\times3$ row vectors of $3 \times 3$ homography matrix and $P_w$ is $3\times1$ column vector i.e. \begin{equation} \hat{h_i^T} = \begin{bmatrix} \hat{h_{i1}}& \hat{h_{i2}}& \hat{h_{i3}} \end{bmatrix} and \; P_w = \begin{bmatrix} P_x \\ P_y \\ P_z \end{bmatrix} \end{equation} In the above equation, $\hat{x_d}, \hat{y_d}, \hat{r_d}, \hat{h_2^T}, \hat{h_2^T}, P_w$ are known and $ k_1, k_2, k_3, \hat{h_3^T}$ are unknowns. Now I will present my derivation of the above modified equation (Polynomial model) in order to solve it with OpenCV Solve() function.
To use OpenCV solve() function, I need to get the above equation to $A\times X = B$ format. I tried the below steps to get the equation into the desired format (However I am not able to get proper results).
Since, $s = \hat{h_3}^TP_w$, I divided s on LHS and $\hat{h_3}^TP_w$ on RHS and Thus I removed the third row of the equation. i.e. \begin{equation} \begin{bmatrix}\hat{x_d}\cdot L(\hat{r_d},k)\\\hat{y_d}\cdot L(\hat{r_d},k)\end{bmatrix} = \begin{bmatrix}\frac{\hat{h_1}^TP_w}{\hat{h_3}^TP_w}\\\frac{\hat{h_2}^TP_w}{\hat{h_3}^TP_w}\end{bmatrix} \end{equation}
After multiplying $\hat{h_3}^TP_w$ on both sides, I will get the below matrix \begin{equation} \begin{bmatrix}\hat{x_d}\cdot L(\hat{r_d},k)\cdot\hat{h_3}^TP_w\\\hat{y_d}\cdot L(\hat{r_d},k)\cdot \hat{h_3}^TP_w\end{bmatrix} = \begin{bmatrix}\hat{h_1}^TP_w\\\hat{h_2}^TP_w\end{bmatrix} \end{equation}
Now I wanted to convert LHS into $A\times X$ format. Since $k_1,k_2, k_3$ and $\hat{h_3}^T$ (=($h_{31} \;h_{32}\; h_{33}$) are unknowns, I tried to take them out by expanding $L(\hat{r_d},k)$ and $\hat{h_3}^T$ out as below: \begin{equation} \begin{bmatrix}\hat{x_d}\cdot (1+k_1r^2+k_2r^4+k_3r^6) \cdot \hat{h_3}^TP_w\\\hat{y_d} \cdot (1+k_1r^2+k_2r^4+k_3r^6) \cdot \hat{h_3}^TP_w\end{bmatrix} \end{equation} Further taking $\hat{x_d}$ inside, I get \begin{equation} \begin{bmatrix}(\hat{x_d}+\hat{x_d}\cdot k_1r^2+\hat{x_d}\cdot k_2r^4+\hat{x_d}\cdot k_3r^6)\cdot \hat{h_3}^TP_w\\(\hat{y_d}+\hat{y_d}\cdot k_1r^2+\hat{y_d}\cdot k_2r^4+\hat{y_d}\cdot k_3r^6)\cdot \hat{h_3}^TP_w\end{bmatrix} \end{equation} Then taking $\hat{h_3}^TP_w$ inside, I get, \begin{equation} \begin{bmatrix}\hat{x_d}\cdot \hat{h_3}^TP_w+\hat{x_d}\cdot k_1r^2\cdot \hat{h_3}^TP_w+\hat{x_d}\cdot k_2r^4\cdot \hat{h_3}^TP_w+\hat{x_d}\cdot k_3r^6\cdot \hat{h_3}^TP_w \\ \hat{y_d}\cdot \hat{h_3}^TP_w+\hat{y_d}\cdot k_1r^2\cdot \hat{h_3}^TP_w+\hat{y_d}\cdot k_2r^4\cdot \hat{h_3}^TP_w+\hat{y_d}\cdot k_3r^6\cdot \hat{h_3}^TP_w \end{bmatrix} \end{equation}
Since $\hat{h_3}^T$ is a $1\times 3$ vector and $P_w$ is $3\times 1$ vector, we can write $\hat{h_3}^TP_w$ as $P_w^T\hat{h_3}$
Using $\hat{h_3}^TP_w$ = $P_w^T\hat{h_3}$ and rearranging the terms I get, \begin{equation} \begin{bmatrix}\hat{x_d}P_w^T\cdot \hat{h_3}+\hat{x_d}r^2P_w^T\cdot k_1\hat{h_3}+\hat{x_d}r^4P_w^T\cdot k_2\hat{h_3}+\hat{x_d}r^6P_w^T\cdot k_3\hat{h_3} \\ \hat{y_d}P_w^T\cdot \hat{h_3}+\hat{y_d}r^2P_w^T\cdot k_1\hat{h_3}+\hat{y_d}r^4P_w^T\cdot k_2\hat{h_3}+\hat{y_d}r^6P_w^T\cdot k_3\hat{h_3} \end{bmatrix} \end{equation} Now, I tookout $ \hat{h_3} $ and $k_1, k_2, k_3 $ and got equation in $A\times X = B$ format as below \begin{equation} \begin{bmatrix}\hat{x_d}P_w^T&\hat{x_d}r^2P_w^T&\hat{x_d}r^4P_w^T&\hat{x_d}r^6P_w^T \\ \hat{y_d}P_w^T&\hat{y_d}r^2P_w^T&\hat{y_d}r^4P_w^T&\hat{y_d}r^6P_w^T \end{bmatrix} * \begin{bmatrix} \hat{h_3}\\ k_1\hat{h_3}\\ k_2\hat{h_3}\\ k_3\hat{h_3} \end{bmatrix} = \begin{bmatrix}\hat{h_1}^TP_w\\\hat{h_2}^TP_w\end{bmatrix} \end{equation}
Here I just want to remind that \begin{equation} P_w^T = \begin{bmatrix} P_x \; P_y \; P_z \end{bmatrix} \end{equation} and \begin{equation} \hat{h_3} = \begin{bmatrix} \hat{h_{31}}\\ \hat{h_{32}}\\ \hat{h_{33}} \end{bmatrix} \end{equation} and since $k_1, k_2, k_3$ are just a scalar quantities, I have taken, \begin{equation} k_1\hat{h_3} = \begin{bmatrix} k_1\hat{h_{31}}\\ k_1\hat{h_{32}}\\ k_1\hat{h_{33}} \end{bmatrix} \end{equation} and so on.
So in my above final equation in $A\times X$ format, A is of $2n \times 12$ matrix where n is number of pixels (1800 in my case) and X is of $12 \times 1$ matrix.
I have formulated the A, B matrix as I have described above and I am using OpenCV solve function with SVD (cv::solve(a, b, x, cv::DECOMP_SVD);)to solve the system of equations.
However, I am not getting a valid results i.e. for example, $k_1$ should be < 1 but I am getting $k_1$ in the order of ten thousands.
Can someone please tell me if my formulation is wrong. If it is wrong, please tell me the correct way of solving such equations.