Visual-Inertial Odometry

Visual-inertial odometry (VIO) estimates the motion of a camera that is rigidly attached to an inertial measurement unit (IMU), by fusing the images with the IMU’s gyroscope and accelerometer readings. It is visual odometry with a second, very different sensor: the camera sees the scene but, alone, gives a single-camera trajectory only up to an unknown scale; the IMU measures rotation rate and acceleration at hundreds of samples per second, in metric units and relative to gravity, but drifts quickly when integrated on its own. The robotics literature also calls the problem vision-aided inertial navigation (Mourikis and Roumeliotis, 2007; Hesch et al., 2014) or a visual-inertial navigation system (VINS), the name behind VINS-Mono. Forster et al. (2017) summarize why the pairing works: both sensors are cheap and complementary, and the IMU makes the metric scale of monocular vision and the direction of gravity observable.

This article covers the task: the sensor models, IMU preintegration, filter-based estimators such as the multi-state constraint Kalman filter (MSCKF) against optimization-based ones such as OKVIS and VINS-Mono, what an IMU can and cannot make observable, initialization, and the benchmarks with their protocols. The Kalman filter and the extended Kalman filter are assumed; factor graphs describe the optimization side in general.

Problem Definition

What an IMU measures

An IMU combines a three-axis gyroscope and a three-axis accelerometer. Forster et al. (2017, Eqs. 27–28) model the readings in the body frame $B$ of the IMU as

$$\tilde{\boldsymbol\omega}(t) = \boldsymbol\omega(t) + \mathbf{b}_g(t) + \boldsymbol\eta_g(t), \qquad \tilde{\mathbf{a}}(t) = R^\top(t)\big(\mathbf{a}(t) - \mathbf{g}\big) + \mathbf{b}_a(t) + \boldsymbol\eta_a(t),$$

where $R$ rotates body coordinates into the world frame, $\mathbf{a}$ is the body’s acceleration in the world frame, $\mathbf{g}$ is gravity, $\boldsymbol\eta$ is white noise and $\mathbf{b}_g$, $\mathbf{b}_a$ are slowly varying biases. Two features of this model drive everything that follows. First, the accelerometer measures acceleration minus gravity, so a resting IMU reads about $9.81\ \text{m/s}^2$ upward; gravity is not noise to be removed but a signal that reveals which way is down. Second, the biases are unknown and change: MSCKF and Forster et al. model them as random walks (integrated white noise), and OKVIS adds a slow decay of the accelerometer bias towards zero, a bounded random walk with a time constant (Leutenegger et al., 2015, Eq. 12), which its experiments set to infinity (Table I). MSCKF also models the Earth’s rotation, choosing an Earth-centred, Earth-fixed world frame, while Forster et al. and Leutenegger et al. neglect it.

The motion follows from integrating the readings: $\dot R = R\,\boldsymbol\omega^\wedge$, $\dot{\mathbf{v}} = \mathbf{a}$, $\dot{\mathbf{p}} = \mathbf{v}$ (Forster et al., 2017, Eq. 29), where $\boldsymbol\omega^\wedge$ is the skew-symmetric matrix of $\boldsymbol\omega$. It follows from the second integration that an uncorrected constant accelerometer bias $\mathbf{b}$ displaces the position by $\tfrac12 \mathbf{b}\,t^2$, and that a small tilt error $\delta\theta$ leaks gravity into the horizontal acceleration as about $9.81\,\delta\theta\ \text{m/s}^2$; this is why an IMU integrated alone drifts within seconds, and why its biases and attitude must be estimated rather than assumed.

The estimation problem

A VIO system estimates, at each time, the IMU orientation $R$, position $\mathbf{p}$ and velocity $\mathbf{v}$ in a gravity-aligned world frame, together with the biases $\mathbf{b}_g$ and $\mathbf{b}_a$. In MSCKF this is the 16-dimensional IMU state (a unit quaternion, two biases, velocity and position), with a 15-dimensional error state because orientation errors are represented by three angles (Mourikis and Roumeliotis, 2007, Eqs. 1–3). The camera adds observations: a static landmark seen from camera pose $(R_C, \mathbf{p}_C)$ projects to the image through the pinhole model, and the camera pose follows from the IMU pose through the camera–IMU extrinsic transformation $(R_{BC}, \mathbf{p}_{BC})$, which MSCKF treats as known (Eq. 14) and later systems can estimate online. The output is a camera pose trajectory with velocity and biases, available at the camera rate or, by propagating the latest estimate with new IMU samples, at the IMU rate, which the VINS-Mono technical report uses as high-frequency feedback for control.

Like visual odometry, VIO lets error accumulate; it does not recognize revisited places. Systems that add loop closure and a global map are visual-inertial SLAM.

Input and Output

• Input: synchronized, time-stamped images from one or more calibrated cameras; gyroscope and accelerometer samples, typically 100 Hz to 1 kHz according to Forster et al. (2017); the camera intrinsics, the camera–IMU extrinsics and the time offset between the two clocks, from calibration or estimated online; IMU noise and bias random-walk parameters.
• Output: the pose, velocity and bias estimates for each frame or keyframe, often with their covariance, and usually a sparse set of landmarks.

Variants

• Sensor configuration. Monocular-inertial, stereo-inertial and multi-camera-inertial. OKVIS is written for $N \ge 1$ cameras, and its authors note that a stereo rig behaves increasingly like a monocular one as the ratio of baseline to scene depth shrinks (Leutenegger et al., 2015).
• Coupling. Loosely coupled systems estimate the pose from vision alone and fuse it with the IMU in a separate step; tightly coupled systems put camera and IMU measurements into one estimator with all their correlations (Leutenegger et al., 2015). Leutenegger et al. (2013, 2015) report that the correlations matter for high-precision estimation; in the Delmerico and Scaramuzza (2018) benchmark, the only loosely coupled pipeline, SVO+MSF, was the most efficient and the least accurate.
• Estimator. Forster et al. (2017) classify estimators along three axes: how many past states they keep (filters keep the latest, fixed-lag smoothers a window, full smoothers the whole history), whether they represent uncertainty by a covariance or an information matrix, and how often they relinearize measurements (once in a standard EKF, repeatedly in a smoother). They note that the iterated EKF is equivalent to Gauss–Newton, so the families are closer than their names suggest.
• Front end. Most systems track sparse corners or descriptors; ROVIO, in the description of Delmerico and Scaramuzza (2018), is an EKF that tracks 3D landmarks with image patch features, and Forster et al. (2017) use the semi-direct SVO front end, which aligns image patches photometrically (see direct methods).
• Online calibration. The camera–IMU extrinsics and the time offset can be estimated with the motion. OKVIS’s journal version estimates the extrinsics online (Leutenegger et al., 2015), and the VINS-Mono README lists online extrinsic and temporal calibration. Kelly and Sukhatme (2011, not read) studied this self-calibration problem.

Challenges

• What the sensors cannot see. For generic motion, a camera–IMU system has four unobservable degrees of freedom: the three of global position and the rotation about gravity (yaw). Roll, pitch, metric scale, velocity and the biases become observable (Hesch et al., 2014; Martinelli, 2014). Observable is a statement about generic motion: Martinelli (2014) shows that the closed-form visual-inertial problem can have one, two or infinitely many solutions depending on the trajectory, the number and layout of point features and the number of images, that motion at constant acceleration generally gives two solutions even without bias, and that biased inertial data need more images and more restrictive conditions on the trajectory.
• Consistency. An estimator is consistent when its errors match its reported covariance. Hesch et al. (2014) show that, because of linearization errors, the linearized model used by a standard EKF has fewer unobservable directions than the nonlinear system; in Forster et al.’s (2017) summary, linearizing at the wrong estimate leaves three instead of four. The filter then acquires nonexistent information about yaw and becomes inconsistent. Their observability-constrained VINS (OC-VINS) modifies the state-transition and measurement Jacobians so that the linearized system keeps the four unobservable directions. Forster et al. (2017) summarize the remedies of that line of work, first-estimates Jacobians and observability-constrained filters, and check consistency of their own smoother with the normalized estimation error squared (NEES).
• Initialization. A monocular tightly coupled problem is highly nonlinear and has no metric scale until the IMU supplies it, so it needs good initial values for scale, gravity direction, velocity and biases (VINS-Mono technical report, Section III). The section on initialization below and the example show how this is done and where it fails.
• Rates. The IMU produces many measurements per frame. Adding a state per IMU sample to an optimization is impractically slow, and integrating between frames from the current state estimate has to be repeated each time that estimate changes (Forster et al., 2017, Section II). IMU preintegration addresses this.
• Timing and calibration. The fusion assumes that camera and IMU timestamps agree and that the extrinsics are right; the VINS-Mono README asks for accurate time stamps, images above 20 Hz and IMU data above 100 Hz, and ranks a rolling-shutter camera with an unsynchronized low-frequency IMU last (as vision for drones and UAVs describes).
• Computation. Accuracy can be bought with computation, which is scarce on small robots. Delmerico and Scaramuzza (2018) conclude that there is no free lunch: the most accurate and robust pipeline in their benchmark also used the most resources (Datasets and Benchmarks below).

Approaches

IMU preintegration

Between two keyframes $i$ and $j$, the IMU’s readings can be summarized once into a relative-motion measurement that does not depend on the state at $i$. Forster et al. (2017, Eqs. 33 and 35–37) define, for IMU samples $k$ spaced by $\Delta t$ and biases held constant over the interval,

$$\Delta \tilde R_{ij} = \prod_{k=i}^{j-1} \operatorname{Exp}\big((\tilde{\boldsymbol\omega}_k - \mathbf{b}_g)\Delta t\big), \qquad \Delta \tilde{\mathbf{v}}_{ij} = \sum_{k=i}^{j-1} \Delta \tilde R_{ik}\,(\tilde{\mathbf{a}}_k - \mathbf{b}_a)\,\Delta t,$$ $$\Delta \tilde{\mathbf{p}}_{ij} = \sum_{k=i}^{j-1} \Big[\Delta \tilde{\mathbf{v}}_{ik}\,\Delta t + \tfrac12 \Delta \tilde R_{ik}\,(\tilde{\mathbf{a}}_k - \mathbf{b}_a)\,\Delta t^2\Big],$$

where $\operatorname{Exp}$ maps a rotation vector to a rotation matrix (see rotation representations). Up to noise, these quantities equal functions of the two keyframe states (their Eq. 38):

$$\Delta \tilde R_{ij} \approx R_i^\top R_j, \qquad \Delta \tilde{\mathbf{v}}_{ij} \approx R_i^\top(\mathbf{v}_j - \mathbf{v}_i - \mathbf{g}\,\Delta t_{ij}), \qquad \Delta \tilde{\mathbf{p}}_{ij} \approx R_i^\top\big(\mathbf{p}_j - \mathbf{p}_i - \mathbf{v}_i\,\Delta t_{ij} - \tfrac12 \mathbf{g}\,\Delta t_{ij}^2\big).$$

The left-hand sides are computed from the IMU alone, so when the optimizer moves the states the measurements stay fixed. As Forster et al. stress, $\Delta\tilde{\mathbf{v}}_{ij}$ and $\Delta\tilde{\mathbf{p}}_{ij}$ are not physical changes of velocity and position: gravity and the initial velocity have been moved to the other side. The noise of the preintegrated measurement is, to first order, zero-mean Gaussian, and its covariance is propagated sample by sample (their Eqs. 42–43). The one remaining dependence is on the bias: when the bias estimate changes by $\delta\mathbf{b}$, the measurements are corrected to first order, for example $\Delta\tilde R_{ij}(\bar{\mathbf{b}}_g + \delta\mathbf{b}_g) \approx \Delta\tilde R_{ij}(\bar{\mathbf{b}}_g)\operatorname{Exp}\big(\tfrac{\partial \Delta\bar R_{ij}}{\partial \mathbf{b}_g}\delta\mathbf{b}_g\big)$, with Jacobians precomputed during integration (their Eq. 44). In their Monte Carlo test, integrating 100 random IMU measurements and perturbing both biases by 0.04 to 0.2, the error of this correction against re-integration was negligible (Fig. 11). The preintegrated measurement becomes one factor between keyframes, with residuals for rotation (through $\operatorname{Log}$), velocity and position.

The idea is older than the on-manifold form. Forster et al. build on what they call the seminal work of Lupton and Sukkarieh (2012, not read), who showed that the repeated integration can be avoided by reparametrizing the relative-motion constraints, which is what is called IMU preintegration; that work used Euler angles. They list three drawbacks of that parameterization that their formulation avoids: Euler-angle integration is exact only to first order, the resulting likelihood is not invariant to the choice of world frame, and noise propagation degrades near the gimbal-lock singularity, which their simulation shows as a growing divergence between propagated and sampled rotation covariance as the pitch approaches 90°. The RSS 2015 paper introduced the theory; the 2017 journal version adds the derivations and Jacobians in an appendix, consistency and bias experiments, and comparison with full batch estimation, and the authors released the preintegrated IMU factor in GTSAM 4.0.

VINS-Mono uses its own version of preintegration with quaternions and continuous-time noise propagation (technical report, Section II-B), and its code integrates with the mid-point rule (integration_base.h). When the bias estimate moves far from the value used for integration, more than 0.10 m/s² for the accelerometer or 0.01 rad/s for the gyroscope, the code re-integrates instead of relying on the first-order correction (imu_factor.h).

Filtering: the multi-state constraint Kalman filter

In the review of Mourikis and Roumeliotis (2007), one family of visual-inertial filters follows the EKF-SLAM pattern and keeps the 3D landmarks in the state, which accounts for the pose–landmark correlations but is costly when thousands of features are seen. The multi-state constraint Kalman filter keeps no landmarks. It propagates the IMU state with each IMU sample, using fifth-order Runge–Kutta integration and the error-state covariance of the extended Kalman filter, and, each time an image arrives, appends a copy of the current camera pose to the state, up to a maximum of $N_{\max}$ poses.

The measurement model is the paper’s contribution. A feature $j$ seen from $M_j$ camera poses gives $2M_j$ image coordinates, whose linearized residual depends on the pose errors and on the error of the feature’s position, estimated beforehand by Gauss–Newton triangulation in inverse-depth form:

$$\mathbf{r}^{(j)} \approx H_X^{(j)}\,\tilde{X} + H_f^{(j)}\,{}^G\tilde{\mathbf{p}}_{f_j} + \mathbf{n}^{(j)}.$$

Because the feature estimate was computed from the state, its error is correlated with $\tilde{X}$ and the residual cannot be used directly. MSCKF projects it onto the left null space of $H_f^{(j)}$, a matrix $A$ with $A^\top H_f^{(j)} = 0$:

$$\mathbf{r}_o^{(j)} = A^\top \mathbf{r}^{(j)} \approx A^\top H_X^{(j)}\,\tilde{X} + A^\top \mathbf{n}^{(j)},$$

which removes the feature and leaves $2M_j - 3$ constraints among all the poses that saw it (Mourikis and Roumeliotis, 2007, Eqs. 22–24). The authors call the update optimal up to linearization error, because these constraints carry all the information the measurements give about the poses. Updates happen when a feature leaves the field of view, or when the pose window is full: then $N_{\max}/3$ poses, evenly spaced and starting from the second oldest, are removed after using their features, while the oldest pose is kept because its longer baseline carries more information. A QR decomposition shrinks the stacked residual before the Kalman update, and the cost is linear in the number of features and at most cubic in the number of poses in the state. Delayed linearization (each feature is processed after all its observations) and the inverse-depth parameterization make the filter more robust to linearization error.

The ICRA paper’s experiment drove a car with a 640 × 480 camera at 3 Hz and an IMU at 100 Hz for about 9 minutes and 3.2 km through a residential area, with up to 30 poses in the state. Without loop closure or motion constraints, the final position error was about 10 m, which the authors give as 0.31% of the distance travelled; there was no GPS, so this figure comes from the known initial and final parking spots, and the rest of the trajectory was judged against a street map. The filter processed the recorded data at 14 Hz on one core of an Intel T7200 (2 GHz).

Hesch et al. (2014) applied their observability constraints to MSCKF. ROVIO, another EKF, appears in the benchmark below. MSCKF’s measurement model is also the structureless idea that Forster et al. (2017) carried into smoothing.

Optimization: OKVIS, preintegrated smoothing and VINS-Mono

Optimization-based estimators write VIO as nonlinear least squares over a set of states, as in bundle adjustment, and relinearize every measurement at each iteration; the factor graph article describes the machinery and the marginalization that keeps a window bounded.

OKVIS. Leutenegger et al. (2013, RSS; 2015, IJRR, read in the authors’ revised manuscript) minimize a cost that combines weighted reprojection errors of landmarks with one IMU error term between consecutive frames. Each IMU term predicts the next state from the previous one by Runge–Kutta integration of all IMU samples in between and compares it with the estimated state, weighted by the propagated covariance (journal version, Eqs. 16–17). The window holds $S$ recent frames and $M$ keyframes that may lie far in the past, so a stationary camera does not lose its visual constraints; older states and landmarks are marginalized into a linear prior. The journal version used $M = 7$ and $S = 3$ in all its datasets, added a monocular variant and online extrinsic estimation, and compared it with a reference MSCKF implementation that kept five pose clones. The authors report that their method consistently gave lower errors at higher computational cost. Forster et al. (2017) point out that OKVIS repeats the IMU integration whenever the linearization point changes, which preintegration avoids.

Preintegrated full smoothing. Forster et al. (2017) build their estimator as a factor graph of preintegrated IMU factors and structureless vision factors that, as in MSCKF, eliminate the 3D points, and they solve it with the incremental smoother iSAM2 without marginalizing past states. On a 430 m indoor sequence recorded with an ADIS16448 IMU at 800 Hz and a 20 Hz camera, with Vicon ground truth, they report an average drift of 0.3 m over 360 m segments, against 0.7 m for both OKVIS and MSCKF; they caution that the three systems use different front ends, so the comparison is of whole systems.

VINS-Mono. The T-RO paper (Qin, Li and Shen, 2018) was not read for this article; the description follows the technical report in the official repository, which lists a fourth author, Zhenfei Yang, and the code at commit 90dabb5. (The repository’s BibTeX entry gives the T-RO paper three authors, while the README’s list of papers also names Yang.) The front end tracks 100 to 300 corners per image with KLT and rejects outliers with a fundamental-matrix RANSAC test. A frame becomes a keyframe when the average parallax of its tracked features, after compensating for rotation with the gyroscope, exceeds a threshold, or when too few features remain tracked. The back end optimizes, over a sliding window of keyframes (WINDOW_SIZE is 10 in the code), the positions, velocities, orientations and biases of the frames, the camera–IMU extrinsics and the inverse depths of the features, by minimizing with Ceres Solver the sum of a prior from marginalization, the preintegrated IMU residuals and visual residuals defined on the unit sphere (report, Section IV and Eq. 20). If the second-newest frame is a keyframe, the oldest keyframe is marginalized into the prior; otherwise the second-newest frame’s visual measurements are dropped and its IMU measurements kept, which preserves sparsity. A separate thread detects loops with DBoW2 and BRIEF descriptors and optimizes a pose graph in only four degrees of freedom, position and yaw, because roll and pitch are observable from gravity; loop edges use a Huber cost against false detections.

How the families compare. The sources describe a trade-off rather than a winner. Forster et al. (2017) write that filters are fast but lose accuracy as linearization errors accumulate, that fixed-lag smoothers inherit some of the filters’ consistency problems through marginalization, and that full smoothing is the most accurate but grows with the trajectory unless incremental methods are used. The VINS-Mono technical report compares its odometry with OKVIS on EuRoC MH03 and MH05 and calls the two “hard to distinguish” in pure odometry accuracy, while noting that OKVIS estimated roll and pitch more stably on MH05; its claimed advantages are completeness: initialization, failure recovery and loop closure. Loop closure and pose graph optimization, which turn a VIO system into visual-inertial SLAM, go beyond odometry.

Initialization

Every estimator needs initial values for attitude, velocity, biases and, with one camera, scale. The VINS-Mono technical report describes the usual practice as assuming the device is stationary and taking the average of the first seconds of accelerometer data as gravity, and criticizes it: the assumption breaks with non-trivial bias or when the device accelerates at the start. Its procedure is loosely coupled (Section III):

1. Visual structure from motion in the window. When a previous frame has more than 30 tracked features and more than 20 pixels of average parallax with the newest frame (the code at commit 90dabb5 requires more than 20 correspondences and tests the parallax against 30 pixels at a focal length of 460), the five-point algorithm gives their relative pose up to scale, the common features are triangulated, perspective-n-point poses the other frames, and bundle adjustment refines everything. The result is a set of camera poses $(R^{c_0}_{b_k}, \bar{\mathbf{p}}^{c_0}_{c_k})$ in the first camera frame $c_0$, with translations in unknown units.

2. Gyroscope bias. For consecutive frames, the visual relative rotation is compared with the preintegrated rotation, linearized in $\mathbf{b}_g$; a linear least-squares solve gives the bias, and the IMU data are re-integrated with it (Eq. 14). The code does this as one linear solve (solveGyroscopeBias).

3. Velocity, gravity and scale. With rotations known, the preintegrated position and velocity relations are linear in the unknowns $X_I = [\mathbf{v}_0, \dots, \mathbf{v}_n, \mathbf{g}^{c_0}, s]$: one velocity per frame in its body frame, gravity in the $c_0$ frame and the scale $s$ that maps the visual translations to metres (Eqs. 15–18). In the notation of this article, with gravity pointing down and frames $k$, $k+1$ separated by $\Delta t_k$,

   $$\Delta\tilde{\mathbf{p}}_{k,k+1} = R_k^\top s\,(\bar{\mathbf{p}}_{k+1} - \bar{\mathbf{p}}_k) - \mathbf{v}_k\,\Delta t_k - \tfrac12 R_k^\top \mathbf{g}\,\Delta t_k^2, \qquad \Delta\tilde{\mathbf{v}}_{k,k+1} = R_k^\top R_{k+1}\mathbf{v}_{k+1} - \mathbf{v}_k - R_k^\top \mathbf{g}\,\Delta t_k,$$

   where $R_k$ is the body-to-$c_0$ rotation and the velocities are in body coordinates.

4. Gravity refinement. Since $|\mathbf{g}|$ is known, gravity has two degrees of freedom; it is re-parameterized on the tangent plane of the sphere around the current estimate, and the linear solve is repeated until it converges (four iterations in the code). Rotating the refined gravity onto the world $z$ axis fixes the world frame.

The accelerometer bias is not estimated in this procedure: the code re-integrates with a zero accelerometer bias and leaves $\mathbf{b}_a$ to the subsequent nonlinear optimization. The code rejects the solution if the norm of the first gravity estimate differs from the configured gravity magnitude (9.81007 m/s² in the EuRoC configuration) by more than 1 m/s², or if the scale is negative. It also computes the spread of the preintegrated $\Delta\tilde{\mathbf{v}}/\Delta t$ over the window as an “IMU excitation” check against 0.25 m/s², but at commit 90dabb5 the line that would abort on low excitation is commented out, so only a message is logged.

Martinelli (2014) derived a closed-form solution of the same visual-inertial structure-from-motion problem, motivated by this initialization issue, and used it to analyse when the problem has a finite number of solutions, including the biased case.

Example: Scale, Gravity and Bias at Initialization

The script follows the VINS-Mono initialization above on synthetic data and asks how much the scale and gravity depend on estimating the biases. A camera–IMU rig (camera frame equal to the IMU frame, an identity extrinsic, f = 460 px, 752 × 480 images) moves past 300 random points while the IMU runs at 200 Hz and frames arrive at 10 Hz. The true trajectory is produced by integrating chosen angular rates and accelerations with the same discrete scheme the preintegration uses, so noise-free preintegration is exact. Gravity is $(0, 0, -9.81)$ in a $z$-up world frame, rotations map body to world coordinates, and the camera axes are $x$ right, $y$ down and $z$ forward. The IMU has constant biases $\mathbf{b}_g = (0.02, -0.01, 0.015)$ rad/s and $\mathbf{b}_a = (0.10, -0.08, 0.12)$ m/s², values chosen for the example; the noisy runs add per-sample noise of 0.08 m/s² and 0.004 rad/s, the acc_n and gyr_n values of VINS-Mono’s EuRoC configuration (which the file labels standard deviations), and 0.5 px of pixel noise.

Each window is first reconstructed from its image points by a small monocular SfM: an eight-point essential matrix between the first and last frame, triangulation, PnP and bundle adjustment, with the first-to-last camera distance fixed to 1, so the true scale is that distance in metres. Three estimators then align the IMU to it:

• E0 assumes zero biases.
• E1 estimates the gyroscope bias with one linear solve and then velocities, gravity and scale, as VINS-Mono does.
• E2 adds the accelerometer bias as three more unknowns of the same linear system; this extension is ours, possible because the preintegrated velocity and position are linear in $\mathbf{b}_a$.

Two motions are compared, a handheld-like “exciting” one (±0.4 m sinusoids and ±15° attitude oscillations on top of 0.5 m/s forward motion) and a “gentle” one with accelerations and rotations ten times smaller, each over a 1 s window of 11 frames (VINS-Mono’s window size plus one) and a 3 s window of 31 frames. The script prints the scale error, the bias errors and the angle between estimated and true gravity, first on noise-free data and then over noise seeds 0–9, and counts the runs that pass VINS-Mono’s acceptance test on the first gravity estimate and the scale.

"""VIO initialization example: scale, gravity and IMU biases from a synthetic monocular camera and IMU."""
import numpy as np
from scipy.optimize import least_squares
from scipy.sparse import lil_matrix
from scipy.spatial.transform import Rotation

Exp = lambda w: Rotation.from_rotvec(w).as_matrix()
Log = lambda R: Rotation.from_matrix(R).as_rotvec()
G = np.array([0.0, 0.0, -9.81])           # gravity in the world frame (z up)
F_PX, W, H = 460.0, 752, 480              # pinhole camera; camera frame = IMU (body) frame
IMU_HZ, CAM_HZ = 200, 10
R_BODY0 = np.array([[0, 0, 1], [-1, 0, 0], [0, -1, 0]], float)     # columns: body x,y,z in world


def simulate(motion, T, bg, ba, sig_a, sig_g, sig_px, seed):
    """Discrete-time truth (same scheme as the preintegration), IMU samples and image points."""
    rng = np.random.default_rng(seed)
    dt, n = 1.0 / IMU_HZ, int(round(T * IMU_HZ))
    k = 1.0 if motion == "exciting" else 0.1
    A, f, ph = 0.4 * k * np.ones(3), np.array([0.6, 0.75, 0.9]), np.array([0.3, 1.1, 2.0])
    rot_amp, rf = np.radians(15) * k, np.array([0.7, 0.5, 0.8])
    t = np.arange(n + 1) * dt
    acc_w = -(A * (2 * np.pi * f) ** 2) * np.sin(2 * np.pi * f * t[:, None] + ph)
    omg_b = rot_amp * 2 * np.pi * rf * np.cos(2 * np.pi * rf * t[:, None])
    R, v, p = [R_BODY0], [np.array([0.5, 0, 0]) + A * 2 * np.pi * f * np.cos(ph)], [np.zeros(3)]
    for i in range(n):
        R.append(R[-1] @ Exp(omg_b[i] * dt))
        p.append(p[-1] + v[-1] * dt + 0.5 * acc_w[i] * dt ** 2)
        v.append(v[-1] + acc_w[i] * dt)
    acc_m = np.array([R[i].T @ (acc_w[i] - G) for i in range(n)]) + ba + sig_a * rng.normal(size=(n, 3))
    gyr_m = omg_b[:n] + bg + sig_g * rng.normal(size=(n, 3))
    cam_idx = np.arange(0, n + 1, IMU_HZ // CAM_HZ)
    pts = np.column_stack([rng.uniform(3, 8, 300), rng.uniform(-3, 3, 300), rng.uniform(-2, 2, 300)])
    obs = []                                        # per frame: (point ids, pixel coords)
    for c in cam_idx:
        Xc = (pts - p[c]) @ R[c]                    # world -> camera
        uv = F_PX * Xc[:, :2] / Xc[:, 2:3] + [W / 2, H / 2]
        vis = (Xc[:, 2] > 0.3) & (uv[:, 0] > 0) & (uv[:, 0] < W) & (uv[:, 1] > 0) & (uv[:, 1] < H)
        ids = np.flatnonzero(vis)
        obs.append((ids, uv[ids] + sig_px * rng.normal(size=(len(ids), 2))))
    truth = dict(R=[R[c] for c in cam_idx], p=np.array([p[c] for c in cam_idx]),
                 v=np.array([v[c] for c in cam_idx]))
    return truth, acc_m, gyr_m, cam_idx, obs


def preintegrate(acc, gyr, bg, ba):
    """Forster et al. (2017) Eq. 33: dR, dv, dp between two keyframes, in the first body frame."""
    dt = 1.0 / IMU_HZ
    dR, dv, dp = np.eye(3), np.zeros(3), np.zeros(3)
    for a, w in zip(acc, gyr):
        dp = dp + dv * dt + 0.5 * dR @ (a - ba) * dt ** 2
        dv = dv + dR @ (a - ba) * dt
        dR = dR @ Exp((w - bg) * dt)
    return dR, dv, dp, len(acc) * dt


# ---------------- monocular SfM up to scale ----------------
def to_bearing(uv):
    return np.column_stack([(uv - [W / 2, H / 2]) / F_PX, np.ones(len(uv))])


def triangulate(Ps, xs):
    A = np.vstack([np.vstack([x[0] * P[2] - P[0], x[1] * P[2] - P[1]]) for P, x in zip(Ps, xs)])
    X = np.linalg.svd(A)[2][-1]
    return X[:3] / X[3]


def essential_pose(x1, x2):
    """Normalized 8-point essential matrix, decomposed with the cheirality check."""
    A = np.column_stack([x2[:, [0]] * x1, x2[:, [1]] * x1, x1])
    E = np.linalg.svd(A)[2][-1].reshape(3, 3)
    U, _, Vt = np.linalg.svd(E)
    if np.linalg.det(U) < 0: U = -U
    if np.linalg.det(Vt) < 0: Vt = -Vt
    Wm = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]])
    best = None
    for R in (U @ Wm @ Vt, U @ Wm.T @ Vt):
        for t in (U[:, 2], -U[:, 2]):
            P1, P2 = np.hstack([np.eye(3), np.zeros((3, 1))]), np.hstack([R, t[:, None]])
            Xs = np.array([triangulate([P1, P2], [a, b]) for a, b in zip(x1[:30], x2[:30])])
            good = np.sum((Xs[:, 2] > 0) & ((Xs @ R.T + t)[:, 2] > 0))
            if best is None or good > best[0]:
                best = (good, R, t)
    return best[1], best[2]


def project(R, t, X):
    Y = X @ R.T + t
    return Y[:, :2] / Y[:, 2:3]


def pnp(X, x, R0, t0):
    R, t = R0, t0
    for _ in range(10):
        r = (project(R, t, X) - x).ravel()
        J = np.zeros((len(r), 6))
        for m in range(6):
            d = np.zeros(6); d[m] = 1e-7
            J[:, m] = ((project(Exp(d[:3]) @ R, t + d[3:], X) - x).ravel() - r) / 1e-7
        dx = np.linalg.lstsq(J, -r, rcond=None)[0]
        R, t = Exp(dx[:3]) @ R, t + dx[3:]
    return R, t


def sfm(obs):
    """Up-to-scale poses (camera-to-c0 rotation, camera centre in c0) of all frames."""
    nF = len(obs)
    b = [dict(zip(ids, to_bearing(uv))) for ids, uv in obs]
    common = sorted(set(b[0]) & set(b[-1]))
    Rl, tl = essential_pose(np.array([b[0][i] for i in common]), np.array([b[-1][i] for i in common]))
    Rs, ts = [None] * nF, [None] * nF
    Rs[0], ts[0], Rs[-1], ts[-1] = np.eye(3), np.zeros(3), Rl, tl
    P = lambda k: np.hstack([Rs[k], ts[k][:, None]])
    X = {i: triangulate([P(0), P(nF - 1)], [b[0][i], b[-1][i]]) for i in common}
    for k in range(1, nF - 1):
        ids = [i for i in b[k] if i in X]
        Rs[k], ts[k] = pnp(np.array([X[i] for i in ids]), np.array([b[k][i][:2] for i in ids]),
                           Rs[k - 1], ts[k - 1])
    for i in set().union(*[set(d) for d in b]) - set(X):
        ks = [k for k in range(nF) if i in b[k]]
        if len(ks) >= 2:
            X[i] = triangulate([P(k) for k in ks], [b[k][i] for k in ks])
    pid = sorted(X)
    col = {i: j for j, i in enumerate(pid)}
    meas = [(k, col[i], b[k][i][:2]) for k in range(nF) for i in b[k] if i in col]
    x0 = np.concatenate([np.concatenate([Log(Rs[k]), ts[k]]) for k in range(1, nF)] + [X[i] for i in pid])
    nP = 6 * (nF - 1)

    mk = np.array([m[0] for m in meas]); mj = np.array([m[1] for m in meas])
    mz = np.array([m[2] for m in meas])

    def resid(x):
        poses = x[:nP].reshape(-1, 6); pts = x[nP:].reshape(-1, 3)
        Rk = np.concatenate([np.eye(3)[None], Rotation.from_rotvec(poses[:, :3]).as_matrix()])
        tk = np.vstack([np.zeros(3), poses[:, 3:]])
        Y = np.einsum("mab,mb->ma", Rk[mk], pts[mj]) + tk[mk]
        return (Y[:, :2] / Y[:, 2:3] - mz).ravel() * F_PX

    S = lil_matrix((2 * len(meas), len(x0)), dtype=int)
    for m, (k, j, _) in enumerate(meas):
        if k > 0:
            S[2 * m:2 * m + 2, 6 * (k - 1):6 * k] = 1
        S[2 * m:2 * m + 2, nP + 3 * j:nP + 3 * j + 3] = 1
    x = least_squares(resid, x0, jac_sparsity=S, x_scale="jac", method="trf").x
    poses = x[:nP].reshape(-1, 6)
    Rk = [np.eye(3)] + [Exp(q[:3]) for q in poses]
    tk = [np.zeros(3)] + [q[3:] for q in poses]
    R_c0 = [R.T for R in Rk]                                  # camera-to-c0 rotations
    c = np.array([-R.T @ t for R, t in zip(Rk, tk)])         # camera centres in c0
    return R_c0, c / np.linalg.norm(c[-1] - c[0])            # first-to-last distance = 1


# ---------------- visual-inertial alignment ----------------
def num_jac(f, x0, eps=1e-6):
    f0 = f(x0)
    return np.column_stack([(f(x0 + eps * e) - f0) / eps for e in np.eye(len(x0))])


def initialize(R_c0, cbar, acc, gyr, cam_idx, est_bg, est_ba):
    n = len(R_c0)
    seg = [(cam_idx[k], cam_idx[k + 1]) for k in range(n - 1)]
    bg = np.zeros(3)
    if est_bg:                                  # one linearized solve, as solveGyroscopeBias
        A, rhs = np.zeros((3, 3)), np.zeros(3)
        for k, (i0, i1) in enumerate(seg):
            dR0 = preintegrate(acc[i0:i1], gyr[i0:i1], bg, np.zeros(3))[0]
            J = num_jac(lambda b: Log(dR0.T @ preintegrate(acc[i0:i1], gyr[i0:i1], b, np.zeros(3))[0]), bg)
            r = Log(dR0.T @ R_c0[k].T @ R_c0[k + 1])
            A += J.T @ J; rhs += J.T @ r
        bg = bg + np.linalg.solve(A, rhs)
    pre, Jba = [], []
    for i0, i1 in seg:
        dR, dv, dp, T = preintegrate(acc[i0:i1], gyr[i0:i1], bg, np.zeros(3))
        pre.append((dv, dp, T))
        Jba.append(num_jac(lambda b: np.concatenate(preintegrate(acc[i0:i1], gyr[i0:i1], bg, b)[1:3]),
                           np.zeros(3)))

    def solve(gdir=None):
        """Unknowns: body-frame velocities, gravity (3, or 2 tangent coords), scale, [b_a]."""
        ng = 3 if gdir is None else 2
        if gdir is not None:
            B = np.linalg.svd(gdir[None, :])[2][1:].T          # tangent basis of the gravity sphere
        nx = 3 * n + ng + 1 + (3 if est_ba else 0)
        rows, rhs = [], []
        for k, (dv, dp, T) in enumerate(pre):
            Ri, Rj = R_c0[k], R_c0[k + 1]
            Ap, Av = np.zeros((3, nx)), np.zeros((3, nx))
            Ap[:, 3 * k:3 * k + 3] = -T * np.eye(3)
            Av[:, 3 * k:3 * k + 3] = -np.eye(3)
            Av[:, 3 * k + 3:3 * k + 6] = Ri.T @ Rj
            gp, gv = -0.5 * T ** 2 * Ri.T, -T * Ri.T
            bp, bv = dp.copy(), dv.copy()
            if gdir is None:
                Ap[:, 3 * n:3 * n + 3], Av[:, 3 * n:3 * n + 3] = gp, gv
            else:
                g0 = 9.81 * gdir
                Ap[:, 3 * n:3 * n + 2], Av[:, 3 * n:3 * n + 2] = gp @ B, gv @ B
                bp, bv = bp - gp @ g0, bv - gv @ g0
            Ap[:, 3 * n + ng] = Ri.T @ (cbar[k + 1] - cbar[k])
            if est_ba:                          # dp(b_a) = dp(0) + J b_a (exact: linear in b_a)
                Ap[:, -3:], Av[:, -3:] = -Jba[k][3:], -Jba[k][:3]
            rows += [Ap, Av]; rhs += [bp, bv]
        x = np.linalg.lstsq(np.vstack(rows), np.concatenate(rhs), rcond=None)[0]
        return x, (B if gdir is not None else None)

    x, _ = solve()
    g_lin, s_lin = x[3 * n:3 * n + 3], x[3 * n + 3]
    vins_ok = abs(np.linalg.norm(g_lin) - 9.81) <= 1.0 and s_lin > 0   # VINS-Mono's acceptance check
    gdir = g_lin / np.linalg.norm(g_lin)
    for _ in range(4):                          # gravity refinement on the sphere, 4 iterations
        x, B = solve(gdir)
        g = 9.81 * gdir + B @ x[3 * n:3 * n + 2]
        gdir = g / np.linalg.norm(g)
    s = x[3 * n + 2]
    ba = x[-3:] if est_ba else np.zeros(3)
    return dict(s=s, g=9.81 * gdir, bg=bg, ba=ba, vins_ok=vins_ok)


def evaluate(est, truth, bg, ba):
    s_true = np.linalg.norm(truth["p"][-1] - truth["p"][0])
    g_true = truth["R"][0].T @ G
    cosang = est["g"] @ g_true / (np.linalg.norm(est["g"]) * 9.81)
    return (100 * (est["s"] / s_true - 1), np.linalg.norm(est["bg"] - bg),
            np.linalg.norm(est["ba"] - ba), np.degrees(np.arccos(np.clip(cosang, -1, 1))))


ESTIMATORS = {"E0 no bias": (False, False), "E1 gyro bias": (True, False), "E2 gyro+acc bias": (True, True)}
BG, BA = np.array([0.02, -0.01, 0.015]), np.array([0.10, -0.08, 0.12])

if __name__ == "__main__":
    SIG_A, SIG_G, SIG_PX = 0.08, 0.004, 0.5       # VINS-Mono's EuRoC acc_n, gyr_n; pixel noise
    print("noise-free measurements, both biases: scale error % / accel-bias error m/s^2 / gravity error deg")
    for motion in ("exciting", "gentle"):
        for T in (1.0, 3.0):
            truth, acc, gyr, cam_idx, obs = simulate(motion, T, BG, BA, 0, 0, 0, seed=0)
            R_c0, cbar = sfm(obs)
            cells = []
            for eg, ea in ESTIMATORS.values():
                e = evaluate(initialize(R_c0, cbar, acc, gyr, cam_idx, eg, ea), truth, BG, BA)
                cells.append(f"{e[0]:+6.2f} / {e[2]:.3f} / {e[3]:.2f}")
            print(f"  {motion:8s} {T:.0f} s: " + " | ".join(f"{n.split()[0]} {c}" for n, c in zip(ESTIMATORS, cells)))

    print("\nnoisy measurements, seeds 0-9: median [min, max] over seeds")
    for motion in ("exciting", "gentle"):
        for T in (1.0, 3.0):
            res = {n: [] for n in ESTIMATORS}
            for seed in range(10):
                truth, acc, gyr, cam_idx, obs = simulate(motion, T, BG, BA, SIG_A, SIG_G, SIG_PX, seed)
                R_c0, cbar = sfm(obs)
                for n, (eg, ea) in ESTIMATORS.items():
                    est = initialize(R_c0, cbar, acc, gyr, cam_idx, eg, ea)
                    res[n].append(evaluate(est, truth, BG, BA) + (est["vins_ok"],))
            print(f"{motion}, {T:.0f} s window")
            for n, r in res.items():
                r = np.array(r)
                f = lambda c, fmt: f"{np.median(r[:, c]):{fmt}} [{r[:, c].min():{fmt}}, {r[:, c].max():{fmt}}]"
                print(f"  {n:17s} scale % {f(0, '+.2f')}  bg {f(1, '.4f')}  ba {f(2, '.3f')}  "
                      f"g deg {f(3, '.2f')}  check passed {int(r[:, 4].sum())}/10")

Output:

noise-free measurements, both biases: scale error % / accel-bias error m/s^2 / gravity error deg
  exciting 1 s: E0  +0.84 / 0.175 / 1.03 | E1  +0.87 / 0.175 / 1.07 | E2  -0.00 / 0.000 / 0.00
  exciting 3 s: E0  +0.24 / 0.175 / 0.93 | E1  +0.23 / 0.175 / 0.95 | E2  +0.00 / 0.000 / 0.00
  gentle   1 s: E0  +8.66 / 0.175 / 1.06 | E1  +8.68 / 0.175 / 1.08 | E2  -0.00 / 0.000 / 0.00
  gentle   3 s: E0  +2.51 / 0.175 / 0.90 | E1  +2.51 / 0.175 / 0.91 | E2  -0.00 / 0.000 / 0.00

noisy measurements, seeds 0-9: median [min, max] over seeds
exciting, 1 s window
  E0 no bias        scale % +0.70 [+0.10, +1.61]  bg 0.0269 [0.0269, 0.0269]  ba 0.175 [0.175, 0.175]  g deg 1.02 [0.97, 1.16]  check passed 10/10
  E1 gyro bias      scale % +0.74 [+0.13, +1.64]  bg 0.0009 [0.0004, 0.0018]  ba 0.175 [0.175, 0.175]  g deg 1.06 [1.00, 1.20]  check passed 10/10
  E2 gyro+acc bias  scale % +0.01 [-0.97, +0.88]  bg 0.0009 [0.0004, 0.0018]  ba 0.441 [0.085, 1.005]  g deg 2.62 [0.17, 6.00]  check passed 10/10
exciting, 3 s window
  E0 no bias        scale % +0.26 [-0.23, +0.42]  bg 0.0269 [0.0269, 0.0269]  ba 0.175 [0.175, 0.175]  g deg 0.91 [0.86, 0.97]  check passed 10/10
  E1 gyro bias      scale % +0.25 [-0.24, +0.41]  bg 0.0003 [0.0001, 0.0004]  ba 0.175 [0.175, 0.175]  g deg 0.93 [0.87, 0.99]  check passed 10/10
  E2 gyro+acc bias  scale % +0.02 [-0.46, +0.17]  bg 0.0003 [0.0001, 0.0004]  ba 0.035 [0.016, 0.073]  g deg 0.19 [0.08, 0.45]  check passed 10/10
gentle, 1 s window
  E0 no bias        scale % +11.69 [+4.14, +24.67]  bg 0.0269 [0.0269, 0.0269]  ba 0.175 [0.175, 0.175]  g deg 1.15 [0.96, 1.45]  check passed 10/10
  E1 gyro bias      scale % +11.71 [+4.16, +24.70]  bg 0.0009 [0.0002, 0.0014]  ba 0.175 [0.175, 0.175]  g deg 1.17 [0.99, 1.48]  check passed 10/10
  E2 gyro+acc bias  scale % -1.80 [-12.26, +9.25]  bg 0.0009 [0.0002, 0.0014]  ba 5.754 [3.143, 9.271]  g deg 33.13 [18.61, 56.45]  check passed 5/10
gentle, 3 s window
  E0 no bias        scale % +1.68 [+0.92, +3.21]  bg 0.0269 [0.0269, 0.0269]  ba 0.175 [0.175, 0.175]  g deg 0.88 [0.84, 0.92]  check passed 10/10
  E1 gyro bias      scale % +1.67 [+0.91, +3.20]  bg 0.0003 [0.0001, 0.0004]  ba 0.175 [0.175, 0.175]  g deg 0.90 [0.85, 0.93]  check passed 10/10
  E2 gyro+acc bias  scale % -0.74 [-1.72, +0.47]  bg 0.0003 [0.0001, 0.0004]  ba 0.344 [0.067, 0.720]  g deg 2.01 [0.33, 4.24]  check passed 10/10

In the noise-free rows the scale and gravity errors of E0 and E1 come from the accelerometer bias alone: it is about 1.8% of gravity, and without an estimate of it the alignment tilts gravity by about a degree and stretches the scale by 0.2% to 8.7%, most on the gentle motion and the short window; in our reading, because there the bias is large compared with the motion’s own accelerations. E2, whose model contains the bias, recovers everything exactly. E0 and E1 give nearly the same scale throughout, so in this set-up the gyroscope bias, which E1 recovers to within 0.002 rad/s in every run, matters little for scale; the accelerometer bias drives it.

With noise, estimating the accelerometer bias pays off only when the motion constrains it. On the exciting 3 s window E2 brings the median scale error from about 0.25% to 0.02%, gravity within 0.5° and the bias within 0.08 m/s² in every run. On the exciting 1 s window and the gentle 3 s window its median scale error is still the smallest, but its median bias error is larger than that of assuming zero, and its median gravity error is larger than E1’s. In our reading, the three extra unknowns absorb noise that the short or gentle motion cannot separate from a tilt of gravity. On the gentle 1 s window nothing works. E0 and E1 overestimate the scale by 4% to 25%; E2’s gravity estimates are 19° to 56° wrong, with accelerometer-bias errors of 3 to 9 m/s², and VINS-Mono’s own acceptance test on the first gravity estimate rejects half of them while passing every E0 and E1 solution, whose scale errors it cannot see. This agrees with Martinelli’s (2014) finding that biased inertial data need more images and more restrictive conditions on the trajectory, and, in our reading, with VINS-Mono’s choice to leave the accelerometer bias out of its linear initialization and to estimate it later in the sliding-window optimization. The example simplifies the real problem: the extrinsics are known and the camera coincides with the IMU, the biases are constant, correspondences are perfect and the visual structure comes from an exact pinhole model.

The script takes about 7 to 8 minutes on our machine. It was written for NumPy 2.3.5 and SciPy 1.16.3, and two runs gave identical output. After the first run, the bundle-adjustment residual was vectorized for speed with the same mathematics; this changed one printed line (gentle motion, 1 s window, E2) in its last digits through floating-point differences in the optimization, and the output above is that of the vectorized script.

Datasets and Benchmarks

• EuRoC MAV: eleven indoor flights of a hexacopter with 20 Hz global-shutter stereo images and a 200 Hz IMU, with laser-tracker or Vicon ground truth; the usual VIO test set. It has no evaluation server, so every paper chooses its own alignment and subset; the EuRoC article describes the data, the ground-truth limits and the published conventions.
• The Delmerico and Scaramuzza (2018) benchmark ran six public monocular VIO pipelines, MSCKF, OKVIS, ROVIO, VINS-Mono (with and without loop closure), SVO+MSF and SVO+GTSAM, on all EuRoC sequences, on a laptop, an Intel NUC and two single-board computers, with parameters recommended by each method’s authors and kept fixed across trials, except for reductions needed to run in real time on the small boards. Each trajectory was aligned to the ground truth with a Sim(3) transformation (Umeyama’s method) before the RMS position error was computed. Their conclusions: VINS-Mono was the most consistently accurate and robust across platforms, SVO+GTSAM gave the most accurate trajectories in many combinations among methods without loop closure, OKVIS was accurate on all platforms but slow on the embedded ones, ROVIO was a good compromise with tightly bounded resource use, although on the UP Board no parameter setting made its filter converge, and SVO+MSF was the most efficient and least accurate. They did not search parameters per sequence.
• Later datasets. The handheld TUM VI benchmark was introduced partly to address EuRoC’s limits, as the EuRoC MAV article describes, and the racing-drone UZH-FPV dataset adds much faster flight, as vision for drones and UAVs summarizes.
• Driving data. KITTI was recorded with a GPS/IMU navigation system on the car, which also supplies its odometry ground truth, as the KITTI article describes.

The protocols in the sources differ, so numbers rarely compare across papers. The VINS-Mono technical report discards the first 100 outputs and aligns on the next 150 because a monocular system needs time to converge; Forster et al. (2017) align only the first frame and report relative errors over segments of 10 to 360 m; Delmerico and Scaramuzza align the whole trajectory in Sim(3).

Evaluation Metrics

• Absolute trajectory error after alignment. Because a VIO estimate has four unobservable degrees of freedom, the observability argument calls for aligning with a yaw rotation and a translation only; an SE(3) or Sim(3) alignment would also remove real tilt and scale errors, as the EuRoC article’s example shows. The alignment and the poses it is fitted to must be reported with the number.
• Relative error over trajectory segments of fixed lengths, which measures drift independent of where the error started.
• Consistency, by NEES averaged over Monte Carlo runs in simulation (Forster et al., 2017) or by checking that errors stay within the reported $3\sigma$ bounds (Hesch et al., 2014).
• Resources: CPU load, memory and time per frame on the target computer, as in Delmerico and Scaramuzza (2018).

Applications

• Flying robots, where cameras and IMUs fit the payload and power budget; the VINS-Mono README says it is primarily designed for drone state estimation and control (vision for drones and UAVs). Safe operation is outside the scope of this article.
• Phones and headsets for augmented reality: VINS-Mono has an iOS port, and Forster et al. (2017) compared their system with a Google Tango device on walking trajectories.
• Ground robots and vehicles, as in MSCKF’s car experiment; robot perception describes how an IMU combined with a camera serves a robot’s state estimation.
• Front ends of visual-inertial SLAM, which add loop closure and maps.

Open Problems

Drift in position and yaw remains, by the observability argument, unless loop closures or absolute sensors are added. Initialization under weak excitation, poorly known biases and fast starts is still fragile, as Martinelli’s analysis and the example show. Consistency of linearized estimators, time synchronization and calibration of low-cost sensors, rolling shutters, and the trade-off between accuracy and onboard computation all limit what the benchmark numbers mean for a given device. Benchmark comparisons themselves depend on alignment and tuning choices that differ between papers.