Self-Calibrated In-Process Photogrammetry for Large Raw Part Measurement and Alignment before Machining

Abstract

Photogrammetry methods are being used more and more as a 3D technique for large scale metrology applications in industry. Optical targets are placed on an object and images are taken around it, where measuring traceability is provided by precise off-process pre-calibrated digital cameras and scale bars. According to the 2D target image coordinates, target 3D coordinates and camera views are jointly computed. One of the applications of photogrammetry is the measurement of raw part surfaces prior to its machining. For this application, post-process bundle adjustment has usually been adopted for computing the 3D scene. With that approach, a high computation time is observed, leading in practice to time consuming and user dependent iterative review and re-processing procedures until an adequate set of images is taken, limiting its potential for fast, easy-to-use, and precise measurements. In this paper, a new efficient procedure is presented for solving the bundle adjustment problem in portable photogrammetry. In-process bundle computing capability is demonstrated on a consumer grade desktop PC, enabling quasi real time 2D image and 3D scene computing. Additionally, a method for the self-calibration of camera and lens distortion has been integrated into the in-process approach due to its potential for highest precision when using low cost non-specialized digital cameras. Measurement traceability is set only by scale bars available in the measuring scene, avoiding the uncertainty contribution of off-process camera calibration procedures or the use of special purpose calibration artifacts. The developed self-calibrated in-process photogrammetry has been evaluated both in a pilot case scenario and in industrial scenarios for raw part measurement, showing a total in-process computing time typically below 1 s per image up to a maximum of 2 s during the last stages of the computed industrial scenes, along with a relative precision of 1/10,000 (e.g., 0.1 mm error in 1 m) with an error RMS below 0.2 pixels at image plane, ranging at the same performance reported for portable photogrammetry with precise off-process pre-calibrated cameras.

1. Introduction

Large parts machining is performed on a near-to-shape raw part, obtained by processes like casting or welding. These raw parts very often do not have any reliable surface or feature reference that can be used for in-machine alignment. However, initial alignment of the part at the machine is a critical process, since an incorrect alignment will give rise to material shortage, which is associated with spoiling the part or with a costly recovery process. Due to the high cost associated with the rejection of a part, the initial alignment is usually done by way of long, time consuming manual processes.

Photogrammetry methods are being used more and more as a 3D technique for large scale applications [1,2,3,4]. To overcome the aforementioned limitations of the large raw part manual alignment, a photogrammetry based process was developed in previous work [5]. The process consists of two steps (Figure 1): (a) Out-of-machine measurement by photogrammetry of optical targets on raw part surfaces, and its mathematical orientation (fitting) to the ideal designed part frame, where the optimal location of optical targets is obtained to use them as feature references in the machine, and (b) In-machine measurement of such references by a machine-integrated measurement system, a special purpose machine vision or a spindle-integrated contact probe, where the in-machine location and orientation of the raw part is measured, assisting the machine operator through an efficient alignment and part fixturing process prior to part machining.

Two main limitations were observed for the out-of-machine portable photogrammetry presented in [5]: First, since a post-process ray-net bundle adjustment [6] was adopted, long computing times were needed for solving the 3D target coordinates after taking an image set, being in the range of some minutes. Due to the post-process computing approach, the lack of adequate images for having a solvable and consistent multiple view geometry or the diagnosis of pending optical targets to measure, could only be performed after all images were taken [7,8]. This resulted in a time consuming and user dependent iterative review and re-process measuring procedure until an adequate set of images was completed. As a result, the measuring process was clearly far away from the potential that photogrammetry may bring to easy, low-cost and fast industrial metrology for large components [9]. Better in-process computational efficiency was required for solving the bundle adjustment, so that images taken together thus far with solved 3D target coordinates could be checked-out in quasi real time for a reliable measuring process, making the most of the potential of portable photogrammetry in large scale applications.

Second, the system was thoroughly tested in a pilot case scenario showing measuring errors in the range of 1 mm for 1 m long parts (1/1000 relative precision), mainly due to the lack of precise calibration of the camera. Better accuracy was necessary for gaining performance when measuring raw parts larger than 1 m where the errors would increase proportionally to the size, as well as for applications where tighter raw material clearances may be expected. Indeed, camera calibration is a relevant uncertainty contributor in photogrammetry [2,10,11], and precise measurements require precise camera model calibration (so called intrinsic parameters). Continuous and recent improvements regarding the mathematical analysis of the data taken for calibration [11,12,13] and image distortion correction [14,15] can be found in the literature.

In the present work, a new efficient procedure is presented for solving the bundle adjustment problem in portable photogrammetry. The presented development demonstrates efficient in-process bundle computing capability on a consumer grade desktop PC, enabling quasi real time 2D image and 3D scene diagnosis so that a reliable measuring procedure can be conducted by portable photogrammetry, avoiding user dependent and inefficient measuring procedures. Additionally, based on the efficient computational approach, an in-process self-calibration method is implemented using the solved 3D point cloud scene itself as calibration geometry, avoiding the need for specialized off-process pre-calibrated cameras or the use of special purpose calibration artifacts in the scene for precise measurements in an industrial metrology application, where traceability is set only by the scale bars available in the measuring scene.

The developed self-calibrated in-process photogrammetry has been evaluated both in the pilot case scenario of [5] (1.5 m long reference part) and at two large scale industrial scenarios (up to 15 m long raw parts), showing successful results with an in-process computing time ranging around 1 s per image and a relative precision of 1/10,000, which is at the precision level reported for portable photogrammetry with precise pre-calibrated cameras.

Section 2 below describes the raw part measuring process by photogrammetry, along with the adopted method for evaluating measurement accuracy both in the pilot case and industrial scenarios. Section 3 introduces the multiple view geometry in portable photogrammetry, along with the non-linear multivariable optimization problem to compute the bundle adjustment. Section 4 presents the implemented in-process image and target computing procedure, and Section 5 describes the self-calibration functionality development, where computing efficiency and precision performance are demonstrated for the pilot case under study, respectively. Section 6 shows the evaluation results for large scale industrial scenarios. Finally, main conclusions and further steps are described in Section 7.

2. Materials and Methods

Regarding precision in portable photogrammetry, the uncertainty determination of optical measurement methods is still an issue [16]. Analytical calculated uncertainty budget of photogrammetric systems according to GUM [17] using error propagation theory is a challenge and consistent approaches have not been found in the literature. The main reason is the large number of undefined error sources affecting the process, especially when not working in laboratory conditions [1,2]. Variability simulation is an alternative to uncertainty budgeting for complex measurement tasks. However, it also requires deep knowledge of the measurement chain and the statistical distribution of each influence quantity [18], which is also a challenging issue in photogrammetry [19]. Therefore, despite the existence of well-accepted guidelines and standards, the definition of accuracy performance is still very heterogeneous [20]. However, industrially accepted assessment guidelines such as the VDI 2634 [21] offer an alternative for evaluating photogrammetric measurement performance, where the recommended procedure is to arrange a set of calibrated scale bars (control bars) around the object scene and through the principal coordinate axes. Lengths measured by photogrammetry are compared to their calibrated values and corresponding length measuring errors (LME) are then evaluated.

Photogrammetry is the core technique used in the procedure presented in [5] for raw part alignment into a machine tool prior to its machining (Figure 1). The raw part is measured by photogrammetry, using retroreflective coded and non-coded optical targets. Non-coded targets (Figure 2) are used to measure the raw part surfaces to be machined. Auxiliary coded targets are used to assist the measuring process by portable photogrammetry so that consistent image-taking is performed during measurement, enabling the correspondence matching between the information given at consecutive images. Images are taken around the part by a digital camera, and transferred in-process to a desktop PC via wireless communication. Images are computed and optical target 3D coordinates are obtained. After measuring, non-coded target coordinates are fitted to the desired nominal part geometry (defined in a CAD file in this work), using a specific purpose computing criteria for fitting [5], enabling even and positive overstock distribution over the surfaces to be machined. As a result, available overstock over all raw part surfaces can be controlled prior to part set-up into the machine. After out-of-machine measurement and fitting finishes, a subset among the optical target 3D coordinates are selected and used as a reference for raw part in-machine alignment. An in-machine measuring system is then used for measuring the selected target coordinates in the machine frame, and the raw part fixturing is adjusted so that its alignment to the machine axes corresponds with its computed optimal by fitting.

Figure 3 shows the verification scenario adopted in this work in order to evaluate the precision of the developed self-calibrated portable photogrammetry. A reference part is adopted as a pilot case (Figure 3a), resembling the measuring scenario at [5] where limited relative precisions of 1/1000 were reported for portable photogrammetry, with four prismatic steel subelements screwed to a mechano-welded structure and milled to a nominal geometry (Figure 3b). A length measuring error (LME) approach similar to VDI 2634 has been adopted for evaluating measuring procedure uncertainty in the developed self-calibrated portable photogrammetry, verifying the performance of the corresponding measurands (i.e., optical target 3D coordinates). Six scale bars (L1 to L6) were set around the scene (1.5 m × 1 m × 0.5 m) in two groups of three (L1 to L3, and L4 to L6), along with a scale bar in the center of the scene (L0), up to a total of seven bars in order to analyze uncertainty in non-coded target measuring (Figure 3c). Each group of three was located in opposite corners of the scene and each scale bar in a group was aligned such that the length components were set in all three spatial directions. For avoiding the influence of length calibration uncertainty between different scale bars in the LME analysis, only one and same scale bar, a calibrated carbon fiber bar with two non-coded targets separated a length L_scale = 1340.099 mm, was used and moved around the scene to each location. Additionally, as an alternative way of evaluating measuring performance, prismatic subelements (Figure 3d) have been measured placing non-coded targets on each milled surface. 3D coordinates of the computed non-coded targets were fitted to its nominal geometry, and the errors corresponding to each nominal reference surface were evaluated.

Eventually, the self-calibrated in-process photogrammetry was evaluated in two industrial scenarios with up to 15 m long raw parts, having as a reference a spindle integrated contact probe measuring process (repeatability ranging 1 μm) executed in relatively accurate milling machines (assumed uncertainty ranging at a typical value of 0.01 mm). Raw part surfaces to be machined were measured and non-coded target coordinates fitted [5] to its nominal. A subset among the fitted target coordinates (control points) were then used as a reference for raw part in-machine alignment, measured by a machine-spindle integrated contact probe during the alignment process. Once the part was properly aligned, machine coordinates of optical targets not used as control points for alignment (i.e., check points) were measured by the contact probe. In order to evaluate the measuring performance of the developed self-calibrated photogrammetry, errors were evaluated between the gauged values by the contact probe and the corresponding 3D coordinates obtained after measurement and fitting.

3. Measurement by Portable Photogrammetry

Photogrammetry has been defined by the American Society for Photogrammetry and Remote Sensing (ASPRS) [23] as the art, science, and technology of obtaining reliable information about physical objects and the environment through processes of recoding, measuring and interpreting photographic images and patterns of recorded radiant electromagnetic energy and other phenomena. In precise industrial metrology applications, it calculates the 3D coordinates of a discrete number of optical targets on an object [20,24] based on their 2D image coordinates. Several photographs of an object are obtained from different positions of the camera. By detection of optical target 2D coordinates at different images it is possible to calculate the positions and orientations of the camera for each of the photographs (so called extrinsic parameters), and also to calculate the 3D coordinates of the target points by multiple view triangulation. The measuring frame is given by a predefined set of coded optical targets with initially known 3D coordinates (Figure 2) and process traceability is set by appropriate scale bars (Figure 2) with calibrated lengths between corresponding pairs of optical targets.

3.1. Multiple-View Geometry

The geometry of a general photogrammetric problem is illustrated in Figure 4. A set of i = 1 … N targets with 3D coordinates $X_i$ are projected in the image plane p_ij from j = 1 … M different points of view. Each camera position has associated extrinsic parameters, given by the rotation matrix R_j defined by the rotations around the reference frame axes x, y, z-based on the α_j, β_j and γ_j Euler angles- and the translation vector d_j. The main goal of portable photogrammetry is the joint determination of the 3D coordinates of the N targets ($X_i$) along with the location and orientation (camera extrinsic parameters, R_j and d_j) of the M image camera frames.

Each target 3D coordinates $X_i={\left[x_i y_i z_i\right]}^T$ can be expressed as $U_{ij}={\left[u_{ij} v_{ij} w_{ij}\right]}^T$ in each camera frame depending on its extrinsic parameters R_j and d_j as:

$$U_{ij}=R_j{X}_i+d_j$$

(1)

Each target 3D coordinate $U_{ij}$ can be projected into the corresponding camera 2D image plane as $p_{ij}$ and $q_{ij}$ coordinates (Figure 5), following the widely assumed pin-hole conic projection model in machine vision [24] as:

$$\begin{array}{c}\left[\begin{array}{c}{p}_{ij}\\ q_{ij}\end{array}\right]\end{array}=f\left[\begin{array}{c}\raisebox{1ex}{$u_{ij}$}\!\left/ \!\raisebox{-1ex}{$w_{ij}$}\right.\\ \raisebox{1ex}{$v_{ij}$}\!\left/ \!\raisebox{-1ex}{$w_{ij}$}\right.\end{array}\right]$$

(2)

being f the focal distance of the camera lens.

The widely assumed Brown’s model [1,16,25] is adopted for modelling camera and lens distortions, which can be expressed as follows:

$$\begin{gathered} \left[\begin{array}{c}\widetilde{p_{ij}}\\ \widetilde{q_{ij}}\end{array}\right]=\left[\begin{array}{c}{h}_{ij}\left(1+k_1{r}_{ij}^2+k_2{r}_{ij}^4\right)+{\pi }_1\left(r_{ij}^2+2h_{ij}^2\right)+2{\pi }_2{h}_{ij}{v}_{ij}\\ v_{ij}\left(1+k_1{r}_{ij}^2+k_2{r}_{ij}^4\right)+{\pi }_2\left(r_{ij}^2+2v_{ij}^2\right)+2{\pi }_1{h}_{ij}{v}_{ij}\end{array}\right] \\ {r}_{ij}=\sqrt{{\left(h_{ij}-h_0\right)}^2+{\left(v_{ij}-v_0\right)}^2} \end{gathered}$$

(3)

where ($h_{ij}, v_{ij}$) are the 2D detected target coordinates in the image plane and $\left(\begin{array}{c}\widetilde{p_{ij}, }\end{array}\widetilde{q_{ij} }\right)$ are the undistorted or corrected image coordinates, given $k_1$ and $k_2$ as the first and second order coefficient for modelling the radial distortion, ($h_0$,$v_0$) as the center of distortion and ${\pi }_1$ and ${\pi }_2$ as the first order tangential distortion coefficients.

The camera and lens distortion model, along with the focal distance, form the so called camera intrinsic parameters characterizing the camera condition for multiple view geometry. The camera intrinsic parameters can be obtained by previous off-line calibration processes [26,27].

Reprojection residual errors $r_{p_{ij}}=\left(p_{ij}-\widetilde{p_{ij}}\right)$ and $r_{q_{ij}}=\left(q_{ij}-\widetilde{q_{ij}}\right)$ can then be defined for every target observed at every image (Figure 5) as the difference between the corrected target coordinates ${\tilde{p}}_{ij}$ and ${\tilde{q}}_{ij}$ and the projected target coordinates $p_{ij}$ and $q_{ij}$, which directly depend on the 3D coordinates of the targets ($X_i$) and camera extrinsic parameters (R_j and d_j) to be solved. As a result, given by the camera intrinsic parameters and the multiple view geometry, target 3D coordinates and camera extrinsic parameters can be jointly computed so that residual errors are minimized for every target at every image, leading to the so-called bundle adjustment [24] for solving the reprojection optimization problem in portable photogrammetry.

3.2. Optimization Problem

The problem is thus defined as a nonlinear multivariate overdetermined system, which can be solved as a least square optimization problem or bundle adjustment by numerical methods such as Levenberg-Marquardt or Gauss-Newton [6,20,24,25]. An initial approximation is defined for the set of variables to solve (α_j, β_j, γ_j, and d_j from j = 1 … M camera views, and $X_i$ from i = 1 … N optical targets) and a numerical iteration procedure is applied to compute them so that residual errors ($r_{p_{ij}}$ and $r_{q_{ij}}$) are minimized.

Residual errors can be grouped into a single residual vector r according to the following structure:

$$\begin{gathered} r_j=\left[\begin{array}{c}{r}_{11}\\ r_{11}\\ ⋮\\ r_{1Nj}\\ r_{2Nj}\end{array}\right]=\left[\begin{array}{c}{p}_1-\widetilde{p_1}\\ q_1-\widetilde{q_1}\\ ⋮\\ p_{Nj}-\widetilde{p_{Nj}}\\ q_{Nj}-\widetilde{q_{Nj}}\end{array}\right] \\ r=\left[\begin{array}{c}{r}_1\\ ⋮\\ r_M\end{array}\right] \end{gathered}$$

(4)

where each r_j vector corresponds to the reprojection error distribution at the j image, being N_j the number of targets detected in each image.

There are two main sets of parameters to be solved-extrinsic and target coordinates, represented by θ_e and θ_x, respectively, which can be grouped in a common vector θ.

The extrinsic parameters vector θ_e is represented as follows,

$${\theta }_e={\left[{\alpha }_1 {\beta }_1 {\gamma }_1 d_1^T \dots {\alpha }_M {\beta }_M {\gamma }_M d_M^T\right]}^T$$

(5)

and the target coordinates vector θ_x is like

$${\theta }_x={\left[x_1 y_1 z_1 \dots x_N y_N z_N \right]}^T$$

(6)

Based on these definitions, the complete vector θ is defined as

$$\theta =\left[\begin{array}{c}{\theta }_e\\ {\theta }_x\end{array}\right]$$

(7)

Following the above definitions, the process consists in the iteration of vector $\theta$ (${\theta }_k\leftarrow {\theta }_{k-1}+{∆}_{\theta }$) towards an optimal value $\widehat{\theta }$ which minimizes the residual vector ${\Vert \overrightarrow{r}\Vert }^2$ norm. According to Gauss-Newton method [28], each iteration ${∆}_{\theta }$ can be computed as a linear system by

$$J^{T}J{∆}_{\theta }=-J^{T}r$$

(8)

where J is the Jacobian matrix containing the partial derivatives of each component of the residual vector with respect to the parameters to optimize in $\theta$.

$$J\left(\theta \right)=\left[\begin{array}{ccc}\frac{\partial r_1}{\partial {\theta }_1}& \dots & \frac{\partial r_1}{\partial {\theta }_{6M+3N}}\\ ⋮& \dots & ⋮\\ \frac{\partial r_{2m}}{\partial {\theta }_1}& \dots & \frac{\partial r_{2N}}{\partial {\theta }_{6M+3N}}\end{array}\right]$$

(9)

with $m={\displaystyle \sum }_{j=1}^M{N}_j$ being the total number of targets detected at all images.

4. In-Process Computing Procedure for Time Efficiency

The computing performance of the bundle adjustment method plays a relevant role in the development of a time-efficient in-process strategy for portable photogrammetry. The convergence time of the numerical method depends on the size of the Jabobian matrix $J$, due to matrix element allocation, assignation and computation times for managing relatively large and sparse $J^{T}J$ and $J^T$ matrices, and due to the conditioning and size of the ${\left(J^{T}J\right)}^{-1}$ inverse matrix calculation so that ${∆}_{\theta }$ is determined in each iteration. As an example, given a measuring scenario with 100 images and 100 targets, assuming all targets are detected at every image, $J^{T}J$ matrix size is 900 × 900, with $J$ being 40,000 × 900. Indeed, this problem increases with the number of images and targets included in the joint bundle during the measuring process. An alternative for reducing this computational work is the decomposition of the linear system to solve ${∆}_{\theta }$ (Equation (8)) to a set of a lower range and individually solved linear subsystems. In this work, the reprojection error partial derivatives forming the Jacobian matrix $J$ (Equation (9)) are analytically expressed and system decomposition is adopted so that ${∆}_{{\theta }_e}$ and ${∆}_{{\theta }_x}$ are individually solved for interdependent extrinsic and target coordinate iteration, avoiding direct computation of ${∆}_{\theta }$ in order to increase computational efficiency.

On the other hand, the total number of iterations (k) and the corresponding convergence time depends on the initial adopted approximation for every variable in $\theta$ to optimize in the joint bundle minimization problem. Furthermore, the numerical method may diverge if a too-inaccurate initial approximation is adopted. So, along with the former system decomposition approach, the adoption of adequate and robust initial approaches for new camera extrinsic and target coordinates every time a new image is taken can avoid unnecessary joint computation effort every time a new image is taken.

Figure 6 shows a schematic view of the in-process procedure developed. The following steps are conducted each time a new image is taken during the measuring process by portable photogrammetry:

• Image processing proceeds for optical target detection (h_ij, v_ij) and correction $(\widetilde{p_{ij}},\widetilde{ q_{ij}}$), along with decoding of coded targets [24,29].
• Computation of an initial approach is performed for the camera extrinsic parameters of the new image (α_j, β_j, γ_j, and d_j), according to already solved coded optical target 3D coordinates detected on the image (Figure 6a).
• In case a minimum set of three coded optical targets with known coordinates is not available in the new image, the camera extrinsic is not computed and the procedure stops asking for a new image having a minimum set of targets to proceed back in step 1.
• Given by the new camera extrinsic, code assignation is performed to non-coded targets with still unsolved correspondences in all images captured so far, following the so called Hungarian method [30,31].
• Given by the new camera extrinsic parameters, computation of an initial approach is performed for new target 3D coordinates $X_i$ unsolved so far but coded, provided that each one is jointly observed by a minimum set of two camera views with known extrinsic parameters (Figure 6b).
• An intermediate joint bundle adjustment (Figure 6c) is conducted for the camera extrinsic and target coordinates solved according to all images so far. Given by the initial approaches in steps 2 and 5, only one bundle iteration is performed, so that a sufficiently accurate and consistent epipolar net construction is obtained every time a new image is included in the minimization problem, ensuring a reliable correspondence solving for non-coded targets in step 4, and avoiding unnecessary computational work until joint bundle convergence at this step. The measuring process can now continue with the acquisition of new images, computed in-process from step 1 to 6 every time a new image is taken.
• Finally, once the measuring process finishes, the post-process joint bundle of camera extrinsic parameters and target coordinates is computed until convergence and measuring process traceability is set by calibrated scale bar distances available at the scene, where measuring frame target coordinates are also included into the bundle adjustment.

Regarding process reliability, step 3 enables in-process control of the information provided by each image, so that measuring process by portable photogrammetry can be guided ensuring a set of images that will form a consistent and geometrically determined ray-net, avoiding inefficient post-process evaluation and iterative processes [7,8].

The main steps of the presented in-process strategy are described below. Results are shown for the pilot case under study (Figure 3). Section 4.1 and Section 4.2 show the initial computation approach of the camera’s extrinsic parameters (step 2, Figure 6a) and target coordinates (step 5, Figure 6b) as independent computing problems, respectively. Section 4.3 and Section 4.4 describe the joint bundle adjustment problem both before (step 6, Figure 6c) and after, including scale bars for traceability (step 7), respectively. Special attention is paid to the particular analytic expression of Equation (9) corresponding to each step of the in-process procedure. Although linear closed form expressions can be given in Section 4.1 and Section 4.2 for the independent camera extrinsic parameters and target coordinate computing [24], nonlinear approaches are shown as intermediate steps towards presenting the submatrix decomposition of Equation (8) for the nonlinear joint bundle computation described in Section 4.3 and the following. Optical target image coordinate detection ($h_{ij}$,$v_{ij}$) and decoding of coded targets is performed following [29].

4.1. Camera Extrinsic Parameters Initial Approach Computation

Computation of the extrinsic parameters of a j-th camera view (Figure 6a) means determining vector ${\theta }_{e_j}$, where camera principal frame location ($d_x$, $d_y$, and $d_z$) and orientation ($\alpha$, $\beta$, and $\gamma$ Euler angles defining a corresponding rotation matrix $R_j$ as $R_{\gamma }{R}_{\beta }{R}_{\alpha }$) are defined with respect to a measuring frame (Figure 4) determined by a pre-set of optical coded targets. The detection of a minimum set of 3 targets is required in order to geometrically determine image extrinsic computing.

$${\theta }_{e_j}={\left[\alpha \beta \gamma d_x d_y d_z\right]}^T$$

(10)

As inputs for the method, X_i coordinates of the set of coded reference targets are known, as well as camera intrinsic parameters (focal distance, and distortion model parameters for Equation (3)). As a result, partial derivatives (Equation (9)) in a Jacobian matrix $E_j$ can be defined as follows (Equation (11)) for solving camera extrinsic parameters for the j-th image as a nonlinear reprojection minimization problem.

$${\left(E_j\right)}_{2N_j\times 6}=\left[\begin{array}{c}{J}_{e_{1j}}\\ ⋮\\ J_{e_{Nj}}\end{array}\right]$$

(11)

where each submatrix $J_{e_{ij}}$ contains the partial derivatives of the projection errors $r_{p_{ij}}=\left(p_{ij}-\widetilde{p_{ij}}\right)$ and $r_{q_{ij}}=\left(q_{ij}-\widetilde{q_{ij}}\right)$ of the i-th target X_i coordinates with respect to the j-th image camera extrinsic, to a total of N_j optical targets detected in the that j-th image. Each $J_{e_{ij}}$ submatrix can be expressed as

$${\left(J_{e_{ij}}\right)}_{2\times 6}={\left(E_{ij}\right)}_{2\times 6}=D_P{D}_{U_E}$$

(12)

where D_P contains the partial derivatives of the i-th projected target coordinates $p_{ij}$ and $q_{ij}$ in the j-th image with respect to $U_{ij}$ target coordinates (Equation (1)) in the corresponding j-th camera principal frame, given as follows according to Equation (2)

$$D_P=\left[\begin{array}{ccc}\frac{1}{w_{ij}}& 0& -\frac{u_{ij}}{w_{ij}^2}\\ 0& \frac{1}{w_{ij}}& -\frac{v_i}{w_{ij}^2}\end{array}\right]$$

(13)

and where $D_{U_E}$ expresses the partial derivatives of the $U_{ij}$ target coordinates with respect to the j-th camera extrinsic in ${\theta }_{e_j}$ as

$$\begin{gathered} D_{U_E}=\left[ D_A{X}_i D_B{X}_i D_C{X}_i I_{3\times 3}\right] \\ {D}_A=R_{\gamma }{R}_{\beta }{D}_{\alpha } \\ {D}_B=R_{\gamma }{R}_{\beta }{R}_{\alpha } \\ {D}_C=D_{\gamma }{R}_{\beta }{R}_{\alpha } \end{gathered}$$

(14)

given the partial derivatives $D_A$, $D_B$, and $D_C$ of the camera principal frame rotation matrix R with respect to each Euler angle $\alpha$, $\beta$, and $\gamma$, respectively.

At the beginning of the measuring process by portable photogrammetry, first images are taken to the set of coded targets determining the measuring frame. Given that a minimum number of three targets are detected in an image, the corresponding individual iteration of $∆{\theta }_{e_j}$ and independent extrinsic computing can proceed for each image. Figure 7 shows an example of computation of the extrinsic initial approach for the first three images of a measuring process conducted on the pilot case evaluation scenario (Figure 3). Nominal values are adopted as an initial approach for the intrinsic parameters, with focal distance f being 24.0 mm and first order radial distortion coefficient $k_1$ being 5.0 × 10⁻⁹ pixel⁻², where higher order radial distortion coefficients are neglected, along with decentering and tangential distortion coefficients in Equation (3). At this stage, the initial adopted approach for the intrinsic parameters is not intended to be accurate but only a nominal value close enough to the expected precise one, to be computed by the in-process self-calibration (see Section 5) so that the final bundle converges in a stable and efficient way.

Table 1. Input data for extrinsic computing: 3D coordinates (x, y, and z in mm, for X₁ to X₆) for coded targets on reference frame, and image coordinates (h_ij and v_ij in pixels, for image j = 1 to 3).

	Reference Frame Targets (Xi)			Image Coordinates (${\mathit{h}}_{\mathit{i}\mathit{j}}$,${\mathit{v}}_{\mathit{i}\mathit{j}}$)
	x	y	z	hi1	vi1	hi2	vi2	hi3	vi3
X1	0	0	0	−51.652	21.593	−21.713	−24.392	24.592	37.326
X2	−169.963	2.650	−0.356	−696.361	27.686	−10.666	−574.886	18.178	637.234
X3	170.036	0	0	594.253	7.982	−26.418	528.265	23.846	−561.086
X4	−1.742	−169.186	0	−52.039	541.249	−448.374	−41.592	504.455	37.142
X5	−0.162	26.998	145.558	−64.312	−473.212	403.334	−25.024	−431.115	32.518
X6	0.109	26.590	28.314	−54.508	−125.784	104.393	−22.956	−112.206	36.936

presents the known 3D coordinates (X_i, i = 1 ..6) of the detected coded targets on each image, along with the detected image coordinates for each target at each image ($h_{ij}$,$v_{ij}$, j = 1 .. 3).

Table 2. Computed results for independent camera extrinsic approaches ($d_x$, $d_y$, and $d_z$ in mm, and $\alpha$, $\beta$, and $\gamma$ in radians), given by the input data in Table 1, along with the optimization quality index for each minimized reprojection error vector (RMS in pixels).

	dX	dY	dZ	α	β	γ	RMS
Image 1	−13.552	5.620	1145.020	−2.375	0.005	0.020	0.482
Image 2	−6.593	−7.545	1340.136	−2.348	−0.009	−1.580	0.454
Image 3	6.894	10.494	1233.812	−2.378	0.023	−4.712	0.471

shows the computed extrinsic parameters for each camera view, according to Gauss-Newton method following Equation (8) for ${\theta }_{e_j}$ iteration. The extrinsic computing performance can be observed in the resulting RMS value of the reprojection error vector after convergence, ranging at 0.5 pixels. Regarding computing efficiency, a mean time lower than 0.1 ms is observed for each iteration.

4.2. Target 3D Coordinate Initial Approach Computation

Once the camera extrinsic parameters are available (Table 2), 3D coordinates of new targets other than those defining the measuring frame can be computed. Analog to extrinsic computing in Section 4.1, the computation of the 3D coordinates of a i-th new target (Figure 6b) means determining vector ${\theta }_{x_i}$, where X_i coordinates are defined in the same measuring frame at which camera extrinsic parameters are known. The detection of the target in a minimum set of 2 images is required in order to geometrically determine target coordinate computing by triangulation.

$${\theta }_{x_i}=X_i={\left[x y z\right]}^T$$

(15)

As inputs for the method, camera intrinsic and extrinsic parameters (α_j, β_j, γ_j, and d_j) of a minimum set of images are known. As a result, partial derivatives (Equation (9)) in a Jacobian matrix $D{X}_i$ can be defined as follows (Equation (16)) for solving 3D coordinates for the i-th target.

$${\left(D{X}_i\right)}_{2M_i\times 3}=\left[\begin{array}{c}{J}_{x_{1i}}\\ ⋮\\ J_{x_{Mi}}\end{array}\right]$$

(16)

where each submatrix $J_{x_{ij}}$ contains the partial derivates of the projection errors $r_{p_{ij}}$ and $r_{q_{ij}}$ of the i-th target with respect to its X_i coordinates, to a total of M_i images in which the i-th target is detected. Each $J_{x_{ij}}$ submatrix can be expressed as

$$J_{x_{ij}}={\left(X_{ij}\right)}_{2\times 3}=D_P{D}_{U_X}$$

(17)

where D_P was defined in Equation (13), and $D_{U_X}$ expresses the partial derivatives of the $U_{ij}$ target coordinates at the j-th camera frame with respect to its $X_i$ coordinates at the common measuring frame as

$$D_{U_X}=R_j$$

(18)

being $R_j$ the rotation matrix corresponding to the j^th camera frame.

Following the example in Section 4.1, given the computed extrinsic parameters of the set of three camera views (Table 2), the iteration of $∆{\theta }_{x_i}$ and the computation of the X_i coordinates proceeds for coded targets detected at least in two camera views. Figure 8 shows an example of computation of a coded target detected in all three images (Figure 8a–c). Again, nominal values are adopted for the camera intrinsic parameters. Figure 8d depicts the 3D location computed for that new target in the same measuring frame as camera views are given.

Along with computing the new target X_i coordinates, new camera extrinsic parameters also enable code assignation to non-coded targets with pending correspondence to solve between detected coordinates in different images. The set of 3D epipolar lines given by the non-coded targets detected in all images can be expressed as 2D epipolar lines in each image. The distances from each 2D epipolar line to each non-coded target detected in the image can be calculated, so that the subset of 2D epipolar lines closest to a specific non-coded target in that image are likely to correspond to the epipolars of that 3D target detected in the other images. All possible 2D epipolar-to-target distances can be computed according to all possible correspondences between non-coded targets detected in three consecutive images, and the Hungarian method [24,29] can be adopted for solving the correspondences as the optimal combination which minimizes epipolar-to-target distance distribution among all 2D epipolars and non-coded target coordinates. As a result, codes can be assigned to non-coded targets in each image, and their 3D coordinate can be correspondingly computed.

Figure 9a shows all the detected targets in the image presented in Figure 8c, including both coded targets and non-coded targets with solved correspondences between the three images. Figure 9b shows the resulting scene including the 3D coordinates of the computed new targets, along with the 3D epipolar net used to solve them. Computed target coordinates are listed in

Table 3. Computed results for independent target coordinate initial approaches, given by the input data in Table 2, along with the optimization quality index for each minimized reprojection error vector (RMS in pixels).

Target Coordinates (Xi)				Image Coordinates (${\mathit{h}}_{\mathit{i}\mathit{j}}$,${\mathit{v}}_{\mathit{i}\mathit{j}}$)
id	x	y	z	hi1	vi1	hi2	vi2	hi3	vi3	RMS
11	−171.283	−118.827	−4.814	−745.019	400.151	−316.781	−623.685	362.860	679.588	0.334
81	123.935	−80.850	−0.931	446.849	247.991	−221.605	392.106	243.357	−423.092	0.085
150	−389.999	271.547	62.835	−1370.495	−757.428	673.595	−1141.890	−734.434	1274.263	0.691
194	254.725	−340.006	−97.490	1086.097	1369.686	−1175.381	899.456	1284.152	−994.748	1.617
942	48.080	393.260	65.336	84.859	−1027.708	869.982	127.453	−958.663	−95.253	1.144
943	37.041	366.774	50.754	52.101	−940.542	795.293	96.428	−877.274	−63.611	0.970
944	73.083	376.546	50.491	165.633	−959.464	808.956	195.930	−892.928	−170.040	1.008
947	17.904	395.752	53.643	−10.528	−998.670	847.403	43.948	−934.252	−4.878	1.039
948	2.228	387.889	20.125	−57.610	−895.363	759.110	0.405	−839.240	42.036	0.878

, along with corresponding code identification and assignation (id) for coded and non-coded targets, respectively. The reprojection error minimization performance is also shown for computing each individual target, with a RMS value ranging from 0.085 pixels to 1.617 pixels. Again, the mean computing time ranges below 0.1 ms per iteration.

4.3. Joint Bundle Adjustment

Once the new target coordinates are solved, other than those defining the measuring frame, the bundle adjustment for joint computation ($\theta$ in Equation (7)) of i = 1 .. N target 3D coordinates (${\theta }_{x_i}$) and j = 1 ..M camera extrinsic parameters (${\theta }_{e_j}$) proceeds, given by their interdependency through the epipolar net multiple view geometry (Figure 6c). Now, partial derivatives (Equation (9)) in a Jacobian matrix $J$ can be defined (Equation (19)) for minimizing a joint reprojection error vector of 2 m elements, being m as defined at the end of Section 3.

$$J\left(\theta \right)=\left[E X\right]$$

(19)

with ${\left(E\right)}_{2m\times 6M}$ and ${\left(X\right)}_{2m\times 3N}$ containing the partial derivatives of each reprojection error with respect to all camera extrinsic (α_j, β_j, γ_j, and d_j) and target coordinates (X_i), respectively. Reprojection error vector elements can be arranged so that $E$ can be defined as a non-square diagonal matrix with ${\left(E_j\right)}_{2N_j\times 6}$ submatrices in its diagonal with the partial derivative to each image extrinsic parameters ${\theta }_{e_j}$ (Equation (11)) per the set of N_j reprojection errors of all targets detected in each j-th image, and $X$ can be correspondingly expressed as a sparse matrix with ${\left(X_{ij}\right)}_{2\times 3}$ submatrices with partial derivative to each target coordinate ${\theta }_{x_i}$ (Equation (17)) per the set of M_i reprojection errors of each i^th target detected in its subset of images per $X$ column and zeros in the rest of elements. Given this arrangement, $J^{T}J$ in Equation (8) can be expressed as follows

$${\left[E X\right]}^T\left[E X\right]=\left[\begin{array}{cc}A& B\\ B^T& C\end{array}\right]$$

(20)

where ${\left(A\right)}_{6M\times 6M}$ is a square diagonal matrix in which its diagonal is composed of ${(E_j{}^T{E}_j)}_{6\times 6}$ square submatrices accounting for each j-th camera extrinsic contribution to the minimization problem in the least square sense, given $E_j$ as defined in Equation (11), and ${\left(C\right)}_{3N\times 3N}$ is also a square diagonal matrix in which its diagonal is composed of ${({\displaystyle \sum }_{j=1}^{M_i}{X}_{ij}^T{X}_{ij})}_{3\times 3}$ square submatrices accounting for the contribution of the 3D coordinate of each corresponding i-th target detected in its M_i image subtset, given $X_{ij}$ as defined in Equation (17).

${\left(B\right)}_{6M\times 3N}$ is a non-square sparse matrix where the interdependency between extrinsic parameters and target coordinate computation through the joint epipolar net is taken into account, with elements ${(E_{ij}^T{X}_{ij})}_{6\times 3}$ contributing along with the joint product of the partial derivatives to image extrinsic and target coordinates for a i-th target detected in a j-th image, and zero when a i-th target is not seen in a j-th image, given $E_{ij}$ as defined in Equation (12).

Thus, according to Equation (20), Equation (8) can be decomposed as

$$\left[\begin{array}{cc}A& B\\ B^T& C\end{array}\right]\left[\begin{array}{c}{∆}_{{\theta }_e}\\ {∆}_{{\theta }_x}\end{array}\right]=-J^{T}r$$

(21)

where the residual term $J^{T}r$ can be redefined according to Equation (19) as

$$\left[\begin{array}{c}{ϵ}_e\\ {ϵ}_x\end{array}\right]=-J^{T}r=-\left[\begin{array}{c}{E}^T\\ X^T\end{array}\right]r$$

(22)

Given Equations (21) and (22), the linear system in Equation (8) can now be decomposed to a set of two lower range subsystems where ${∆}_{{\theta }_e}$ and ${∆}_{{\theta }_x}$ can be individually solved for interdependent extrinsic and target coordinate iteration, where extrinsic iteration ${∆}_{{\theta }_e}$ can be obtained as

$${∆}_{{\theta }_e}={\left(A-B{C}^{-1}{B}^T\right)}^{-1}\left({ϵ}_e-B{C}^{-1} {ϵ}_x\right)$$

(23)

and target coordinate iteration ${∆}_{{\theta }_x}$ can correspondingly be given as

$${∆}_{{\theta }_x}=C^{-1}\left({ϵ}_x-B^T {∆}_{{\theta }_e}\right)$$

(24)

Following the examples in Section 4.1 and Section 4.2, given the computed initial approaches for the three first image extrinsic parameters (Table 2) and for the corresponding target coordinates (Table 3) in the evaluation scene shown in Figure 9b, their joint computation can be conducted according to Equations (23) and (24). As inputs for the method, 3D coordinates of the set of coded targets on the measuring frame are known (Table 1), as well as nominal camera intrinsic parameters, along with image coordinates of coded and non-coded target with solved correspondences between images (Table 1 and Table 3).

Table 4. Camera extrinsic parameters ($d_x$, $d_y$, and $d_z$ in mm, and $\alpha$, $\beta$, and $\gamma$ in radians) after joint bundle of epipolar net scene at Figure 9b.

	dX	dY	dZ	α	β	γ
Image 1	−13.561	5.492	1145.880	−2.378	0.005	0.020
Image 2	−6.635	−7.598	1339.838	−2.347	−0.008	−1.580
Image 3	7.005	10.447	1232.964	−2.389	0.023	−4.712

and

Table 5. Coded (id) target 3D coordinates (mm) after joint bundle of scene at Figure 9b.

id	x	y	z
11	−170.704	−121.223	−1.921
81	123.995	−80.446	−1.397
150	−390.866	274.147	60.943
194	254.332	−340.750	−96.037
942	49.131	422.060	46.501
943	37.779	392.934	32.957
944	74.603	403.023	32.620
947	18.272	425.879	33.690
948	2.222	417.722	−0.569

show the results of the intermediate joint bundle computing for extrinsic parameters and targets (step 6 of the in-process procedure), respectively. The RMS of the joint reprojection error is optimized to 0.802 pixel, with an iteration time of 0.17 ms.

The measuring process may now continue with the acquisition and in-process computing of new images, following the in-process procedure described in the introduction of Section 4 (steps 1 to 6). An example of a complete measuring set is shown in Figure 10a for the pilot case under study (Figure 3), where non-coded targets were placed on the milled surfaces and coded targets were placed around the scene enabling a consistent epipolar net construction during measurement, along with the reference frame with its corresponding coded targets. A total number of 68 images were taken, solving the 3D coordinates of 80 non-coded (20 at each milled prism) and 340 coded targets. A maximum joint bundle computation time ranging 150 ms was observed for the last images of the measuring process (see Section 4.4).

So far, the measuring scale and traceability depend on the coded target 3D coordinates defining the reference frame (Table 1), along with the adopted nominal camera intrinsic parameters. An alternative for increasing accuracy in portable photogrammetry is the adoption of appropriate scale bars with calibrated lengths between corresponding pairs of optical targets (highlighted in orange in Figure 10a), so that measuring process traceability is given by precise scale bar distances. Then, once image-taking finishes (steps 1 to 6), the final post-process joint bundle of camera extrinsic parameters and target coordinates can proceed until convergence (step 7), imposing calibrated relative distances between corresponding pairs of target coordinates in each scale bar available in the scene. A free network adjustment is then carried out. For this, the 3D coordinates of the coded targets on the reference frame (Table 1) have to be also computed, so that their assumed coordinates so far (in steps 1 to 6) do not influence measuring process traceability, and the measuring frame is correspondingly redefined according to their computed coordinates following the widely adopted 3 (to define a reference plane) −2 (to define a reference axis) −1 (to define the origin point) rule in metrology, setting back a determined origin and orientation to the measuring coordinate system.

To do so, a corresponding error vector $r_t$ can be expressed (Equation (25)) and its corresponding joint minimization can be accomplished along with the reprojection error vector $r$ in Equation (4).

$${\left(r_t\right)}_{\left(B+6\right)\times 1}=\left[\begin{array}{c}{r}_b\\ r_x\end{array}\right]$$

(25)

where $r_b$ expresses (Equation (26)) the square distance between each pair of target coordinates minus the corresponding scale square length, given by ${\left|x_k^b-x_k^e\right|}^2-b_k^2$, from k = 1 ..B bars with $b_k^{}$ scales in the scene, with $x_k^b$ and $x_k^e$ being the 3D computed coordinates for each pair of targets of the k^th scale bar, which can be expressed as

$${\left(r_b\right)}_{B\times 1}=\left[\begin{array}{c}\begin{array}{c}{\left|x_1^b-x_1^e\right|}^2-b_1^2\\ ⋮\end{array}\\ \begin{array}{c}{\left|x_k^b-x_k^e\right|}^2-b_k^2\\ ⋮\end{array}\\ {\left|x_B^b-x_B^e\right|}^2-b_B^2\end{array}\right]$$

(26)

and $r_x$ contains (Equation (27)) the 3D coordinate errors of a set of reference targets selected for determining the measuring coordinate system according to the 3-2-1 rule, which can be expressed as

$${\left(r_x\right)}_{6\times 1}={\left[x_0 y_0 z_0 y_1 z_1 z_2\right]}^T$$

(27)

where a first reference target X₀ is set and constrained to be the coordinate system origin (X₁ in Table 1), so that $X_0={\left[x_0 y_0 z_0\right]}^T$ coordinates have to be minimized to zero as first elements in $r_x$ setting 3 restrictions, a second reference target X₁ is constrained (X₃ in Table 1) to have its Y and Z coordinates to zero, so that $y_1$ and $z_1$ coordinates in $X_1={\left[x_1 y_1 z_1\right]}^T$ are included as the next elements in $r_x$ setting 2 restrictions determining X coordinate axis, and a last reference target X₂ is constrained (X₄ in Table 1) to have its Z coordinate to zero, so that $z_2$ coordinate in $X_2={\left[x_2 y_2 z_2\right]}^T$ is included as the last elements in $r_x$ setting 1 restriction determining XY plane of the measuring frame. As a result, 6 constraints are included and a coordinate system is determined, the rest of the initially assumed reference target coordinates (Table 1) being unconstrained and correspondingly computed, so that measuring traceability is now set by the scales imposed in Equation (26).

A Jacobian matrix $G$ corresponding to $r_t$ minimization problem (as in Equation (8)) can be defined as follows

$${\left(G\right)}_{\left(B+6\right)\times \left(3N\right)}=\left[\begin{array}{c}{G}_b\\ G_x\end{array}\right]$$

(28)

where ${\left(G_b\right)}_{B\times \left(3N\right)}$ contains (Equation (29)) the partial derivatives of the k = 1 ..B square distance errors to the each pair of corresponding target coordinates $x_k^b$ and $x_k^e$ given as

$${\left(G_b\right)}_{B\times \left(3N\right)}=\left[\begin{array}{c}{G}_1\\ ⋮\\ \begin{array}{c}\begin{array}{c}{G}_k\\ ⋮\end{array}\\ G_B\end{array}\end{array}\right]$$

(29)

being each $G_k$

$${\left(G_k\right)}_{1\times \left(3N\right)}=\left[\cdots \cdots \cdots D{x}_k^b \cdots \cdots \cdots D{x}_k^e \cdots \cdots \cdots \right]$$

(30)

with derivatives $D{x}_k^b$ and $D{x}_k^e$ following Equations (31) and (32) and zeros for the rest of the elements in $G_k$, being N the new total number of computed targets,

$$D{x}_k^b=2\left[\left(x_k^b-x_k^e\right) \left(y_k^b-y_k^e\right) \left(z_k^b-z_k^e\right)\right]$$

(31)

$$D{x}_k^e=-D{x}_k^b$$

(32)

and ${\left(G_x\right)}_{6\times \left(3N\right)}$ expresses the partial derivatives of $r_x$ to the constrained reference target coordinates, equal to 1 at the columns corresponding to each computed coordinate, and zeros in the rest of the elements.

Thus, a joint error vector $r_s$ can be defined as follows (Equation (33)) integrating reprojection errors in $r$ along with error vector $r_t$

$$r_s=\left[\begin{array}{c}r\\ r_t\end{array}\right]$$

(33)

and joint bundle can be conducted (as in Equation (8)) with the corresponding redefined joint Jacobian for numerical iteration as

$${\left(J\right)}_{\left(2m+B+6\right)\times \left(6M+3N\right)}=\left[\begin{array}{cc}E& X\\ 0& G\end{array}\right]$$

(34)

Following the same system decomposition approach as in Equations (21)–(24), it can be demonstrated that the main iteration equations given by Equation (21) can be redefined to a similar form as

$$\left[\begin{array}{cc}A& B\\ B^T& C_t\end{array}\right]\left[\begin{array}{c}{∆}_{{\theta }_e}\\ {∆}_{{\theta }_x}\end{array}\right]=\left[\begin{array}{c}{ϵ}_e\\ {ϵ}_{x_t}\end{array}\right]$$

(35)

where terms $C$ and ${ϵ}_x$ are redefined to $C_t$ and ${ϵ}_{x_t}$ according to $G$ and $r_t$ as

$$\begin{gathered} C_t=C+G^{T}G \\ {ϵ}_{x_t}={ϵ}_x-G^T{r}_t \end{gathered}$$

(36)

so that expressions given at Equations (23) and (24) for ${∆}_{{\theta }_e}$ and ${∆}_{{\theta }_x}$, respectively, can be used for computing to convergence the post-process joint bundle (step 7) of the scene extrinsic and target coordinates integrating scale and 3-2-1 rule restrictions. As input for the post-process joint bundle, results after processing the last image at step 6 are known. Figure 10b shows the corresponding histogram of the minimized reprojection error vector $r$ after joint bundle with $r_t$, resulting in a RMS value of 2.378 pixel.

4.4. Computing Performance of the In-Process Approach

The procedure described above was developed in C++ language, using the OpenCV library for image processing and the Eigen library for matrix management and processing, on a desktop PC Intel Core i7-5600U 2.6 Ghz, with 16 Gb RAM, running on Windows 7 with 64 bits. Wireless image transmission was used, observing a transmission time of 0.5–1 s per image for 12.2 Mpixel raw images. Image processing time (step 1) was observed to reach 0.5 s per image, depending on the number of segmented and decoded targets on the image. Computing time for the initial extrinsic and target approach (steps 2 and 5) for new images and targets prior to their first bundle could be relatively neglected, ranging below 1 ms. Code assignation to non-coded targets (step 4) was observed to range up to 0.1–0.3 s per image, depending on the number of non-coded targets and correspondences to solve between images.

Figure 11 shows the dependence of the computation time of the intermediate bundle (step 6) on the number of cumulated images and solved targets so far during the measuring process, with a time ranging up to 0.15 s per image for the last images of the measuring scenario shown in Figure 10. As a result, total in-process computation time (steps 1 to 6) ranged at a maximum of 1 s per image, in the same order of magnitude of the wireless image transmission time itself.

However, an average increase of 3 ms per additional image can be estimated from Figure 11 for step 6, which would lead to relevant in-process computation time contribution if larger measuring scenarios were adopted with more images and targets. If a better in-process computing performance was required for step 6, further optimization could be conducted following the system decomposition approach presented in this paper, taking advantage of the symmetry and characteristics of the involved submatrices in Equation (21) (diagonal A and C submatrices, dominating presence of zeros in sparse B submatrix, etc.) along with the development of analytic expressions for the involved inverse matrices in Equations (23) and (24).

Finally, the post-process joint computation time (step 7) reached 3 s in the same scenario, with 5 iterations to convergence, assuming convergence criteria of ${∆}_{{\theta }_e}$ and ${∆}_{{\theta }_x}$ being below a minimum iteration value of 10⁻⁶ mm and 10⁻⁶ radians for all location and orientations, respectively.

5. Camera Model Self-Calibration for Precision

Along with scale bars in the scene, camera intrinsic parameters determine measuring process traceability. Nominal camera intrinsic parameters have been adopted for the results shown so far. As introduced in Section 2, a length measuring error (LME) approach similar to VDI 2634 can be adopted for evaluating measuring procedure uncertainty in portable photogrammetry. Figure 12 shows the same measuring scenario as in Figure 10a, but additional six scale bars are included covering the scene so that just one bar is used for imposing scale to the measurement and LME errors can be evaluated on the rest of measured scales.

A maximum LME error of 1.2 mm was observed in the photogrammetric measurement given calibrated scale lengths ranging 1340 mm, corresponding to a limited relative precision of approximately 1/1000 similar to the reported in [5], one order of magnitude below the expected performance of portable photogrammetry using precise pre-calibrated cameras, where a relative precision ranging better than 1/10,000 is typically reported [20] for 1 m long scenes.

Other than assuming nominal values for intrinsic parameters, precise camera calibration is required for precise measuring. Accurate camera calibration is necessary in photogrammetry systems and that is why continuous and recent improvements regarding mathematical analysis of the data taken for self-calibration [11,12,13] and image distortion correction [14,15] can be found in the literature. Nowadays there are several software applications for close range photogrammetry in industry using optical targets (i.e., Aicon 3D, Geodetic V-Stars, Creaform MaxShot3D, GOM Tritop, etc.) that can automatically perform camera self-calibration. These photogrammetry specific commercial systems present a main limitation when talking about their self-calibration since they only provide this option by a specific procedure that must be performed out of the measuring process itself since traditional self-calibration methods need precise calibration pattern and extensive post-processing work. This does not enable the option of in-process self-calibration. Such self-calibration capabilities integrated into the measuring process would be essential to compensate for main uncertainty contributors, especially by determining the intrinsic parameters of the camera model on the fly. In addition, such a self-calibrated machine vision would enable the use of low-cost cameras—and not only specific photogrammetry cameras—for metrology applications.

Some approaches for the calibration of commercial digital cameras in different applications can be found in the literature. In [32] an accuracy comparison is carried out between three low-cost consumer grade digital cameras and a specific photogrammetry proven camera. In [33] a study of different commercial software options for the calibration and self-calibration of cameras mounted on unmanned aerial vehicles (UAVs) is presented. A similar approach for the self-calibration of commercial small action cameras applied to photogrammetric purpose is proposed in [34] by using a 2D external calibrated reference and a specifically developed software. However, the options presented do not allow in-process self-calibration since they need extensive post processing.

To overcome the time limitations of off-process calibrations, some approaches propose diagnosis methods for camera internal parameters in order to prevent the need to stop the measuring process so often to check them: in [35] a diagnostic method for internal parameters based on multivariate control charts is proposed in order to provide a comprehensive stability control over all the performed calibrations for systems used for regular monitoring of production lines. However, this cannot be considered a self-calibration, but a stability control of the intrinsic parameters.

Special attention should be paid to texture-based software for 3D reconstruction and modelling (i.e., Agisoft PhotoScan, Photomodeler Scanner, Bundler, Pix4Dmapper, VisualSFM, iWitness, MicMac, 3DF Zephir, etc.), with a broad range of applications (e.g., mapping, architecture, mining, construction, agriculture, heritage recording, archeology, forensic 3D modelling, etc.). Other than using optical retroreflective targets, measuring results are obtained on object surface textures available in the images. Texture-based techniques (SIFT [36] and SURF [37] image descriptors, SfM [38,39,40], etc.) are used, enabling dense 3D surface meshing and reconstruction. In-process camera model self-calibration techniques are integrated, demonstrating the capability of 3D surface reconstruction and modelling using consumer grade low-cost cameras [41,42,43]. However, strong limitations can be found for taking these texture-based approaches to precise industrial metrology. Textures available in the images depend on object characteristics and light conditions, limiting measurement precision and reliability compared to the use of physical optical targets. Additionally, computationally hard texture-based approaches are required, with, again, extensive post-processing times, limiting its application for efficient measuring processes in industry by photogrammetry.

The in-process self-calibration method proposed in this work overcomes these limitations by allowing its computationally efficient application for precise industrial metrology applications using optical targets. The method for the self-calibration of camera and lens distortion has been integrated due to its potential for highest precision and accuracy level [20,44] when using low cost non specialized digital cameras, where the solved 3D point cloud scene itself is used as calibration geometry. Measuring traceability is set only by the scale bars in the measuring scene, avoiding uncertainty contributors from off-process camera calibration processes. The camera can be continuously controlled and compensated on the fly, bringing the potential of getting independent measurement results to changes in camera condition over time, such as with thermo-mechanically unstable low-cost cameras. As a step forward to off-process camera calibration and in-process stability control, the approach adopted in this work consist in taking advantage of redundant information available in the portable photogrammetry scene, so that, along with extrinsic and target coordinates, camera intrinsic parameters are also included into a final bundle adjustment computing, having as the only inputs the measuring images themselves and the scale bar distances for precise traceability.

5.1. Including Camera Model into Bundle Adjustment

Computation of the camera intrinsic parameters means determining vector ${\theta }_i$ (Equation (37)) along with camera extrinsic ${\theta }_e$ (Equation (5)) and target coordinates in ${\theta }_x$ (Equation (6)) in a joint bundle minimizing error vector $r_s$ (Equation (33)) including scale bars restrictions.

$${\left({\theta }_i\right)}_{7\times 1}={\left[f h_0 v_0 k_1 k_2 {\pi }_1 {\pi }_2 \right]}^T$$

(37)

where the focal distance $f$, the coordinates of the distortion center ($h_0$,$v_0$), the radial distortion coefficients ($k_1$ and $k_2$) and the tangential distortion coefficients (${\pi }_1$ and ${\pi }_2$) are included.

Hence, $\theta$ vector to compute in the joint bundle is redefined as

$$\theta =\left[\begin{array}{c}{\theta }_e\\ {\theta }_x\\ {\theta }_i\end{array}\right]$$

(38)

and following Equation (34), a corresponding complete joint Jacobian matrix can be expressed as

$${\left(J\right)}_{\left(2m+B+6\right)\times \left(6M+3N+7\right)}=\left[\begin{array}{ccc}E& X& F\\ 0& G& 0\end{array}\right]$$

(39)

where $F$ expresses the partial derivatives of the reprojection error vector $r$ to the camera intrinsic parameters in ${\theta }_i$, so that

$${\left(F\right)}_{2m\times 7}=\left[\begin{array}{c}{F}_{11}\\ ⋮\\ \begin{array}{c}\begin{array}{c}{F}_{ij}\\ ⋮\end{array}\\ F_{2m}\end{array}\end{array}\right]$$

(40)

being $F_{ij}$ the partial derivatives of the reprojection error $r_{ij}$ of the i-th target in the j-th camera view, according to the projection model in Equation (2) and the distortion model in Equation (3), as

$$F_{ij}=\left[\begin{array}{ccccccc}\frac{\partial r_{ij}}{\partial f}& \frac{\partial r_{ij}}{\partial h_0}& \frac{\partial r_{ij}}{\partial v_0}& \frac{\partial r_{ij}}{\partial k_1}& \frac{\partial r_{ij}}{\partial k_2}& \frac{\partial r_{ij}}{\partial {\pi }_1}& \frac{\partial r_{ij}}{\partial {\pi }_2}\end{array}\right]$$

(41)

Given the partial derivatives in $F$ and following again the same system decomposition approach as in Equation (35) for the join bundle including scale bars, the main iteration equations including camera self-calibration can be redefined as

$$\left[\begin{array}{ccc}A& B& E_f\\ B^T& C_t& D_f\\ E_f^T& D_f^T& F_f\end{array}\right]\left[\begin{array}{c}{∆}_{{\theta }_e}\\ {∆}_{{\theta }_x}\\ {∆}_{{\theta }_f}\end{array}\right]=\left[\begin{array}{c}{ϵ}_e\\ {ϵ}_{x_t}\\ {ϵ}_f\end{array}\right]$$

(42)

where

$$F_f=F^{T}F$$

(43)

$$E_f=E^{T}F$$

(44)

$$D_f=X^{T}F$$

(45)

being $F$, $E$ and $X$ defined as in Equations (40) and (19), respectively, where ${\left(F_f\right)}_{7\times 7}$ is a square symmetric matrix accounting for each camera intrinsic parameter direct contribution to the minimization problem, and ${\left(E_f\right)}_{6M\times 7}$ and ${\left(D_f\right)}_{3N\times 7}$ are non-square matrices accounting for the corresponding joint product of partial derivatives integrating the interdependency between extrinsic and target coordinate computing with camera intrinsic parameters through the joint epipolar net, respectively, and ${ϵ}_f$ is expressed as

$${ϵ}_f=F^{T}r$$

(46)

analog to ${ϵ}_e$ and ${ϵ}_x$ in Equation (22).

From the linear system in Equation (42), it can be obtained

$${∆}_{{\theta }_f}=F_f{}^{-1}\left({ϵ}_f-E_f{}^T{∆}_{{\theta }_e}-D_f{}^T{∆}_{{\theta }_x}\right)$$

(47)

and then, Equation (42) can be reduced to a similar system as Equation (35)

$$\left[\begin{array}{cc}{A}_f& B_{f1}\\ B_{f2}& C_f\end{array}\right]\left[\begin{array}{c}{∆}_{{\theta }_e}\\ {∆}_{{\theta }_x}\end{array}\right]=\left[\begin{array}{c}{ϵ}_{e_f}\\ {ϵ}_{x_f}\end{array}\right]$$

(48)

where ${∆}_{{\theta }_e}$ and ${∆}_{{\theta }_x}$ can be obtained following analog equations to Equations (23) and (24), substituting $A$, $B$, $B^T$, $C$, ${ϵ}_e$ and ${ϵ}_x$, by $A_f$, $B_{f1}$, $B_{f2}$, $C_f$, ${ϵ}_{e_f}$ and ${ϵ}_{x_f}$, given them as

$$A_f=A-E_f{F}_f{}^{-1}{E}_f{}^T$$

(49)

$$B_{f1}=B-E_f{F}_f{}^{-1}{D}_f{}^T$$

(50)

$$B_{f2}=B^T-D_f{F}_f{}^{-1}{E}_f{}^T$$

(51)

$$C_f=C_t-D_f{F}_f{}^{-1}{D}_f{}^T$$

(52)

$${ϵ}_{e_f}={ϵ}_e-E_f{F}_f{}^{-1}{ϵ}_f$$

(53)

$${ϵ}_{x_f}={ϵ}_{x_t}-D_f{F}_f{}^{-1}{ϵ}_f$$

(54)

5.2. Computing Efficiency and Precision Performance Evaluation for Self-Calibrated Photogrammetry

Bundle computing that includes self-calibration can be conducted following Equations (47) and (48) for the joint iteration of ${\theta }_e$, ${\theta }_x$, and ${\theta }_i$ to convergence, having as initial approach the previous joint bundle computation given constant camera intrinsic parameters. Following the examples in Figure 10 and Figure 12, post-process joint bundle including self-calibration was conducted including the same scale bar (L3) for traceability.

Figure 13 shows the histogram of the optimized reprojection error distribution, where an order of magnitude lower RMS value is obtained compared to that obtained by assuming nominal camera instrinsics (histogram in Figure 10b).

Table 6. Camera intrinsic parameters after self-calibration.

f (mm)	$\mathit{c}{\mathit{l}}_{\mathit{0}}$ (pixel)	$\mathit{r}{\mathit{w}}_{\mathit{0}}$ (pixel)	${\mathit{k}}_{\mathit{1}}$ (pixel−2)	${\mathit{k}}_{\mathbf{2}}$ (pixel−4)	${\mathit{\pi }}_{\mathit{1}}$ (pixel−1)	${\mathit{\pi }}_{\mathbf{2}}$ (pixel−1)
24.557	−7.652	31.396	4.866 × 10−09	−2.150 × 10−16	−8.636 × 10−08	−1.236 × 10−08

shows the computed camera intrinsic parameters for this example, optimized to those assumed to be nominal in Section 4. Post-process joint computation including self-calibration took up to 8 s, with 11 iterations to convergence. As a result, the total post-processing after measuring was limited to 8 + 3 = 11 s, enabling on the fly camera self-calibration and measuring self-compensation, assuming a constant camera condition during measurement.

Again, the length measuring error (LME) approach was adopted for evaluating measuring uncertainty in self-calibrated portable photogrammetry. A set of 10 consecutive photogrammetric measurements were taken on the scene. Each measurement was solved imposing L3 scale bar (Figure 12), and remaining scales were used as control bars for LME evaluation. Maximum LME between targets on control scale bars obtained was 121.4 μm in L0, with a standard deviation for all the set of measurements of 45.1 μm. However, except for L0, the values obtained showed relatively lower LME results for lengths measured in the XY plane (ranging below 40 μm). Nevertheless, that maximum observed LME was assumed as a conservative estimation to characterize optical target measurement process uncertainty in all spatial directions, corresponding to a relative precision ranging 1/10,000.

According to the maximum observed LME value of 121.4 μm, the spatial uncertainty for the target coordinates can be estimated in 70.1 μm ($\sigma$), given that $\sigma =LME$/$\surd 3$ assuming a rectangular distribution. As an alternative way of evaluating measuring performance, the 3D coordinates of the computed non-coded targets in each measurement, properly placed on relatively precise milled prismatic elements of the pilot case (see Figure 3), can be fitted and compared with its nominal geometry [5], and the corresponding errors to each nominal surface can be evaluated for all targets in all 10 measurements. Figure 14 shows an example of fitting of a measurement result to the nominal reference geometry for error evaluation. A maximum standard deviation ($\sigma$) of 14 μm was observed for the measuring errors of non-coded targets to all surfaces in the XY plane in all measurements, but up to 60 μm for the measuring errors of the targets laying in lateral surfaces, close to the estimated $\sigma$ of 70.1 μm given by the LME evaluation. This result points out again the conservative character of the spatial uncertainty ($\sigma$) value previously estimated by the LME approach, due to the observed anisotropy on the non-coded target measurement accuracy. In any case, a relative precision ranging 1/10,000 by LME analysis is confirmed in the worst-case scenario for all spatial directions, demonstrating the adequate accuracy of the developed in-process self-calibrated photogrammetry in the pilot case scenario, comparing to previously reported limited performance of 1/1000 relative precision [5], ranging now at the same precision level as the performance expected when using precise pre-process camera calibration.

6. Evaluation at Industrial Scenarios

Finally, the self-calibrated in-process photogrammetry system has been evaluated in an industrial scenario with two different parts (Figure 15a,c) demonstrating a fast, reliable and precise raw part measuring process. The system was applied for the measurement and alignment of up to 15 m long raw parts prior to their machining, at least one order of magnitude larger than the considered pilot case (Figure 3), with scene sizes ranging now up to 200 m³.

The measurement process took an overall time of 1 h per part in both industrial scenarios, determining the 3D coordinates of the set of non-coded target placed on the surfaces to be machined. An approximated total number of 800 targets and 200 images were computed in both scenarios (Figure 16a,c). Performance of the in-process computing procedure was observed, ranging now at a total maximum in-process computing time of almost 2 s per image for the last images taken during the measuring processes, higher than the 1 s per image observed at the pilot case scenario, but still enabling practical quasi real time diagnosis and control of correct images taken in an industrial scenario, so that every time an incorrect images was acquired (due to an image not contributing to a consistent epipolar net, pending target to be solved in a particular zone of the scene, etc.) measuring could be properly guided to locally take new adequate images if necessary before continuing with the process.

The higher intermediate bundle computing time was the main contributor to the increase of the in-process computation time in the last stages of the measuring process, as expected by the dependence of the bundle computational work on the number of images and targets to solve (as shown in Figure 11). Along with the alternatives pointed out in Section 4.4 for increasing its computing performance, a strategy could also be adopted for limiting bundle adjustment only to a minimum subset of the scene measured so far, so that continuous computing of the joint bundle of the complete scene could be avoided or only periodically executed.

Accordingly, the post-process computing time including self-calibration of camera intrinsic parameters ranged a maximum of 15 s, slightly higher than the observed one in the pilot case scenario as a result of the higher number of images and targets, but meeting efficient operation at an industrial scenario.

Regarding precision, according to the LME evaluation results reported in the literature [20], typical scale dependent LME errors could be estimated as 50 μm + 20 μm/m for portable photogrammetry with precise pre-calibrated cameras. For the considered industrial cases, with up to 15 m long parts, an estimated LME performance of 0.35 mm could be expected as a reference for precise operation, corresponding to 0.20 mm ($\sigma$) spatial uncertainty for the measured targets, given that $\sigma =LME$/$\surd 3$ assuming again a rectangular distribution.

After raw part measurement by photogrammetry, the complete process (Figure 1) was conducted towards in-machine raw part alignment in both industrial scenarios, so that gauging on the optical targets by a touch probe integrated in the machine (with machine axes typically ranging at 0.01 mm accuracy) was adopted as a reference for evaluating measuring precision. Optimal target coordinates were computed by fitting (Figure 17). Positive and even overstock was applied as in [5] as best-fitting criteria. A subset of targets was used as control points for alignment by a machine-integrated contact probe. In-machine fixturing of the part was adjusted to properly align it to machine axes so that probed relative coordinates between reference targets matched those optimally computed by fitting, to a difference below 0.1 mm so that alignment process uncertainty could be considered relatively neglected comparing to the measuring uncertainty of 0.2 mm ($\sigma$) estimated for the portable photogrammetry.

Once the part was precisely aligned, the measuring process performance was evaluated by contact probe gauging of a minimum set of 10 check points (i.e., non-coded targets not used as a reference for alignment), distributed in all 3 machine coordinate axes directions and at extreme and opposite surfaces of the part. A maximum probing error interval of ±0.6 mm (3$\sigma$) could be expected according to the above estimated target spatial uncertainty by LME ($\sigma$ = 0.20 mm). Probing errors were evaluated between the overstock values resulting from the fitting process and the actual ones by the in-machine gauging to nominal surface coordinates. In both industrial scenarios (Figure 15a,c), all the probing errors ranged below ±0.5 mm, demonstrating the adequate accuracy of the developed self-calibrated in-process photogrammetry for large raw part measuring and overstock control.

7. Conclusions

A new efficient procedure has been presented in this work for solving the bundle adjustment problem in portable photogrammetry. In-process bundle computing capability is demonstrated on a consumer grade desktop PC, enabling quasi real time 2D image and 3D scene computing and diagnosis so that a reliable measuring procedure can be conducted, avoiding inefficient user-dependent post-process iterative procedures that limit the potential of portable photogrammetry for an easy, low-cost and fast solution for industrial metrology of large components.

A method for the self-calibration of camera and lens distortion has been integrated into the in-process approach due to its potential for highest precision and accuracy levels when using a non-specialized consumer grade digital camera (i.e., Nikon D300S), where the solved 3D point cloud scene itself is used as calibration geometry and measurement traceability is set only by the scale bars in the measuring scene, avoiding the need of off-process pre-calibrated cameras or the use of special purpose calibration artifacts in the scene for precise measurement in an industrial metrology application.

The developed self-calibrated in-process photogrammetry has been evaluated in a pilot case scenario (1.5 m long reference part) and at two large scale industrial scenarios (up to 15 m) for raw part measurement and alignment before machining, showing an in-process computing time typically below 1 s per image to a maximum of 2 s at the last stages of the computed industrial scenes, and a relative precision of 1/10,000 with an error RMS below 0.2 pixel at image plane, ranging at the same precision performance reported for portable photogrammetry with precise off-process pre-calibrated cameras. Efficient camera model in-process self-calibration is also demonstrated, with post-processing times ranging 11 s.

Alternatives for increasing computational in-process efficiency have been pointed out, especially focusing on large scale industrial scenarios where a high number of images and optical targets to compute can be expected. Regarding precision, anisotropy has been observed in the spatial uncertainty distribution of the optical target coordinates, pointing to a better potential for the developed system than the figures reported according to the worst-case scenario. Further steps towards the prediction of ray-net-conditioning-induced uncertainty contributions, such as analytical approaches using error propagation theory, may enable to take the most of that potential, still remaining as a relevant and challenging issue in the state-of-the-art. Additionally, further steps could be conducted for enabling precise measurements with low cost but thermo-mechanically unstable consumer grade digital cameras, where the influence of an image dependent focal length due to integrated autofocus optics might be considered for precise self-calibration.

Finally, in the actual context, the need for “smart” systems integrating software procedures and hardware systems for data acquisition, self-diagnostics, set-up, control of tolerances, etc., is becoming more and more pressing. In addition, rapid advances in computing processing, image measurement and characterization algorithms have allowed photogrammetric systems to integrate with CAD/CAM systems offering high accuracy capabilities, and more and more “real-time” (1/25 s and faster) image sequence acquisition. The photogrammetric system presented here is a step towards meeting these requirements, accomplishing accurate measurements with inexpensive, easy-to-use measuring systems and with optimized procedures that make the whole measuring process much less time consuming.