BUNDLE ADJUSTMENT-BASED STABILITY ANALYSIS METHOD WITH A CASE STUDY OF A DUAL FLUOROSCOPY IMAGING SYSTEM

A fundamental task in photogrammetry is the temporal stability analysis of a camera/imaging-system’s calibration parameters. This is essential to validate the repeatability of the parameters’ estimation, to detect any behavioural changes in the camera/imaging system and to ensure precise photogrammetric products. Many stability analysis methods exist in the photogrammetric literature; each one has different methodological bases, and advantages and disadvantages. This paper presents a simple and rigorous stability analysis method that can be straightforwardly implemented for a single camera or an imaging system with multiple cameras. The basic collinearity model is used to capture differences between two calibration datasets, and to establish the stability analysis methodology. Geometric simulation is used as a tool to derive image and object space scenarios. Experiments were performed on real calibration datasets from a dual fluoroscopy (DF; X-ray-based) imaging system. The calibration data consisted of hundreds of images and thousands of image observations from six temporal points over a two-day period for a precise evaluation of the DF system stability. The stability of the DF system −for a single camera analysis− was found to be within a range of 0.01 to 0.66 mm in terms of 3D coordinates root-mean-square-error (RMSE), and 0.07 to 0.19 mm for dual cameras analysis. It is to the authors’ best knowledge that this work is the first to address the topic of DF stability analysis.


INTRODUCTION
Stability analysis in this research is aimed at evaluating the stability of calibration parameters of an imaging system.Given two sets of calibration parameters (e.g., interior and/or relative orientation parameters, IOPs and ROPs, respectively) from different temporal points, the stability analysis goal is to determine whether or not the two sets are equivalent.From a photogrammetric-science perspective, stability analysis can be achieved through two major steps.First, photogrammetric calibration sessions have to be performed to estimate the IOPs and ROPs of the imaging system at different temporal points.Second, the estimated parameters from the different calibration sessions should be compared and qualitatively evaluated.Regarding the first step, the calibration sessions can be performed in short time-intervals such as hours or days (i.e., short-term stability; e.g., Shortis et al., 2000;Detchev et al., 2017) or over long time-intervals such as weeks or months (i.e., long-term stability; e.g., Läbe and Förstner, 2004;Detchev et al., 2015).The choice of how frequent the calibration sessions should be performed is application specific and is primarily dependent on the confidence in the imaging system stability.In the case of having an imaging system that is subject to functional or physical changes, the calibration sessions must be performed more frequently to capture potential variations among the system parameters.One should carefully address the photogrammetric principles of rigorous calibration (Fraser, 1997(Fraser, , 1984;;Granshaw, 1980) to guarantee the accurate estimation of the calibration parameters at different temporal points.Of course, having accurate parameter sets will facilitate the later task of parameters comparison and stability assessment.Regarding the second step, given two sets of calibration parameters of the same imaging system, the evaluation of parameters similarity can be achieved through a wide variety of methodologies, each of which has its own advantages and disadvantages.The stability analysis on the two parameter sets can be established at the very basic level of comparing the individual parameter values (e.g., Harvey and Shortis, 1995;Shortis et al., 2000) or through more advanced approaches that utilize stability metrics in image space (e.g., Habib et al., 2014) or in object space (e.g., Lichti et al., 2009).
Stability analysis can be achieved by using statistical hypothesis testing to check the similarity between the parameter values (Harvey and Shortis, 1995;Shortis et al., 2000).Such a methodology assumes Gaussian-distributed parameters with no consideration regarding potential biases.Further, on its own, it provides no measure regarding the impact of parameter differences on the output photogrammetric product.More importantly, one can argue that a small variation in a parameter value would be accepted through a statistical test, but this variation could have high impact on the photogrammetric products.
Instead of comparisons at the parameters level, stability analysis can be achieved by image simulation to capture any differences among calibration parameter sets in terms of image space metrics.A primary advantage of these methods is the elimination of excessive parameters comparison as in case of statistical hypothesis testing.The output of such methods is usually reported by a single value, if the value is less than a preset threshold, the parameter sets are assumed equivalent.For instance, Habib et al. (2006) and Habib and Morgan (2005) have reported three stability methods that are based on image simulation.Briefly, in their work, a simulated bundle of light rays is established from one set of calibration parameters and then compared to another simulated bundle using the other set of calibration parameters.The distance differences between the two bundles are estimated in image space and the RMSE of all differences is used as the stability measure.A variety of positional and rotational constraints are established for the bundle simulation in the three methods to serve different photogrammetric applications/geo-referencing techniques, and accordingly the three methods vary with respect to their strictness in measuring stability.These methods were intended to be used for the stability analysis of the IOPs of an individual camera rather than the IOPs-and-ROPs of an imaging system with multiple cameras.In another publication Habib et al. (2014) reported an advanced methodology for the stability analysis of a multi-camera imaging system (where two cameras were considered at a time) in which the stability of IOPs and ROPs was addressed and reported in terms of image space metrics.Note that the methodology in (Habib et al., 2014) was reported as being novel to the field of system stability analysis (i.e., IOPs-and-ROPs stability).
The above-listed stability analysis approaches are either concerned with the stability of IOPs only (individual cameras) or the stability of IOP-and-ROP sets (imaging system with multiple cameras).Nevertheless, due to methodological constraints, even the more advanced method by Habib et al. (2014) can only be applied on one camera pair at-a-time rather than simultaneously evaluating the stability of a complete imaging system with multi-cameras.As these methods are based on image space simulation, they provide no direct assessment of the object space reconstruction quality.In general, the ultimate goal of any photogrammetric application is the reconstruction of object space (e.g., 3D coordinates).Accordingly, the evaluation of stability in terms of object space rather than image space metrics would be more representative.Lichti et al., (2009) have reported a stability analysis methodology that is based on object space simulation.This method focuses on aerial photogrammetry applications.Camera position and orientation constraints are specified to produce simulated stereo image pairs.Image intersection is used to verify the similarity of calibration parameter sets.The stability is assessed with respect to the object space shape and the quality of derived 3D coordinates.While this method did perform a direct object space stability analysis, it was limited to a single camera system.
Here a new simulation-based methodology for stability analysis is introduced.This method uses the collinearity model to evaluate the similarity between calibration parameter sets.The proposed advantages of this method are: (1) The ability of reporting stability measures in both image and object space simultaneously; (2) The ability of performing stability analysis of a single camera parameters (i.e., IOPs only), stereo-camera pair or a multi-cameras imaging system (i.e., IOPs and ROPs); (3) The flexibility to analyze reconstruction scenarios with a stationary camera/imaging-system and portable/moving object, or stationary object and portable camera/imaging system; and (4) The flexibility of comparing different distortion models (i.e., different additional distortion parameter sets).
The paper began with a review of relevant literature regarding the common stability analysis approaches.Section 2 provides a brief background on the need for camera calibration and the requirements for stability analysis.Section 3 explains the proposed stability analysis methodology.The proposed methodology is then evaluated in the results section on DF calibration data (Section 4).Finally, the paper presents relevant conclusions and recommendations for future work (Section 5).

BACKGROUND
Nowadays, the derivation of 3D coordinates from imagery can be seen in engineering, medical, industrial and many other disciplines and applications.Regardless of the type of the camera/imaging-system used and the application of interest, the fundamentals of photogrammetry must be applied for 3D coordinates derivation.For instance, a camera/imaging-system must be parameterized by a set of mathematical parameters to establish the geometric relation between image and objectspace.A classical parametrization of a single-camera internalgeometry is established through interior orientation parameters (IOPs), which ideally consist of a principal point (  ,   ) and a principal distance () in addition to lens/lenses distortion parameters.Also, a classical parametrization of an imaging system with multiple cameras can be established by a set of IOPs for each individual camera and relative orientation parameters (ROPs) between the different cameras.Each ROP set consists of rotational and directional parameters (e.g., three rotation angles and a 3D translation vector).
The IOPs/ROPs of a camera/imaging-system may be subject to changes on their values over time due to changes in the system components (e.g., lens components move due to a consistent vibration applied on a camera).Here comes the role of stability analysis to detect any temporal changes of IOP/ROP values and to quantify the impact of the changes −if they exist− on the quality of the derived 3D coordinates.A stable camera/imagingsystem is the one that maintains equivalent values of parameters and accordingly yields similar reconstruction quality at different temporal points.A non-stable camera/imaging-system is the one that temporally exhibits significant variations in the parameters and provides inconsistent quality of photogrammetric products (Lichti, 2008).Many factors determine the stability status of an imaging system; in general, if an imaging system preserves a fixed geometric relation between its' components and maintains the same functional specifications it will hold the 'stable' status.
If any of these were disturbed for any reason, the system will be deemed non-stable.The decision of the degree to which an imaging system is stable/unstable is primarily application dependent and should be governed by the permitted error tolerance in the 3D reconstruction.

METHODOLOGY
This section introduces the proposed stability analysis methodology using the bundle adjustment (i.e., collinearity) model.The methodology can be applied for calibration parameters of a single camera, an imaging system with stereo cameras (i.e., only one camera pair) or an imaging system with multiple cameras.It is based on the assumption that if two sets of calibration parameters are equivalent they should provide similar quality of photogrammetric products (e.g., 3D object space coordinates) under the same image network geometry.In other words, given two equivalent sets of calibration parameters of the same camera/imaging system (e.g., Set 'A' and Set 'B'), and if these sets were used separately to derive a certain 3D object from equivalent image networks (i.e., similar number of images/image points, with same locations and orientations), both sets should provide the same quality of the reconstructed object.A further assumption is: any interchange at the individual parameter level between 'A' and 'B' should have no impact on the resulting reconstruction quality.Accordingly, any significant differences between 'A' and 'B' will be captured through variations in a derived 3D object using 'A' and 'B'.
In this work, imaging simulation is used to capture the differences between the two parameter sets.A simulated image network is established over a defined/simulated 3D object.A parameter set 'A' is used to generate image measurements/observations and to introduce image distortions.
A parameter set 'B' is used to correct the image coordinates from distortions and to run a least squares bundle adjustment.
Having the 3D reconstruction as the target function of the bundle adjustment, only the 3D object space coordinates are treated as the unknown parameters.The other bundle adjustment parameters are fixed (i.e., introduced as known values) to ensure precise comparison between 'A' and 'B'.The variation between 'A' and 'B' is evaluated in 3D by finding the differences between the estimated and the simulated object space coordinates.The image observation residuals resulting from the adjustment can also be used to quantify the variation between 'A' and 'B' in image space.Assuming no difference between 'A' and 'B', both coordinate sets will be identical, and the residuals values will be zero.By having the IOPs and the exterior orientation parameters (EOPs) fixed through the adjustment, all the differences between parameter sets will propagate as errors in the object space reconstruction.If the EOPs and IOPs are not fixed, part of the discrepancies in the different sets may be compensated.This methodology can be applied for stability analysis of a single camera or an imaging system (multiple cameras), with a slight methodological change for the latter case in order to account for the type of differences in the parameters being evaluated in both scenarios.For a single camera, the target function of stability analysis is to evaluate the differences in IOP sets.For a multi-camera imaging system, the stability analysis is focused on the simultaneous effect of two different IOP and ROP sets.

Stability assessment of a single camera
In the case of a photogrammetric application in which a portable single-camera is used to capture images over a stationary object/scene, and having two sets of IOPs (e.g., 'A' and 'B') to be compared, the stability analysis methodology is described through the following steps alongside equations (1) to ( 16) and Figure 1: (1) Simulation of 3D object/scene that is represented as a set of (  ,   ,   ) vertices.The simulated object can be a simple geometric shape such as a plane, cube, sphere or a more complicated shapes/scene.The shape and dimensions of the simulated object are user defined and should be as close as possible to the expected object products of the photogrammetric application of interest.
(2) Simulation of a network of images over the 3D object/scene.Each image is defined by a unique set of EOPs and the camera IOPs (i.e.,   ,   , ).The image network should be defined in a way that imitates a similar network which is regularly used for the photogrammetric application of interest.It is proposed here that a denser network geometry will yield a more precise stability analysis.Having many images with varying locations and orientations should provide rigorous examination of both parameter sets.
(3) Derivation of image measurements/observations (Figure 1.a).The derivation of (  ,   ) image coordinates is achieved using the collinearity model (equations 1 and 2).Each (  ,   ) point is a function of its (  ,   ,   ) homology, the corresponding image EOPs and the camera IOPs (i.e.,   ,   , ) of 'A'.The simulated image measurements should fall within a pre-defined imageformat/dimensions according to the specifications of the camera of interest.
(4) Addition of distortion to the derived image measurements using the IOPs of 'A' (Equations 4 and 5).
(5) Correction of the image measurements from distortion using the IOPs of 'B' (Equations 15 and 16).
(6) Estimation of the 3D object coordinates (  ′ ,   ′ ,   ′ ) through a least-squares bundle adjustment (Mikhail et al., 2003)    where   and   are the distortion correction terms computed from parameters set 'B'.
For completeness, the image network simulation procedure should carefully address the basics of bundle adjustment.For instance, an object point should appear in at least two images since a single point observation is not sufficient to derive 3D coordinates.

Stability assessment of an imaging system (multicamera)
Given an imaging system of multiple cameras and two sets of calibration parameters 'A' and 'B', with each comprising IOPs and ROPs, a similar methodology to steps 1 to 7 (section 3.1) can be used for system stability analysis.The only required modification in the methodology is to update the image locations and orientations prior to the adjustment (i.e., step 6) in order to capture any variations in the ROPs between 'A' and 'B'.For instance, in the case of a stationary imaging system with two fixed cameras, each camera will be associated with a unique set of EOPs, where the EOPs implicitly define the ROPs between the cameras.Accordingly, for stability analysis, the EOPs from 'A' are used to simulate the image observations (step 2, section 3.1) (Figure 2.a), while the EOPs from 'B' are used to run the adjustment (step 6, section 3.1) (Figure 2.b).
Please note that the stability analysis should adapt to serve different imaging/imaging systems scenarios.In general, there are two scenarios of imaging using multiple cameras in photogrammetry; (1) a stationary imaging-system and a portable or stationary object, or (2) a portable imaging-system and a stationary object.The stability analysis should be performed on an image network that closely imitates the resulting network from the aforementioned scenarios.For the first scenario, and given a two cameras imaging-system, the stability analysis should be performed using one image pair (i.e., one image from each camera) of a simulated 3D object.For the second scenario, the image network may contain many images simulating the movement of the imaging system over the given object.

EXPERIMENTAL RESULTS AND DISCUSSION
The proposed methodology in Section 3 was applied on sets of real calibration parameters for a dual fluoroscopy (DF) imaging system.Section 4.1 briefly introduces the DF system and the calibration data used for stability analysis.Section 4.2 presents the stability analysis experiments on individual cameras.Section 4.3 presents the stability analysis on dual cameras.

DF imaging system and calibration data
The DF system (clinically referred to as biplanar videoradiography; Figure 3) is an X-ray based imaging system used for the six degree-of-freedom human/animal motion analysis.The DF system is used for scientific research at the Movement Assessment Laboratory, University of Calgary, Alberta, Canada.The DF system comprises (1) two X-ray sources, (2) two image intensifiers and (3) two high-speed video cameras.These components are used to provide stereoscopic imaging of a human/animal motion on a treadmill.
Briefly, the DF imaging is performed as follows: (1) an X-ray source produces an X-ray pulse, (2) the X-rays interact with the objects within the imaging field; part of the X-ray energy is attenuated by any radiopaque objects, (3) the image intensifier input screen −at the opposite side from the X-rays source− receives the remaining X-rays energy, which is amplified and converted to optical light rays by the intensifier mechanism, (4) the intensifier output screen provides an optical display of the imaged object, (5) the video camera (attached to the intensifier output screen) records the output screen display.Detailed information about DF systems and the imaging process can be found in (Tashman et al., 2010).
In this research work, the calibration of the DF system parameters (i.e., IOPs and ROPs) is achieved using bundle adjustment with self-calibration.The details on the DF system calibration methodology using bundle adjustment can be found in (Al-Durgham et al., 2017, 2016;Lichti et al., 2015).Briefly, the calibration methodology models a single X-ray source, and a single image intensifier and a single camera, as an ideal pinhole camera model (lumped system model) (Lichti et al., 2015).
Based on this, the DF imaging system is represented by two camera models; camera model 1 and camera model 2.
The IOPs of each camera model consist of 34 parameters (i.e.,   ,   ,  in addition to 31 optical and X-ray distortion parameters).The ROPs between the two cameras consist of six parameters.Six calibration datasets were collected in October 2017 over a two-day period to perform stability analysis on the DF system.The DF system components remained fixed/untouched during the data collection to avoid any disturbance of the calibration parameters.Table 1 lists the calibration data attributes as well as the (, ) RMS of image observation residuals resulting from the bundle adjustments.
The reason for having hundreds of images and thousands of image-observations is to guarantee precise estimation of the DF system calibration parameters (having high redundancy and strong calibration network geometry).Table 2 and Table 3 list sample IOPs from the six epochs (only   ,   ,  are displayed for space management) for camera models 1 and 2, respectively.
As can be seen in Table 2 and Table 3, the   ,   ,  values from the 6 epochs were relatively close to each other with a few pixels differences between epochs.A significant difference (i.e., ⁓14 pixels) can be seen between the values of   of camera model 1 (Epoch1 and 5, Table 2).It is speculated that the image intensifiers are subject to slight changes of internal magnetic field due to changes in electrical current, which might be the reason of such parameter differences.epochs were taken two at a time and tested against each other (i.e., epoch 1 against epoch 2, then epoch 1 against epoch 3, and so on for all epoch pairs).Same procedure was followed for the camera model 2 IOPs.The methodology in section 3.1 was applied to each set of camera parameters.A 20 ×20 (cm 2 ) planar object with dense (, , ) points was simulated for step 1 of the methodology (section 3.1).A planar object was chosen over other geometric shapes to simplify the later task of stability analysis.A dense image network (Figure 4) was simulated around the planar object for step 2 of the methodology (section 3.1).The image network covered the planar object from all directions with having image rotation angles of 0° to 360° for (, , ).The IOP pairs from the six epochs were used for steps 3 to 6 methodology (section 3.1).The distance differences between the simulated object points (step 1) and estimated object points (step 7) were computed, and the RMSE for all the distance differences were then determined.4).This is due to significant differences between the IOPs of epoch 1 and 5 (e.g.,   value of both epochs, Table 2).The lowest RMSE values (i.e., 0.01mm) were between epoch pair (4 and 6) of camera model 1 (Table 4) and epoch pairs (3 and 4) and (3 and 6) of camera model 2 (Table 5).In general, the camera model 2 parameters appeared to be more stable than the ones for camera model 1 judging by the lower RMSE values.

Results of dual-camera stability analysis
The goal of this experiment was to evaluate the stability of the IOPs and ROPs of the DF imaging system.Using the methodology in section 3.2, the IOPs and ROPs from the six epochs were tested against each other (taken two epochs at a time).The experiment was performed by having fixed/stationary cameras (i.e., camera model 1 and camera model 2) defined by the IOPs/ROPs and a simulation of a rotating plane.The plane was rotated through 0 to 360° around its X, Y, Z, axes respectively.The stability analysis was performed using two images (one from each camera) of the plane at different rotation angles.Figure 6.a shows the stability analysis resulting from testing calibration parameters of epoch 1 against calibration parameters of epoch 5.As can be seen in Figure 6.a and Figure 6.b, the variation among RMSE was dependent on the value and type of the plane rotation angle.The object rotation around X and Y resulted in larger/smaller and significantly changing RMSE values when compared to rotation around Z. The rotation of the plane around the Z-axis had little impact on the relative distances between the 3D points from the two cameras (i.e., no change in depth observations); accordingly, the RMSE values were slightly changing through the 360°, which was not the case for X and Y rotations.Based on the stability analysis using image pairs -and given two days calibration data− one can conclude that the changes over IOPs and ROPs of the DF system would result in error of 3D coordinates reconstruction between 0.07 to 0.19 mm.

CONCLUSIONS AND FUTURE WORK
This paper has presented a simple and rigorous stability analysis methodology.The methodology has the flexibility to be applied for a single camera or a multi-camera imaging system.Thanks to the collinearity model, which facilitates the scientific analysis over different imaging scenarios and photogrammetric applications.The case study on a dual fluoroscopy imaging system has reported that the system is subject to calibration parameter instability, which leads to a sub-millimetre coordinates reconstruction error.For a single camera analysis, the stability of the DF system was found to be within a range of 0.01 to 0.66 mm in terms of the 3D coordinates RMSE.For dual cameras analysis the stability of the DF system was found to be within a range 0.07 to 0.19 mm.Future work will focus on stability analysis of calibration parameters from different calibration models (e.g., different IOP parametrizations).Further, the presented methodology will be qualitatively evaluated against other existing stability analysis methods.

Figure 3 .
Figure 3. Dual Fluoroscopy imaging system components; (a) time synchronized X-ray sources (G-1086, Varian, USA), (b) instrumented treadmill (Bertec, USA), (c) 406 mm diameter quad-mode image intensifiers (E5876SD-P2A, Toshiba, Japan), and (d) high-speed video cameras (DIMAX, PCO, Germany) Figure 5.a displays the shape of the estimated plane from epochs 1 and 5, and Figure 5.b displays the shape of estimated plane from epochs 4 and 6 (please note the different scale of the Z axis in both figures).The differences in IOPs of the epoch pairs resulted in larger errors at the image observations near to the image boundaries (image points with large radial distance).Accordingly, the errors in point reconstruction were large at the boundaries of the estimated plane.Different error signs were observed at opposite plane boundaries.

Figure 4 .
Figure 4. Simulated image network (108 images) arround a simulated planar object Based on this analysis, one can conclude that the IOPs of each DF camera model are subject to a few pixels change over time, which may cause a sub-millimetre error in 3D coordinates reconstruction.In real DF reconstruction scenarios, both cameras will be involved to derive 3D coordinates from stereo image pairs (e.g., two images with an 80° convergence angle) Figure 5. Surface display of simulated and estimated object space Figure 6.Stability analyis of stationary imaging-system and rotating plane ) are the pixel coordinates of an object point  that appears in image ; (  ,   ) are the pixel coordinates of the principal point,  is the principal point distance; the opposite signs in the rational terms of the collinearity equations are to account for the reflection required to transform from the lefthanded pixel coordinate system into the right-handed one; ( * ,  * ,  * ) are the object space coordinates of the perspective centre (PC); (  ,   ,   ) are the object space coordinates of the object point ; and (, , ) is the object-to-image rotation matrix that is parametrized as a sequence of Euler angles (, , ).

Table 1 .
Six epochs of calibration data

Results of individual cameras stability analysis
The goal of this experiment was to quantify the impact of IOP variations of each camera model separately.Specifically, the IOPs (i.e., 34 parameters) of camera model 1 from the six

Table 4 and
Table 5 list the RMSE values of all epoch pairs, for camera model 1 and camera model 2, respectively.The largest RMSE value (0.66 mm) was between epochs 1 and 5 of camera model 1 (Table

Table 5 .
instead of individual cameras and dense image network as in this experiment.Though, having individual camera stability analysis over dense image network with wide variety of image rotation angles and calibration data from two days, this experiment establishes a reference guide on DF IOPs stability and sets the limits of the DF error budget.Also, it will help future investigation of the reason of individual camera model instability.RMSE of the 3D points differences (mm)