Accurate calibration scheme for a multi-camera mobile mapping system

: Mobile mapping systems (MMS) are increasingly used for many photogrammetric and computer vision applications, especially encouraged by the fast and accurate geospatial data generation. The accuracy of point position in an MMS is mainly dependent on the quality of calibration, accuracy of sensor synchronization, accuracy of georeferencing and stability of geometric conﬁguration of space intersections. In this study, we focus on multi-camera calibration (interior and relative orientation parameter estimation) and MMS calibration (mounting parameter estimation). The objective of this study was to develop a practical scheme for rigorous and accurate system calibration of a photogrammetric mapping station equipped with a multi-projective camera (MPC) and a global navigation satellite system (GNSS) and inertial measurement unit (IMU) for direct georeferencing. The proposed technique is comprised of two steps. Firstly, interior orientation parameters of each individual camera in an MPC and the relative orientation parameters of each cameras of the MPC with respect to the ﬁrst camera are estimated. In the second step the o (cid:11) set and misalignment between MPC and GNSS / IMU are estimated. The global accuracy of the proposed method was assessed using independent check points. A correspondence map for a panorama is introduced that provides metric information. Our results highlight that the proposed calibration scheme reaches centimeter-level global accuracy for 3D point positioning. This level of global accuracy demonstrates the feasibility of the proposed technique and has the potential to ﬁt accurate mapping purposes.


Introduction
Mobile mapping system (MMS) is a photogrammetric mapping agent that is usually defined as a set of navigation (global navigation satellite system (GNSS) and inertial measurement unit (IMU)) and remote sensors-such as cameras, lidar, and odometer sensors-integrated in a common moving platform [1].The importance of an MMS has been widely highlighted based on its cost-effectiveness, high data-capturing rate, and acceptable level of accuracy [2].System calibration is an indispensable part of the process of employing an MMS.Regarding an MMS that is based on multi-cameras and navigational systems, two calibration process could be mentioned: multi-camera calibration and MMS calibration.Multi-camera calibration aims to estimate interior and relative orientations of a multi-camera.MMS calibrations refers to estimating relative position and orientation between a multi-camera and GNSS/IMU sensors.Accurate system calibration ensures high-quality outputs for at-least a minimum period that a mapping system stays relatively still; a periodic system calibration scheme is able to guarantee the correctness of time-dependent parameters.
The first aspect of employing an optical-based MMS relates to the problem of multi-camera calibration [2][3][4][5][6][7][8][9][10][11].Nowadays, many MMSs are equipped with multi-projective cameras (MPC) because of their sturdy design, large field of view (FOV), and promising sensor models.Sensor modeling of an MPC usually consists of interior orientation parameters (IOP) of individual cameras, relative orientation parameters (ROP) between cameras with respect to a reference camera, and a scale factor that connects ROPs to a global framework.An important aspect of integrating a camera system in an MMS relates to the problem of employing a rigorous sensor model that fits to the physics of the camera.The sensor model maps a 3D object point into its corresponding image point [12].The next step is to employ the sensor model as the core of a statistical optimization model to find optimum values and express our uncertainties about the unknowns such as camera parameters or 3D position of object points.
Many MMSs use multiple cameras but do not necessarily generate panoramas.A panorama is a continuous presentation of an environment which is demonstrated by one photo or a series of photos that are merged together by a 'stitching' or a 'geometric non-stitching' approach.This form of photography provides viewers with an unlimited viewing possibility in all directions [12].Because of the all-angle viewing property, panoramic photography has found a wide range of applications; few of those applications are computer vision, robotics, surveillance, virtual reality, indoor/outdoor photography, and historical heritage documentation.
Recently developed photogrammetric models for MPCs has brought a new perspective in panoramic images applications.If these cameras, which are initially designed for all-direction photography, treated with photogrammetric models, numerous potential applications will emerge for surveying and mapping purposes.
A categorization of panoramic cameras can help us to find shared properties and mathematical model for members of each class.Amiri Parian and Gruen [12] categorized panoramic imaging into four groups: stitched, mirror-based rotating-head, near 180, and scanning 360 panoramas.Stitched panoramas are mainly used for non-metric applications where accurate directions are not important; therefore, this class is outside the focus of a metric categorization.We may re-categorize most of commercially available panoramic cameras that are suitable for surveying tasks according to their internal capturing technology into four groups [13]: (1) rotating-head, (2) multi-fisheye and (3) multi-projective camera, and (4) catadioptric systems.The benefits of the latter categorization are twofold.Most of commercially available cameras that are useful in metric photogrammetry tasks fit into this categorization; moreover, it simplifies the sensor modeling since it has founded based on a similarity measure that places cameras to a category according to their common sensor model.
The category (1) rotating head cameras contains cameras based on linear array CCDs that is mounted on a vertically rotating head.Examples of this class include EYESCAN M3 that is jointly developed by DLR and KST [14] and SpheroCam which is manufactured by SpheronVR AG.Equivalently, a projective camera such as a Canon EOS 6D could be mounted on a motorized rotating head such as Phototechnik AG, Roundshot metric, piXplorer 500, or GigaPan EPIC pro.The camera and step motors are synchronously controlled by a central processing Unit (CPU) that plans the movements and takes image shots.If the system has been accurately calibrated (optical center of the projective camera's lens should be precisely located on the center of rotation), high-precision metric information and high-resolution panoramic images can be obtained with least effort [15].Depending on the target applications, one or two step motors are working simultaneously to rotate the central camera around a point/axe.One main configurable parameter in building this class is the option to choose between one or two rotating axes; a post-processing module is responsible to transfer and merge images that have been captured with different orientations.Subsequently, a color blending approach is required to ensure that color tunes between stitched images are consistent.Usually a full resolution shot takes few seconds to even few minutes to compile which makes this class of panoramic cameras relatively unsuitable for mobile-mapping.Some important use-case of this class is cultural heritage recording and classification [16,17], high-precision surveying [12,[18][19][20][21], industrial-level visualization, indoor visualization, 3D modeling, and artistic photography.
The category (2) multi-fisheye cameras comprises a structure of a dual spherical camera module mounted on a rigid frame to construct a light-weight mobile 360 imaging system; Few consumer-grade samples of this class are Samsung Gear 360 [22], SVPro VR 360, MoFast 720, MGCOOL cam 360, SP6000 720, Ricoh Theta S [7].Some recent works suggests multi-camera system based on a dual-fisheye design [23,24].Fisheye camera models for photogrammetric applications were extensively studied and tested [25][26][27][28][29].The simultaneous geometric calibration of multi-camera system based on fisheye images, aiming a 360 • field of view also started to be discussed recently.For instance, a dual fisheye calibration model is proposed by [30] and [7].For this class of cameras, a customized statistical optimization process that involves using weighted observations and initial distributions of unknowns has proved sufficiently accurate for low to medium-accuracy surveying applications.However, some limitation of these systems can be mentioned, such as non-perspective inner geometry, huge scale and illumination variations between scenes, large radial distortion and nonuniformity of spatial resolution.Therefore, the overall image quality of panoramic images that are produced by the cameras of this class are usually worse than the first or third categories, however, some desired aspects-such as data-capturing rate or simplicity-makes them ideal candidates for applications such as low accuracy mobile mapping or 3D virtualization [31][32][33][34][35].
The category (3) multi-projective cameras contains panoramic cameras with arbitrary number of projective cameras mounted fix on a rigid frame (an industrial examples is LadyBug with six cameras [v.3 and v.5; MPC 6 ] by FLIR; a consumer-grade example is Panono 360 with 36 cameras [MPC 36 ]).Cameras of this class (multi-projective cameras or in short MPC) are usually customized for certain imaging, navigation, or surveying tasks.A synchronous shutter mechanism is applied to take simultaneous shots (<1 msec delay).A geometric model for MPC integrated into a statistical adjustment model is proposed by many researchers, e.g., ([9-11]).This model ensures desirable geometric accuracies for many tasks such as 3D mapping and surveying, 3D visualization, and texturing.Panoramas that are initially generated for this class are based on stitching techniques that mainly have visualization and artistic values; it is in contrast to the geometric values of panoramic images that are taken from cameras of the first or the second categories; this class of cameras has rigorously found many applications such as surveying, robotics, visualization, cinema, and artistic photography.Similar to the second class, it has the potential to be employed in applications that needs fast data-capturing rates such as mobile mapping, or navigation.
Early attempts to employ relative orientation constraints among multiple cameras was applied to a stereo camera, e.g., ([3,36,37]).He et al. [36] developed an MMS with a stereo camera and a GPS receiver to measure global coordinates of any point through photogrammetric intersection.King [3] modified the conventional bundle block adjustment (BBA) to accept relative orientation stability constraints.Zhuang [37] employed a fixed-length moving target to calibrate a stereo camera.Later, more complex systems of cameras went under research investigations.For example, Svoboda et al. [38] calibrated an MPC 4 that were installed on a visual room.Lerma et al. [39] calibrated a MPC 3 that consisted of a stereo camera and a thermal camera by employing distance constraints.Habib et al. [40] analyzed variations in IOP/Exterior Orientation Parameter (EOP) of multi-cameras.Detchev et al. [41] presented a system calibration scheme by employing system stability analysis.Some researchers employed a calibration field for multi-camera calibration.For example, Tommaselli et al. [5] employed AURUCO coded targets [42] to design a terrestrial calibration field.They used their proposed photogrammetric field to calibrate fisheye, catadioptric, and multi-cameras.Khoramshahi and Hokovaara [10] employed a customized coded target (CT) to create a calibration room.They employed it to calibrate a complex MPC 36 (Panono), and an MPC 3 (LadyBug v.3).In this work, we follow the calibration model that was proposed by Khoramshahi and Honkavaara [10].
Category (4) contains catadioptric systems that employ complex optical systems.A set of spherical and aspherical lenses, shaped mirrors such as parabolic, hyperbolic, or elliptical mirror [43], and refractors are employed in catadioptric systems to cover a large FOV.An example of this class was proposed by [23] as a prototype dual catadioptric camera.A calibration model for a catadioptric camera consists of a wide-lens camera and a conic mirror was proposed by [44].
The second aspect of employing an MMS relates to the problem of direct geo-referencing.Using multi cameras of third category for mapping applications require accurate determination of orientation and position of each panorama.This orientation can be done by indirect or direct georeferencing techniques, or combining both.Direct geo-referencing is a very useful approach for real-time application such as mobile mapping, or aerial surveying.Direct geo-referencing is the process to find position and orientation of captured images (EOPs) in a global reference frame without employing any ground control point (GCP), which requires the integration of additional sensors such as GNSS receivers and IMU sensors into the camera's mounted frame.This integration either provides initial values for positions and orientations of the camera shots as weighted observations in the BBA, or helps the system to instantly estimate position vectors and Euler angles of shots [45].It is worth noting that a successful direct geo-referencing requires two main aspects.Firstly, GNSS and IMU sensors should be accurately synchronized to an MPC.Secondly, a robust estimation of lever-arm vector and boresight angles should be known.These later parameters are known as mounting parameters [46].Rough estimation of the MMS calibration parameters could be performed by comparing observations of GNSS and IMU with the output of BBA.To enable a direct geo-referencing, a customized sensor model is required to be employed by considering additional parameters of lever-arm and boresight misalignments of an MPC with respect to GNSS and IMU sensors.
A variety of calibration schemes has been overly discussed in the literature; however, as far as the authors concerned, a rigorous, practical, and easy solution for direct georeferencing of an MPC does not exist that completely fit to the configuration of a multi-camera MMS.Moreover, stitching operation is required to compile a panorama from images of the third category cameras.This panoramic presentation usually has just artistic or visual values.A non-stitching panoramic creation scheme is also feasible for this class that adds geometric value to generated panoramas.This paper contributes to a rigorous, easy, and practical scheme for calibrating a multi-camera MMS equipped with an MPC, GNSS and IMU sensors in a terrestrial vehicle.We also present a novel non-stitching algorithm by introducing a panoramic correspondence map.We propose a modified BBA that comprises MPC calibration and MMS calibration.We assess quantitatively the quality of the proposed approach by checking the accuracy of the point intersection by independent checkpoints.Our results demonstrate that sub-decimeter level accuracy is achievable by this technique.

System Calibration
Mathematical theory behind the implementation of this work is discussed in this part.We first start with sensor model of an MPC, then we continue toward discussing a complete sensor model for the underlying MMS.Finally, a modified BBA is discussed to employ the sensor model and statistically express uncertainties about parameters.

Individual Camera Model
The following interior orientation model removes main distortions from image coordinates and outputs undistorted pinhole image coordinates as where (x) t 1 and (y) t 1 are distorted coordinates of a point in pixel unit, PP x and PP y are coordinates of the principal point, K i are radial distortion coefficients, P i are tangential distortion parameters, and δ and λ are scale and shear factor respectively.It is trivial to show that the underlying model (Equation ( 1)) is exchangeable to Brown's model [47,48].
An image of a projective camera is linked to the camera's sensor model, and relates to six orientation and position parameters that uniquely determines it inside a 3D Euclidean space.Therefore, a linear relationship between an object point and its undistorted pinhole coordinates is established by the collinearity model where M i is the i th row of the rotation matrix M 3×3 , X is a vector of (3 × 1) of the object point's coordinates, and X 0 is a (3 × 1) vector of location of image in a 3D cartesian space (bold letters are used for vectors and matrices and normal letters are used for scalers throughout this paper).

Multi-Projective Camera Model
An MPC is presented in this work by the following calibration parameters: where IOP i is interior orientation parameters of an individual camera of an MPC, ROP i is relative orientation of an individual camera with respect to the first camera, (ζ, η, ψ) are Euler angles, ∆ x , ∆ y , ∆ z are displacements, and Λ MPC is the unknown scale factor.The last parameter (Λ MPC ) plays a role if scale of an MPC sensor would be different from scale of the model.It could be ignored if both scales are equal.In overall, there will be a number of nCalib.= n × 10(IOP) + (n − 1) × 6(ROP) + 1(scale) = 16 × n − 5 calibration parameters for an MPC with (n) projective camera; a multi projective image (MPI) is geometrically defines by its corresponding MPC sensor and six EOPs.

Space Resection of Multi-Projective Images
The goal of the space resection is to find approximate locations of a given MPI with respect to a 3D frame under the conditions that at-least (3) image-object correspondences exist and localized in a 3D Cartesian framework.Space resection helps to orient a set of MPIs with respect to a given point set (e.g., a calibration room).The orientations are finally adjusted through a BBA.
An alternative way to resect an MPI is to employ the coplanarity equation.Coplanarity condition connects two normalized pinhole coordinates and image centers through a matrix called Essential, which is a (3 × 3) matrix that is composed of three rotations and three translations as where [(X 0 ) t2 ] x is the matrix presentation of a cross product which is a 3 × 3 matrix of the form Essential matrix is estimated by the following equations where x t 1 is normal image point at time (t 1 ), and x t 2 is the calculated normal location of the same 3D point in a synthesized image.Space resection of an MPI is performed for each of its projective images.At least 5-8 tie points are required based on the selected method (5-point [49], 6-point, 8-point, or normalized 8-point [50]).Then a synthesized pair is formed.The synthesized image is localized such that the resulted stereo pair becomes stable.Next, the retrieved position is translated into approximate position of the parent MPI (X0 i , wpk i ); finally, all positions are averaged to estimate the position and orientation of the MPI

Relative Orientation between Two Multi-Projective Images
Finding ROPs between two MPIs is essential to independently construct a scene.ROPs between two MPIs are expressed by six parameters out of which five are independent.Usually the largest element of a stereo pair's baseline is scaled to one.Three rotations and two displacements are finally solved by relative orientation module.To estimate ROPs between two MPIs, a similar averaging approach to the space resection has been chosen.Therefore, all image pairs between MPIs are oriented.Then the resultant ROPs are transformed to the parent MPI and averaged.
In the process of MMS calibration, initial values of lever-arm and boresight misalignments are estimated by independently compiling a local network of few shots, then connecting the local network to the global frame of GNSS/IMU sensors.

Bundle Block Adjustment
BBA employs non-linear least-square adjustment method to minimize a cost function based on residuals of observational equations.It finally expresses the uncertainties about the unknowns of an MPC calibration as a covariance matrix.The MPC cost function is a combination of all observational equations.It is expressed as the vector function (7) where, (f, PP, K, P, δ, λ ) i=1:n are IOPs for projective cameras (1:n), R (ζ,η,ψ) i=2:n , ∆ i=2:n are structural parameters regarding the MPC, R (ω,φ,κ) t=t 1 :tm , (X 0 ) t=t 1 :t m are orientations and locations of image shots, X (1:nO) , x (1:nI) are position of object points and the corresponding image points, la and bs are lever-arm vector and boresight misalignments respectively.To model the entire system, equations for image observations, GNSS/IMU-observations, and GCPs are introduced.
Two observational equations are considered in the MPC's BBA for each image tie points, as In Equations ( 8) and ( 9), x i j , y i j are pinhole coordinates of object point (i) inside image (j) (pinhole camera, f = 1), (M i ) t 1 j is the (i th ) row of the rotation matrix M t 1 j (j th projective camera [ j : 1 − n] at time t 1 ), (X 0 ) t 1 j is location of j th projective camera at time t 1 , and X i is 3D location of the corresponding object point (i).
In order to combine observations of GNSS and IMU sensors to the BBA, a link between GNSS, IMU and the MPC is established (Figure 1).A number of ( 6) additional parameters are introduced to represent a shift (lever-arm), and orientations (boresight misalignments).The connection is formulated by six non-linear equations (Equations ( 10) and ( 11)).The observational equation for lever-arm and boresight misalignments are where GCPs are contributed to the BBA by three observational equations where is unknown 3D position of a GCP, and is observed position.

Panoramic Compilation
Creating a panorama for an MPC without image stitching is possible under the condition that the MPC is calibrated (known IOPs and ROPs).Therefore, calibration data is sufficient to compile a non-stitching panorama.The algorithm for building a panorama contains looping over a hypothetical sphere or cylinder (optional projection system) that covers an MPI, then for each pixel in the final  GCPs are contributed to the BBA by three observational equations where X g is unknown 3D position of a GCP, and X go is observed position.

Panoramic Compilation
Creating a panorama for an MPC without image stitching is possible under the condition that the MPC is calibrated (known IOPs and ROPs).Therefore, calibration data is sufficient to compile a non-stitching panorama.The algorithm for building a panorama contains looping over a hypothetical sphere or cylinder (optional projection system) that covers an MPI, then for each pixel in the final panoramic compilation the algorithm finds a corresponding camera index (0-n) and an image location (pix).For a better rendering result, an appropriate interpolation technique (bilinear or bicubic) could be considered since the projection from the sphere (or cylinder) to projective images results in subpixel positions.If the correspondence data is saved as a meta-data file (it may be called panoramic correspondence map), then the final panorama will have geometric values.It is because of the fact that collinearity condition between an object point, its corresponding image point, and perspective center of a projective camera of an MPC is established.A second usage of the correspondence map is efficiently compiling a panorama, since indexing is expected to be much quicker than direct calculation.
There are regions in the correspondence map that more than one choice of correct image is possible for every point in the final panoramic compilation.It is because of the existence of overlaps between adjacent cameras.A criterion needs to be considered to build a unique map, for example minimum incident angle could be used.
Depending on the distance of a 3D object points to a panorama, discontinuity on edges is expected to occur; it is because the projection center of different individual cameras does not match to the center of the hypothetical sphere.Discontinuity is expected to be more severe for closer objects.Finally, an important post-processing is the use of color blending on edges.It improves the quality of a compiled panorama.Without color blending, sharp steps are expected to be observed on the edges that pixel labels start to switch from one camera to another.Color blending softens those steps and improves the color consistency of the final panorama.

Mobile Mapping System
The MMS in this works comprises an MPC, a GNSS, and an IMU that are firmly installed on a car's roof.This section describes the experiments to assess the proposed method described in Section 2. The proposed method was implemented in C++ and MATLAB languages and evaluated using a real dataset from an MMS composed by a multi-camera LadyBug 5, and a navigation system installed on a rigid aluminum truss structure on top of a Skoda Superb (Figure 2).
Depending on the distance of a 3D object points to a panorama, discontinuity on edges is expected to occur; it is because the projection center of different individual cameras does not match to the center of the hypothetical sphere.Discontinuity is expected to be more severe for closer objects.Finally, an important post-processing is the use of color blending on edges.It improves the quality of a compiled panorama.Without color blending, sharp steps are expected to be observed on the edges that pixel labels start to switch from one camera to another.Color blending softens those steps and improves the color consistency of the final panorama.

Mobile Mapping System
The MMS in this works comprises an MPC, a GNSS, and an IMU that are firmly installed on a car's roof.This section describes the experiments to assess the proposed method described in Section 2. The proposed method was implemented in C++ and MATLAB languages and evaluated using a real dataset from an MMS composed by a multi-camera LadyBug 5, and a navigation system installed on a rigid aluminum truss structure on top of a Skoda Superb (Figure 2).The GNSS receiver used was a NovAtel PwrPak7, which together with IMU compiles a Synchronized Position Attitude Navigation (SPAN) system.PwrPak7 logged GPS and Glonass observations via NovAtel VEXXIS GNSS-850 satellite antenna.Sensor trajectory was computed by combining GNSS and IMU data with Waypoint Inertial Explorer software.Virtual GNSS base station reference data was acquired from commercial Trimnet service.The maximum distance from trajectory (Figure 3) to the base station was ~350 meters.The 3D RMSE for the GPS positions was 1.6 centimeters.A NovAtel IMU-ISA-100C, manufactured by Northrop-Grumman Litef GMBH, was paired with a NovAtel SPAN GNSS receiver.The IMU-ISA-100C is a near navigation grade sensor containing GNSS-850 satellite antenna (in front of IMU) mounted on a truss structure on top of a Skoda passenger car.NovAtel PwrPak7 GNSS receiver together with operating and logging laptop were inside the car.

MPC Calibration
FGI's calibration room is an empty space with dimension of 519 × 189 × 356 cm 3 that is equipped with 215 CTs (Figure 4).3D positions of all targets are a-priori known based on previous observations.Images taken with a Canon EOS 6D camera (sensor 20 Megapixel, image size 5472 × 3648 pixels, lens Canon EF 24 mm, focal length 20.6 mm) were used to accurately calculate 3D coordinates of the targets [10] in the calibration room with bundle adjustment, before MPC calibration.
In the first step, FGI's camera-calibration room was employed to calibrate the MPC.To do so, 3D positions of the CTs had been accurately estimated by employing a Canon EOS-6D camera achieving a relative 3D positioning precision of std < 1 mm.The automatic target extraction included: A NovAtel IMU-ISA-100C, manufactured by Northrop-Grumman Litef GMBH, was paired with a NovAtel SPAN GNSS receiver.The IMU-ISA-100C is a near navigation grade sensor containing GNSS-850 satellite antenna (in front of IMU) mounted on a truss structure on top of a Skoda passenger car.NovAtel PwrPak7 GNSS receiver together with operating and logging laptop were inside the car.

MPC Calibration
FGI's calibration room is an empty space with dimension of 519 × 189 × 356 cm 3 that is equipped with 215 CTs (Figure 4).3D positions of all targets are a-priori known based on previous observations.Images taken with a Canon EOS 6D camera (sensor 20 Megapixel, image size 5472 × 3648 pixels, lens Canon EF 24 mm, focal length 20.6 mm) were used to accurately calculate 3D coordinates of the targets [10] in the calibration room with bundle adjustment, before MPC calibration.The estimated focal length of the Canon was 3156.10 pixels with a standard deviation of 0.04 pixels.The principal point was estimated as (2752.37,1509.36)pixels with standard deviation of 0.05 pixels.Regarding each of the IOPs, a significance analysis is performed assessing the amount of distortion that is corrected by freezing all the IOPs except for an underlying parameter (or set of parameters) on a grid of 10 × 10 pixels.The largest distortion on the grid was chosen as significance value for a given IOP.The significance values for radial distortions (set), tangential distortion (P1, P2, and P3), scale, and shear were 69.2, 2.4, 1.0, 0.03, and 0.06 pixel respectively.The estimated standard deviations of the CT points were 0.15, 0.03, and 0.02 mm in principal directions of 3D error ellipsoids respectively, which indicates that our target accuracy was achieved.
In the second step, the same CTs were automatically extracted in individual images of the MPC (LadyBug 5).Next, each MPI was resected (Chapter 2.3) using the calibration-room's CTs, and finally, the BBA was performed to optimize the sensor parameters of the MPC.In this processing, the CTs coordinates based on the measurement by the Canon were kept fixed.
Two datasets of 26 and 53 MPIs were captured from the FGI's calibration room to calibrate the MPC.In the first dataset, the camera is attached to a horizontal slide, therefore each few images are coplanar and at the same height level.In the second dataset, the camera is put on a tripod with two different height levels.In total, 474 individual projective images have been captured; then, CTs were automatically extracted in all individual images.

MMS Calibration
The MMS calibration was carried out using outdoor dataset from Inkoo, which is a municipality in southern Finland.A set of 8424 MPIs were captured by the LadyBug 5 camera for Inkoo harbour area and two close-by regions.Most roads have been captured by the MMS in both directions.In the first step, FGI's camera-calibration room was employed to calibrate the MPC.To do so, 3D positions of the CTs had been accurately estimated by employing a Canon EOS-6D camera achieving a relative 3D positioning precision of std < 1 mm.The automatic target extraction included: 1.
Employing adaptive thresholding using first-order statistics to roughly estimate the position of blobs (accuracy of better than 2pixel); 2.
Finding CTs from extracted clusters by employing structural signature of the CTs, and 5.
Reading ID for each target.
The automatic CT detection was implemented in MATLAB.The standard deviation for image-point observations was set as 0.1 pixel.
The estimated focal length of the Canon was 3156.10 pixels with a standard deviation of 0.04 pixels.The principal point was estimated as (2752.37,1509.36)pixels with standard deviation of 0.05 pixels.Regarding each of the IOPs, a significance analysis is performed assessing the amount of distortion that is corrected by freezing all the IOPs except for an underlying parameter (or set of parameters) on a grid of 10 × 10 pixels.The largest distortion on the grid was chosen as significance value for a given IOP.The significance values for radial distortions (set), tangential distortion (P1, P2, and P3), scale, and shear were 69.2, 2.4, 1.0, 0.03, and 0.06 pixel respectively.The estimated standard deviations of the CT points were 0.15, 0.03, and 0.02 mm in principal directions of 3D error ellipsoids respectively, which indicates that our target accuracy was achieved.
In the second step, the same CTs were automatically extracted in individual images of the MPC (LadyBug 5).Next, each MPI was resected (Section 2.3) using the calibration-room's CTs, and finally, the BBA was performed to optimize the sensor parameters of the MPC.In this processing, the CTs coordinates based on the measurement by the Canon were kept fixed.
Two datasets of 26 and 53 MPIs were captured from the FGI's calibration room to calibrate the MPC.In the first dataset, the camera is attached to a horizontal slide, therefore each few images are coplanar and at the same height level.In the second dataset, the camera is put on a tripod with two different height levels.In total, 474 individual projective images have been captured; then, CTs were automatically extracted in all individual images.

MMS Calibration
The MMS calibration was carried out using outdoor dataset from Inkoo, which is a municipality in southern Finland.A set of 8424 MPIs were captured by the LadyBug 5 camera for Inkoo harbour area and two close-by regions.Most roads have been captured by the MMS in both directions.
The GCPs were measured using the Topcon Hiper HR RTK dual band receiver with an accuracy within 5 mm + 0.5 ppm horizontally and 10 mm + 0.8 ppm vertically [1.source:TOPCON].Nearest base station was located in Metsähovi 28 km away from Inkoo, therefore accuracy is vertically 19.0 mm and horizontally 32.4 mm, respectively.
The MMS calibration involved estimating relative parameters of the MPC with respect to the GNSS and IMU (lever-arm and boresight miss-alignments).Observational Equations ( 10) and (11) were added to BBA.Analytical partial derivatives of latter equations were considered for calculating the Jacobian matrix.The calibrated MPC sensor was used to construct an outdoor case.The sensor was kept fixed during the scene reconstruction except for its scale parameter.A number of six MPIs were connected together by the relative orientation module (Section 2.4) to construct a local network (initial network).Model coordinates of all 3D points were accurately determined with corresponding sub-pixel re-projection error and added as observational equations.Few GCPs were observed in individual images and added to the model to establish a link to the global coordinate system (Equation ( 6)).Then, GNSS and IMU observations were added to the initial network.Lever-arm and boresight was initialized from the initial network and finally optimized through the BBA.The MPC sensor kept fixed during MMS calibration to avoid a singularity in the Jacobian matrix.In the optimization process, locations, and orientations of MPIs, 3D positions of GCPs, and six parameters of lever-arm and boresight were set free.

Performance Assessment
In order to assess the accuracy of the MMS calibration and georeferencing, seven check sites approximately covering the area were selected (Figure 3).In each site, the GNSS and IMU readings of the MMS were translated into MPI orientations that were converted into individual image's orientation and rotation by using MPC sensor calibration information from Step 1 (Section 3.2).Then image positions of check points were observed in individual images; consequently, two sets of coordinates were available for each cross-check site: 1-intersecting coordinates from image measurements and collinearity equations; and 2-observed positions of check-points from GPS.The difference was recorded along with the maximum intersecting angle of each 3D check point as a quality-control indicator.It was expected that weaker interesting geometry would result higher object-space residuals.

MPC Calibration
Indoor calibration of the MPC resulted in sub-pixel individual image accuracies (mean image residual 0.4 pixel, std.0.45 pixel).Moreover, internal structure of the final optimized MPC resembled very well to the physical reality of the underlying MPC.Different height levels of MPIs resembled well to the reality of the capturing configuration and played a role as a visual cross-checking.Undistorted individual image, as depicted in Figure 5 had very large radial distortions.The MPC's significance matrix also showed a considerable value of 2183 pixel for radial distortion for camera number (1) that was caused by the lens structure.Significance values for tangential distortions were relatively small (<10 pixel).The calibrated IO parameters of individual cameras with their corresponding standard deviation are given in Table 1.Most of the standard deviation values of camera 6 were higher than average.One reason could be the weaker distribution of CTs on the ceiling in comparison to the walls and the floors distributions.Scale and shear factors appeared meaningful in the BBA since their optimized values were relatively bigger than the standard deviation values.Table 2 demonstrates the calibrated structure of the LadyBug camera.The order of standard deviation values demonstrated an acceptable level of precision in the calibration process.In this table, the values related to the first camera demonstrates a priori distribution of the camera 1's structural parameters.In this table, standard deviation values that are bigger than average are highlighted with a gray color.Figure 6 demonstrates absolute values of normalized cross-correlations between all free parameters (IOPs and ROPs (structure)) of the camera 2 (indoor calibration).In this figure off-diagonal elements are of high importance and need to be carefully treated in the adjustment to avoid singularity.Main off-diagonal elements are marked by red rectangles.In this figure, radial distortions are the first group of highly correlated parameters.This high correlation implies that the collective effect of the group is important, therefore investigating an individual parameter in this group is irrelevant.The second highly correlated group is ROPs {ζ and η}.This group has a strong correlation to PP x .The last high correlation relates to ROP {ψ} that is highly correlated to PP y .Figure 7 demonstrates absolute values of normalized cross-correlations between all free parameters of the MPC.This figure demonstrates a slightly repeated pattern among cameras.Two diagonal parts of this figure are zero that are related to the structure of the first camera and lever-arm vector and boresight misalignment that kept fixed during the indoor calibration.

13
this figure are zero that are related to the structure of the first camera and lever-arm vector and boresight misalignment that kept fixed during the indoor calibration.

Mobile Mapping System Calibration
Table 3 shows the calibrated values for the lever-arm vector and boresight misalignments and their standard deviation values.As expected, standard deviation values were larger than structural parameters (ROPs of the MPC).On average, 0.5-degree standard deviation were obtained for boresight misalignments, and 2 cm for each component of lever-arm vector.
-2.4189 -0.2824 0.7361 0.02 0.02 0.02 In the adjustment processing, on average 9 cm residuals were obtained for the GCPs (Table 4) on the outdoor calibration site (Figure 8).All GCPs had at-least 6 intersecting rays.Table 5 shows the differences between the image positions and orientations estimated by the BBA and observed by the GNSS/IMU.On average, approximately 0.2° difference was observed between GNSS/IMU angles and BBA outputs.For the directly measured positions, an average difference of 11.3, 2.8, 0.6 cm in the X, Y, and Z coordinates, respectively, were observed.

Mobile Mapping System Calibration
Table 3 shows the calibrated values for the lever-arm vector and boresight misalignments and their standard deviation values.As expected, standard deviation values were larger than structural parameters (ROPs of the MPC).On average, 0.5-degree standard deviation were obtained for boresight misalignments, and 2 cm for each component of lever-arm vector.In the adjustment processing, on average 9 cm residuals were obtained for the GCPs (Table 4) on the outdoor calibration site (Figure 8).All GCPs had at-least 6 intersecting rays.Table 5 shows the differences between the image positions and orientations estimated by the BBA and observed by the GNSS/IMU.On average, approximately 0.2 • difference was observed between GNSS/IMU angles and BBA outputs.For the directly measured positions, an average difference of 11.3, 2.8, 0.6 cm in the X, Y, and Z coordinates, respectively, were observed.2) respectively and Tables 6 and 7 shows the errors for check-points, which are differences between their intersected and observed values.On average, error of 5.6 cm and 1.4 cm was observed between intersected and observed positions of check points for cross-check dataset (1) and ( 2), respectively.The RMSEs were   2) respectively and Tables 6  and 7 shows the errors for check-points, which are differences between their intersected and observed values.On average, error of 5.6 cm and 1.4 cm was observed between intersected and observed positions of check points for cross-check dataset (1) and (2), respectively.The RMSEs were 0.69 cm and 0.18 cm, respectively.Relatively comparable GPS errors for close-by observations resulted to a similarity between errors in 3D check points of a dataset.16 0.69 cm and 0.18 cm, respectively.Relatively comparable GPS errors for close-by observations resulted to a similarity between errors in 3D check points of a dataset.16 0.69 cm and 0.18 cm, respectively.Relatively comparable GPS errors for close-by observations resulted to a similarity between errors in 3D check points of a dataset.Figure 11 plots sorted incident angles (degrees) with their corresponding errors.This plot highlights the effect of intersection-geometry's strength on the quality of positioning.In this plot, intersections of rays with maximum incident angle less than 20 degrees resulted in higher fluctuation in 3D-position inaccuracy.The inaccuracy in 3D positions was significantly improved when the incident angles took values over 20 degree.A logical conclusion from this figure is that a stronger intersection geometry is obtained as incident angles come closer to .Figure 12a shows errors for 3D check point.Figure 12b shows the corresponding intersecting angles.In Figure 12a, most errors were less than 7 cm.By comparing these two figures, a negative correlation between incident angles and accuracy was obvious which was also concluded from Figure 11.On average, 4.2 cm error was observed for all 3D check points with RMSE 3.6 cm.

18
Figure 12 (A) shows errors for 3D check point.Figure 12 (b) shows the corresponding intersecting angles.In Figure 12 (a), most errors were less than 7 cm.By comparing these two figures, a negative correlation between incident angles and accuracy was obvious which was also concluded from Figure 11.On average, 4.2 cm error was observed for all 3D check points with RMSE 3.6 cm.

Panoramic Compilation
Panoramic images were compiled using the IOPs and ROPs. Figure 13 demonstrates footprints of individual cameras on the final panoramic compilation as a contribution map of individual cameras.This map demonstrates that each individual image covers a suitable amount of area in the final panoramic compilation.It consequently confirms that the underlying MPC is well-designed.
The 'minimum incident angle' criterion is employed to create a unique correspondence map (Figure 14) from the contribution map (Figure 13).The correspondence map connects every pixel in the final panoramic compilation to a pixel in a camera.Therefore, a vector of 3 numbers (image id, x and y in pixel) were saved for every pixel of the final panorama as the correspondence map.Approximately 7 minutes of computational time was required to build a correspondence map of 7200 × 3600 pixels; then it was efficiently quick and straight-forward to compile any new panorama from individual images (~2.1 sec.).Non-stitching panoramic compilation was performed by employing the correspondence map and resulted into geometrically-accurate panoramas (Figure 15).Closer objects to the MPC had, as expected, larger discontinuity over the edges than far objects.An edge in this context means the boundary that pixel labels change.The discontinuity on edges has been considered by [51] as systematic errors; however, by using the correspondence map this systematic error is totally irrelevant and is completely eliminated and avoided.

Panoramic Compilation
Panoramic images were compiled using the IOPs and ROPs. Figure 13 demonstrates footprints of individual cameras on the final panoramic compilation as a contribution map of individual cameras.This map demonstrates that each individual image covers a suitable amount of area in the final panoramic compilation.It consequently confirms that the underlying MPC is well-designed.The 'minimum incident angle' criterion is employed to create a unique correspondence map (Figure 14) from the contribution map (Figure 13).The correspondence map connects every pixel in the final panoramic compilation to a pixel in a camera.Therefore, a vector of 3 numbers (image id, x and y in pixel) were saved for every pixel of the final panorama as the correspondence map.Approximately 7 minutes of computational time was required to build a correspondence map of 7200 × 3600 pixels; then it was efficiently quick and straight-forward to compile any new panorama from individual images (~2.1 sec.).Non-stitching panoramic compilation was performed by employing the correspondence map and resulted into geometrically-accurate panoramas (Figure 15).Closer objects to the MPC had, as expected, larger discontinuity over the edges than far objects.An edge in this context means the boundary that pixel labels change.The discontinuity on edges has been considered by [51] as systematic errors; however, by using the correspondence map this systematic error is totally irrelevant and is completely eliminated and avoided.

Discussion
In this work, we followed a similar technique as [11] in employing free-network calibration inside an indoor calibration room; however, we pushed this work in three directions.Firstly, our modified BBA model contains offset and misalignment with respect to the GNSS/IMU.Secondly, a non-stitching panorama was introduced along with a correspondence map which adds geometric values to the non-stitching panorama scheme.The new scheme improves the development of [51] by taking the systematic error of edges into equations.Thirdly, we demonstrated error propagations through a statistical model as standard deviation values of the parameters under-investigation.

Discussion
In this work, we followed a similar technique as [11] in employing free-network calibration inside an indoor calibration room; however, we pushed this work in three directions.Firstly, our modified BBA model contains offset and misalignment with respect to the GNSS/IMU.Secondly, a non-stitching panorama was introduced along with a correspondence map which adds geometric values to the non-stitching panorama scheme.The new scheme improves the development of [51] by taking the systematic error of edges into equations.Thirdly, we demonstrated error propagations through a statistical model as standard deviation values of the parameters under-investigation.

Discussion
In this work, we followed a similar technique as [11] in employing free-network calibration inside an indoor calibration room; however, we pushed this work in three directions.Firstly, our modified BBA model contains offset and misalignment with respect to the GNSS/IMU.Secondly, a non-stitching panorama was introduced along with a correspondence map which adds geometric values to the non-stitching panorama scheme.The new scheme improves the development of [51] by taking the systematic error of edges into equations.Thirdly, we demonstrated error propagations through a statistical model as standard deviation values of the parameters under-investigation.
A complete calibration scheme for a multi-camera mobile mapping system (MMS) was presented to calibrate a multi projective camera (MPC) to GNSS and IMU.It was based on two steps: indoor MPC calibration, and outdoor MMS calibration, which calibrated an MPC with respect to the GNSS and IMU.Most of recent MPC calibration schemes are based on reducing the number of interior orientation parameters (IOP) in order to avoid singularity.In this work, the singular situation was addressed by employing a photogrammetric calibration room.Our modified bundle block adjustment model considers offset and misalignment of an MPC with respect to the GNSS/IMU, which enables 3D reconstruction using the MMS either in a direct georeferencing mode, or via integrated bundle block adjustment.Furthermore, a non-stitching panoramic compilation scheme was introduced along with a correspondence map which connects pixels of a compiled panorama to their correct position in the MPC.The new scheme takes the systematic error of discontinuous edges into equations.We finally demonstrated error propagations through a statistical model as standard deviation values of the parameters under-investigation.The statistical modeling was proved as a helpful tool to propagate the uncertainties from observations to unknowns based on a non-linear least-square model.
are 'observed' Euler angles of the INS at time (I t ), (R C ) I is the rotation matrix of an MPC with respect to the INS local coordinate system, and R C t is rotation matrix of the MPC at time (t).In Equation (11), X0 I t is observed position of INS at time (t), X C t is unknown position of the MPC at time (t), and (X C ) I is relative position of the camera with respect to the INS in the local coordinate system of the INS.Here, R indicates the function that maps a rotation matrix into its corresponding Euler angles.( ).In Equation (11), is observed position of INS at time ( ), is unknown position of the MPC at time ( ), and ( ) is relative position of the camera with respect to the INS in the local coordinate system of the INS.Here, ℝ indicates the function that maps a rotation matrix into its corresponding Euler angles.

Figure 1 .
Figure 1.Configuration of GNSS, IMU and the multi-camera.

Figure 1 .
Figure 1.Configuration of GNSS, IMU and the multi-camera.

Figure 2 .
Figure 2. Ladybug 5+ (in front), NovAtel IMU-ISA-100C (black-and-white box in rear) and NovAtel VEXXIS fiber optic gyros and temperature compensated Micro Electromechanical Systems (MEMS) accelerometers.IMU measurements are used by the SPAN receiver to compute a blended GNSS+INS position, velocity, and attitude solution at rates up to 200 Hz.Ladybug V.5 is an MPC with five side-looking cameras and one upward-looking camera.Each individual camera contains a sensor of size 2464 × 2048 pix.Its focal length is 18 mm with sensor size of (35.8 × 29.7 mm) with field of view of 89.3° ± 81′.Panoramas cover an area of 360° horizontally by 145° vertically.Resolution of output panoramas is approximately 8k (7200 × 3600 pixels, aspect ratio 2:1).

Figure 2 .
Figure 2. Ladybug 5+ (in front), NovAtel IMU-ISA-100C (black-and-white box in rear) and NovAtel VEXXIS fiber optic gyros and temperature compensated Micro Electromechanical Systems (MEMS) accelerometers.IMU measurements are used by the SPAN receiver to compute a blended GNSS+INS position, velocity, and attitude solution at rates up to 200 Hz.

Figure 3 .
Figure 3. Orthophoto of central Inkoo area.Blue rectangle is the outdoor calibration dataset.Red rectangles are check sites.Orange dots are locations of captured multi-images.

Figure 3 .
Figure 3. Orthophoto of central Inkoo area.Blue rectangle is the outdoor calibration dataset.Red rectangles are check sites.Orange dots are locations of captured multi-images.

Figure 5 .
Figure 5.A distorted image (a) and its corresponding undistorted image (b) from camera 1 of the Ladybug V5.

Figure 6 .
Figure 6.Absolute values of normalized cross-correlations for 'camera 2' after indoor calibration of the LadyBug camera.

Figure 5 .
Figure 5.A distorted image (a) and its corresponding undistorted image (b) from camera 1 of the Ladybug V5.

Figure 5 .
Figure 5.A distorted image (a) and its corresponding undistorted image (b) from camera 1 of the Ladybug V5.

Figure 6 .
Figure 6.Absolute values of normalized cross-correlations for 'camera 2' after indoor calibration of the LadyBug camera.

Figure 6 . 14 Figure 7 .
Figure 6.Absolute values of normalized cross-correlations for 'camera 2' after indoor calibration of the LadyBug camera.

Figure 7 .
Figure 7. Absolute values of normalized cross-correlations for 'all cameras' after indoor calibration of the LadyBug camera.

Figure 8 .
Figure 8. Configuration of the outdoor calibration dataset.(a) 3D project, (b) orthophoto overlaid with image positions and GCPs.

Figure 8 .
Figure 8. Configuration of the outdoor calibration dataset.(a) 3D project, (b) orthophoto overlaid with image positions and GCPs.

Figure 11
Figure11plots sorted incident angles (degrees) with their corresponding errors.This plot high-lights the effect of intersection-geometry's strength quality of positioning.In this plot, intersections of rays with maximum incident angle less than 20 degrees resulted in higher fluctuation in 3D-position inaccuracy.The inaccuracy in 3D positions was significantly improved when the incident angles took values over 20 degree.A logical conclusion from this figure is that a stronger intersection geometry is obtained as incident angles come closer to π 2 .

Figure 11 .
Figure 11.Sorted re-projection errors and their corresponding incident angles.

Figure 11 .
Figure 11.Sorted re-projection errors and their corresponding incident angles.

Figure 12 .
Figure 12.Error in check points in all areas.(a) Differences between measured GPS positions and direct-georeferencing positions, (b) the corresponding intersecting angles.

Figure 12 .
Figure 12.Error in check points in all areas.(a) Differences between measured GPS positions and direct-georeferencing positions, (b) the corresponding intersecting angles.

Figure 13 .
Figure 13.Contribution and overlaps of all cameras pixels (1-6) in the final panoramic compilation.

Figure 14 .
Figure 14.Correspondence map (min incident angle criteria) of all images (1-6) in the final panoramic compilation.

Figure 13 .
Figure 13.Contribution and overlaps of all cameras pixels (1-6) in the final panoramic compilation.

Figure 13 .
Figure 13.Contribution and overlaps of all cameras pixels (1-6) in the final panoramic compilation.

Figure 14 .
Figure 14.Correspondence map (min incident angle criteria) of all images (1-6) in the final panoramic compilation.

Figure 15 .
Figure 15.Compilation of a non-stitching panorama (no color blending).

Figure 14 .
Figure 14.Correspondence map (min incident angle criteria) of all images (1-6) in the final panoramic compilation.

Figure 13 .
Figure 13.Contribution and overlaps of all cameras pixels (1-6) in the final panoramic compilation.

Figure 14 .
Figure 14.Correspondence map (min incident angle criteria) of all images (1-6) in the final panoramic compilation.

Figure 15 .
Figure 15.Compilation of a non-stitching panorama (no color blending).

Figure 15 .
Figure 15.Compilation of a non-stitching panorama (no color blending).

Table 1 .
Calibrated interior orientation of individual cameras of the LadyBug camera.

Table 2 .
Calibrated structure of the LadyBug camera (orientations in degree and positions in (cm)).

Table 3 .
Calibrated values for lever-arm vector and boresight misalignments

Table 3 .
Calibrated values for lever-arm vector and boresight misalignments.

Table 4 .
Residuals for the ground control points of the outdoor calibration dataset

Table 5 .
Differences between adjusted values from BBA and direct-georeferencing for the outdoor calibration dataset

Table 4 .
Residuals for the ground control points of the outdoor calibration dataset.

Table 5 .
Differences between adjusted values from BBA and direct-georeferencing for the outdoor calibration dataset.

Table 6 .
Differences between observed and intersected position of check points in check1 dataset

Table 6 .
Differences between observed and intersected position of check points in check1 dataset

Table 6 .
Differences between observed and intersected position of check points in check1 dataset.

Table 7 .
Differences between observed and intersected position of check points in check2 dataset.

Table 7 .
Differences between observed and intersected position of check points in check2 dataset