**4. Method**

Our method consists of three stages, which share the same set of calibration images. The first stage is an initial estimation. DIC method [38] is applied on images of speckle pattern calibration targets. Then, using Zhang's approach, a set of control points extracted from DIC result is used for calibration. In the second stage, all parameters and distortions are optimized using a novel object function. In the third stage, distortion rectification mapping is extracted via point-to-point calculation. We will discuss each stage in detail in the following sections.

#### *4.1. Initial Estimation*

Our calibration target is based on a speckle pattern synthesized from Equation (6) [39]. In Equation (6), n and D are the number and radius (unit in pixels) of speckle, respectively. (xk, yk) is the random location of the kth speckle with a random peak intensity of I0k. The synthesized speckle pattern image is shown in Figure 4a, which is denoted as Ir. We printed it as our calibration target. Additionally, we created a mask Im with logical value representation, indicating the scope of the speckle pattern in image Ir, as displayed in Figure 4b.

$$\mathbf{I}(\mathbf{x}, \mathbf{y}) = \sum\_{\mathbf{k}=1}^{n} \mathbf{I}\_{\mathbf{k}}^{0} \exp\left[ -\frac{(\mathbf{x} - \mathbf{x}\_{\mathbf{k}})^{2} + (\mathbf{y} - \mathbf{y}\_{\mathbf{k}})^{2}}{\mathbf{D}^{2}} \right] \tag{6}$$

**Figure 4.** (**a**) Speckle pattern image; (**b**) Mask of speckle pattern image.

With the camera to be calibrated, we captured 15–30 images of this camera calibration target; the ith image is denoted as Idi . We allowed the speckle pattern area to extend beyond the photo's edge. Figure 5 illustrates a calibration target's pose in our calibration image. A rectangle outlines the image with thick solid lines. The array of black points represents control points for initial estimation. It is a noticeable principle that the speckled area can exceed the scope of the image, as shown on the left and bottom of Figure 6, but the array of control points must remain inside the scope of the image.

**Figure 5.** Taking an image of the speckle pattern calibration target.

**Figure 6.** A group of images for DIC calculation in the second stage of our method, containing: (**a**) reprojection of reference image; (**b**) projection of mask; (**c**) image taken by the camera.

For initial estimation, we employ Zhang's camera calibration method. Control points are extracted from the result of the DIC calculation performed on these images. DIC calculation can determine a point-to-point correspondence between points in reference and deformed images. The result of DIC calculation is expressed as displacement of pixels in the deformed image relative to corresponding pixels in the reference image. Equation (7) represents DIC calculation, where Ir is a reference image, Im is the mask of Ir, and Idi is the deformed image. The displacement of all the pixels can be denoted as two mapping matrices, Mxi and Myi , corresponding to X and Y directions, respectively. If we have n deformed images, there are 2n mapping matrices.

$$\mathbf{M}\_{\mathrm{i}}^{\infty}, \mathbf{M}\_{\mathrm{i}}^{\mathrm{y}} = {}^{\circ}\!\!/ \mathrm{d} \mathrm{c} \left(\mathrm{I}^{\mathrm{r}}, \mathrm{I}^{\mathrm{m}}, \mathrm{I}\_{\mathrm{i}}^{\mathrm{d}}\right) \tag{7}$$

By using the DIC approach, we can obtain the displacement of pixels in the speckle pattern area of each Idi relative to the corresponding point in Ir by the DIC method. We took displacement of an array of pixels in Ir and calculated their corresponding subpixel coordinates in a deformed image Idi using displacement and pixel coordinates in Ir, as in Equation (5). These corresponding points were saved as control points.

From initial estimation, we obtained camera parameter **A** and the pose of calibration targets **Ri** and **ti**. Radial and tangential distortion is considered to obtain a more accurate calibration result.

#### *4.2. Optimization with a Novel Objective Function*

At this stage, we performed optimization with a novel objective function, Equations (11)–(13), that is also based on DIC. **A**, **Ri**, and **ti** were used as optimization variables, with initial guess calculated by Zhang's method. Then, we set radial and tangential distortion parameters to zero and reprojected reference image Ir with parameters **A**, **Ri**, and **ti** to obtain "virtual photos" Pr i , as in Equation (8). The mask I m was also projected with the same method, as in Equation (9). Therefore, Pm i is a mask that indicates the scope of the speckle pattern in image Pr i.

$$\mathbf{P}^{\mathbf{r}}\_{\mathbf{i}} = \text{Proj}(\mathbf{A}, \mathbf{R}\_{\mathbf{i}}, \mathbf{t}\_{\mathbf{i}}, \mathbf{I}^{\mathbf{r}}) \tag{8}$$

$$\mathbf{P}\_{\rm i}^{\rm m} = \text{Proj}(\mathbf{A}, \mathbf{R}\_{\rm i}, \mathbf{t}\_{\rm i}, \mathbf{I}^{\rm m}) \tag{9}$$

For every pixel in the projection of reference image Pr i , the DIC method can obtain the displacement of the corresponding point in distorted image I d i taken by the camera, as in Equation (10). A group of these, Pr i , P m i and Id i , is illustrated in Figure 6. It is worth noting that Pr i and P m i share the same estimation of pose corresponding to Id i .

$$\mathbf{M}^{\prime \mathbf{x}}\_{\mathbf{i}\ \mathbf{\prime}}, \mathbf{M}^{\prime \mathbf{y}}\_{\mathbf{i}} = {}^{\prime \mathbf{\prime}} \not\mathbb{C} \left( \mathbf{P}^{\mathbf{r}}\_{\mathbf{i}\ \mathbf{\prime}} \mathbf{P}^{\mathbf{m}}\_{\mathbf{i}}, \mathbf{I}^{\mathbf{d}}\_{\mathbf{i}} \right) \tag{10}$$

Our objective function is set as the square error of these 2n mapping matrices, as in Equations (11)–(13). <sup>Δ</sup>pxu,v,i is element (u, v) of Mx i , meaning *X* direction displacement of point (u, v) in image Id i relative to the corresponding point in image Iproj i . <sup>Δ</sup>pxu,v is average of <sup>Δ</sup>pxu,v,i for every <sup>Δ</sup>pxu,v,i that does not equal to 0, nx is the number of <sup>Δ</sup>pxu,v,i that we took into account. <sup>Δ</sup>pyu,v,iand <sup>Δ</sup>pyu,v are all the same for *Y* coordinate.

$$\min \sum\_{\mathbf{u}, \mathbf{v}} \sum\_{\mathbf{i}} \left[ \left( \Delta \mathbf{p}\_{\mathbf{u}, \mathbf{v}, \mathbf{i}}^{\mathbf{x}} - \overline{\Delta \mathbf{p}\_{\mathbf{u}, \mathbf{v}}^{\mathbf{x}}} \right)^{2} + \left( \Delta \mathbf{p}\_{\mathbf{u}, \mathbf{v}, \mathbf{i}}^{\mathbf{y}} - \overline{\Delta \mathbf{p}\_{\mathbf{u}, \mathbf{v}}^{\mathbf{y}}} \right)^{2} \right], \forall \Delta \mathbf{p}\_{\mathbf{u}, \mathbf{v}, \mathbf{i}}^{\mathbf{x}} \neq 0, \Delta \mathbf{p}\_{\mathbf{u}, \mathbf{v}, \mathbf{i}}^{\mathbf{y}} \neq 0 \tag{11}$$

$$
\forall \Delta \mathbf{p}\_{\text{u,v,i}}^{\text{x}} \neq 0, \Delta \mathbf{p}\_{\text{u,v,i}}^{\text{y}} \neq 0 \tag{12}
$$

$$
\overline{\Delta \mathbf{p}\_{\mathbf{u},\mathbf{v}}^{\mathbf{y}}} = \frac{\sum\_{\mathbf{i}} \Delta \mathbf{p}\_{\mathbf{u},\mathbf{v},\mathbf{i}}^{\mathbf{y}}}{\mathbf{n}\_{\mathbf{y}}}, \forall \Delta \mathbf{p}\_{\mathbf{u},\mathbf{v},\mathbf{i}}^{\mathbf{y}} \neq \mathbf{0} \tag{13}
$$

This objective function means to minimize the difference of displacements between Pr i and I d i . When the optimization process was complete, we obtained new camera parameters and pose of calibration target, namely **A<sup>o</sup>**, **Ro i** , and **t o i** .

#### *4.3. Distortion Rectification Map Extraction*

Distortion rectification maps of *X* and *Y* directions were calculated with reference image Ir, calibration images I d i , and parameters **A<sup>o</sup>**, **Ro i** , and **t o i** , obtained in optimization using a novel objective function.

Using parameters **A<sup>o</sup>**, **Ro i** , and **t o i** , we reprojected reference image Ir and its mask I m to obtain "virtual photos" Tr i and their mask T m i , as in Equations (14) and (15). DIC analysis was performed on every part of Tr i , T m i , and I d i , as in Equation (16). <sup>Δ</sup>po,x u,v,i is element (u, v) of Mo,x i , and <sup>Δ</sup>po,y u,v,i is element (u, v) of Mo,y i . <sup>Δ</sup>po,x u,v is the average of every <sup>Δ</sup>po,x u,v,i that does not equal 0, as in Equation (17). no,x is the number of <sup>Δ</sup>po,x u,v,i that we took into account. <sup>Δ</sup>po,y u,v,i and <sup>Δ</sup>po,y u,v are all the same for *Y* coordinate, as in Equation (18). The result <sup>Δ</sup>po,y u,v, <sup>Δ</sup>po,y u,v is displacement dxu,v, dyu,v used in Equation (5). Therefore, point-to-point mapping of lens distortion is obtained.

$$P\_{\mathbf{i}}^{\mathbf{o},\mathbf{r}} = \text{Proj}(\mathbf{A}^{\mathbf{o}}, \mathbf{R}\_{\mathbf{i}}^{\mathbf{o}}, \mathbf{t}\_{\mathbf{i}}^{\mathbf{o}}, \mathbf{I}^{\mathbf{r}}) \tag{14}$$

$$P\_{\rm i}^{\rm o,m} = \text{Proj}(\mathbf{A}^{\bullet}, \mathbf{R}\_{\rm i}^{\bullet}, \mathbf{t}\_{\rm i}^{\bullet}, \mathbf{l}^{\rm m}) \tag{15}$$

$$\mathbf{M}\_{\rm i}^{\rm o,v}, \mathbf{M}\_{\rm i}^{\rm o,v} = \widetilde{\mathcal{T}}\_{\rm idc}^{\rm r} \left( \mathbf{P}\_{\rm i}^{\rm o,r}, \mathbf{T}\_{\rm i}^{\rm o,m}, \mathbf{I}\_{\rm i}^{\rm d} \right) \tag{16}$$

$$
\overline{\Delta \mathbf{p}\_{\mathbf{u},\mathbf{v}}^{\mathbf{o},\mathbf{x}}} = \frac{\sum\_{\mathbf{i}} \Delta \mathbf{p}\_{\mathbf{u},\mathbf{v},\mathbf{i}}^{\mathbf{o},\mathbf{x}}}{\mathbf{n}\_{\mathbf{o},\mathbf{x}}}, \forall \Delta \mathbf{p}\_{\mathbf{u},\mathbf{v},\mathbf{i}}^{\mathbf{o},\mathbf{x}} \neq \mathbf{0} \tag{17}
$$

$$
\overline{\Delta \mathbf{p}\_{\mathbf{u},\mathbf{v}}^{\mathbf{o},\mathbf{y}}} = \frac{\sum\_{\mathbf{i}} \Delta \mathbf{p}\_{\mathbf{u},\mathbf{v},\mathbf{i}}^{\mathbf{o},\mathbf{y}}}{\mathbf{n}\_{\mathbf{o},\mathbf{y}}}, \forall \Delta \mathbf{p}\_{\mathbf{u},\mathbf{v},\mathbf{i}}^{\mathbf{o},\mathbf{y}} \neq \mathbf{0} \tag{18}
$$
