Direct Sparse Odometry

Direct Sparse Odometry (DSO) is a monocular visual odometry algorithm by Jakob Engel, Vladlen Koltun and Daniel Cremers. It tracks a camera and reconstructs a sparse set of scene points by minimizing a photometric error, the difference between the intensities of the same scene point in different frames, rather than the reprojection error of matched keypoints. It combines two choices that recent systems had kept apart; its authors knew of no recent work combining them. Like other direct methods, it works on pixel intensities and can use any pixel with enough gradient, including edges and smooth shading. Like keypoint-based systems, it is sparse and has no smoothness prior on depth, so it can optimize all its variables jointly: the camera poses, the inverse depth of every point, an affine brightness transfer for every frame and the camera intrinsics, over a sliding window of keyframes. The authors describe this as the photometric equivalent of windowed sparse bundle adjustment (Engel, Koltun and Cremers, 2016).

DSO also models how the camera turns light into pixel values. Each frame is corrected with a calibrated response function and vignetting map, and known exposure times enter the error. A companion paper introduced the photometrically calibrated TUM monoVO dataset and the calibration method (Engel, Usenko and Cremers, 2016). This article follows the arXiv preprint of DSO (v2, October 2016), which is the PDF the authors’ group links as the paper; the journal version in IEEE TPAMI (2018) was not read, and details here are attributed to the preprint and to the official code. The Direct Methods article places DSO among direct and indirect methods and summarizes its comparison with ORB-SLAM.

Problem

Monocular visual odometry estimates the trajectory of a single moving camera from its images, incrementally and without closing loops. A direct formulation inherits two problems that keypoint methods largely avoid: the error is highly non-convex, its linearization is valid only within one or two pixels of the solution, and it is sensitive to everything that changes a pixel’s value other than geometry: exposure changes, the camera’s nonlinear response and lens vignetting (Engel et al., 2016).

Earlier dense and semi-dense direct systems regularized their depth maps with a smoothness prior, which correlates the geometry parameters and makes consistent joint optimization in real time infeasible; Engel et al. (2016) also report having found that priors can introduce a bias and reduce long-term accuracy. DSO drops the prior and the density instead: it samples a limited number of points and optimizes them jointly with everything else.

Inputs and Outputs

Input: a monocular video, the intrinsic calibration of a pinhole camera (DSO undistorts the images first), and, when available, a photometric calibration: the inverse response function, the vignetting map and the exposure time of every frame.
Output: a camera pose for every frame, keyframe poses refined in the window, a sparse point cloud of inverse depths hosted in keyframes, affine brightness parameters per frame and refined intrinsics. The scale of a monocular trajectory is arbitrary and drifts.

The official code reads the inverse response as a file of 256 values, one irradiance value per 8-bit intensity, and the vignetting as a 16-bit or 8-bit image of attenuation factors. It writes the keyframe poses in the TUM RGB-D and TUM monoVO format, as camera-to-world transformations (DSO README, commit 7b0c99f).

Intuition

A DSO point has one unknown, its inverse depth in the frame that hosts it (see Direct Methods), so pixels a detector would never fire on, such as points along an edge, can be used as long as the error around them constrains the depth.

Because points are few and have no prior linking them, the problem has the same sparse structure as bundle adjustment: depths are independent given the poses, the depth block of the Hessian is diagonal, and the depths can be eliminated with the Schur complement. The window is kept small by marginalizing old keyframes and points, which turns their information into a prior on the remaining variables instead of discarding it (see factor graph).

The photometric calibration addresses the other weakness: a direct error rests on brightness constancy, which exposure, response and vignetting break.

Algorithm

The steps below follow Sections 2 and 3 of the preprint, with defaults from the paper.

Photometric correction. Each new frame $I$ is mapped to $I' = G^{-1}(I)/V$ , an image proportional to exposure times irradiance (Section 2.1.2).
Initial frame tracking. All active points are projected into the newest keyframe and slightly dilated, giving a semi-dense depth map. The new frame is aligned against that keyframe alone by two-frame direct image alignment on a multi-scale image pyramid, starting from a constant-motion prediction. If the final RMSE is more than twice that of the previous frame, tracking is assumed to have failed and is retried from up to 27 small rotations on the coarsest level, each try taking about 0.5 ms (Section 3.1, Step 1).
Keyframe decision. A new keyframe is created when a weighted sum of three quantities, measured during initial tracking, exceeds a threshold: the root-mean-square optical flow since the last keyframe, the same flow with rotation removed (which measures translation and hence occlusion), and the relative brightness change between the frames (Section 3.1, Step 2).
Candidate tracking. Non-keyframes refine the depths of candidate points (Section 3.2, Step 2).
New keyframe. Candidates are selected in it, candidates from the window are activated, and the window’s total photometric error is minimized with up to six Gauss–Newton iterations (Sections 2.3 and 3.2).
Marginalization. Keyframes and their points are marginalized so that at most seven keyframes stay active (Section 3.1, Step 3).

The first frame pair needs depths from somewhere. The preprint does not describe its initializer; the code’s README calls it slow and unreliable, and advises a slow initial motion with much translation and little rotation (DSO README).

Mathematical Formulation

Conventions

DSO represents the pose of keyframe $i$ as $T_i \in \mathrm{SE}(3)$ , transforming points from the world frame into the camera frame. Increments $\mathbf{x}_i \in \mathfrak{se}(3)$ are applied on the left, $\mathbf{x}_i \boxplus T_i := e^{\hat{\mathbf{x}}_i} T_i$ (Engel et al., 2016, Eq. 1). The relative pose from host frame $i$ to target frame $j$ is $T_j T_i^{-1}$ . A point $\mathbf{p}$ hosted in frame $i$ with inverse depth $d_\mathbf{p}$ projects into frame $j$ at

\mathbf{p}' = \Pi_c\big(R\,\Pi_c^{-1}(\mathbf{p}, d_\mathbf{p}) + \mathbf{t}\big), \qquad \begin{bmatrix} R & \mathbf{t} \\ \mathbf{0} & 1 \end{bmatrix} := T_j T_i^{-1},

where $\Pi_c$ is the pinhole projection with intrinsics $c$ . The official code stores camera-to-world poses for output (camToWorld) and world-to-camera poses in the optimization, and Stereo DSO’s paper defines $T_i$ the other way round, as camera-to-world (Wang, Schwörer and Cremers, 2017), so formulas from different sources need care.

Photometric camera model and its calibration

DSO uses the image formation model of the TUM monoVO paper:

I_i(\mathbf{x}) = G\big(t_i\, V(\mathbf{x})\, B_i(\mathbf{x})\big),

with irradiance $B_i$ , exposure time $t_i$ , a vignetting map $V: \Omega \to [0, 1]$ and a response function $G: \mathbb{R} \to [0, 255]$ (Engel et al., 2016, Eq. 2). Both $G$ and $V$ are observable only up to a scale factor (Engel, Usenko and Cremers, 2016, Section 2.3). The exposure and image sensor articles describe where these effects come from, and the digital camera article describes the tone curves that make $G$ nonlinear.

Response calibration. For a static scene the vignetting can be absorbed into the irradiance, $B' = V B$ . With $U := G^{-1}$ , images $I_i$ at known exposure times $t_i$ , and Gaussian noise assumed on $U(I_i(\mathbf{x}))$ , Engel, Usenko and Cremers (2016, Eqs. 6–8) minimize

E(U, B') = \sum_i \sum_{\mathbf{x} \in \Omega} \big( U(I_i(\mathbf{x})) - t_i B'(\mathbf{x}) \big)^2

by alternation, since fixing either unknown decouples all the others:

U(k) = \frac{\sum_{(i,\mathbf{x}) \in \Omega_k} t_i B'(\mathbf{x})}{|\Omega_k|}, \qquad B'(\mathbf{x}) = \frac{\sum_i t_i\, U(I_i(\mathbf{x}))}{\sum_i t_i^2},

where $\Omega_k$ is the set of all pixels in all images with value $k$ . Overexposed pixels are removed; the calibration code also marks their neighbours as overexposed by dilating twice (leakPadding = 2) and alternates ten times by default (mono_dataset_code, main_responseCalib.cpp). $U(255)$ is never observed and is extrapolated, and $U$ is scaled so that $U(255) = 255$ . Like the method of Debevec and Malik, this recovers a 256-entry lookup table for the inverse response; unlike it, there is no smoothness prior. Instead the authors use about 1,000 images covering 120 exposure times from 0.05 ms to 20 ms in steps of a factor of 1.05, recorded as a video of a static scene while the exposure changes slowly. They note that a regularized approach would be needed with few images, and that a non-monotonic $U$ would have to be smoothed, which did not happen for their calibration datasets. The high dynamic range imaging article covers response recovery in general.

Vignette calibration. The vignetting is estimated as a dense, non-parametric map from several hundred images of a planar, bright, Lambertian surface, such as a white wall, whose pose is found with an AR marker. With $C$ the unknown irradiance of the plane, discretized into $1000 \times 1000$ points, the energy

E(C, V) = \sum_{i,\, \mathbf{x} \in P} \Big( t_i\, V\big([\pi_i(\mathbf{x})]\big)\, C(\mathbf{x}) - U\big(I_i(\pi_i(\mathbf{x}))\big) \Big)^2

is again minimized by alternating closed-form updates, and $V$ is scaled so that its maximum is 1 (Engel, Usenko and Cremers, 2016, Eqs. 9–11). Unlike earlier dense methods, the surface need not be uniformly coloured. The authors impose no smoothness or radial symmetry; they report that the symmetric, smooth maps they obtained come from the data alone. They argue from random bipartite graphs that 30 images showing the whole plane make the problem almost certainly well defined, and use a short video of about 60 seconds for a good solution. Their calibration code warns that it needs about 16 GB of memory for 1,000 images (mono_dataset_code README).

Correction. DSO applies the model to every frame as its very first step,

I'_i(\mathbf{x}) := t_i B_i(\mathbf{x}) = \frac{G^{-1}\big(I_i(\mathbf{x})\big)}{V(\mathbf{x})},

and works with $I'_i$ from then on (Engel et al., 2016, Eq. 3). The code checks that the inverse response is strictly increasing, rescales it to 0–255, looks up each 8-bit value and multiplies by $1/V$ (DSO, Undistort.cpp). Point selection is the exception: it uses the raw images, before correction (Engel et al., 2016, Section 3.2; setting_gammaWeightsPixelSelect = 1 in the code).

The photometric error

The error of a point $\mathbf{p}$ hosted in frame $i$ and observed in frame $j$ is a weighted sum over a pattern $\mathcal{N}_\mathbf{p}$ of eight pixels around it:

E_{\mathbf{p}j} := \sum_{\mathbf{q} \in \mathcal{N}_\mathbf{p}} w_\mathbf{q} \left\lVert \big(I_j[\mathbf{q}'] - b_j\big) - \frac{t_j e^{a_j}}{t_i e^{a_i}} \big(I_i[\mathbf{q}] - b_i\big) \right\rVert_\gamma ,

where $\mathbf{q}'$ is the projection of pattern pixel $\mathbf{q}$ with the point’s inverse depth (the paper reuses $\mathbf{p}$ for the summation variable), $\lVert\cdot\rVert_\gamma$ is the Huber norm and $a_i, b_i$ are the frame’s affine brightness parameters, which describe the transfer $e^{-a_i}(I_i - b_i)$ (Engel et al., 2016, Eq. 4). The Direct Methods article explains why the gain is parametrized logarithmically. The pattern spreads eight pixels within two pixels of the centre and omits the bottom-right pixel of a nine-pixel pattern so that the computation fits SSE instructions; in the code it is staticPattern[8], the offsets $(0,-2)$ , $(\pm1,-1)$ , $(\pm2,0)$ , $(0,0)$ , $(-1,1)$ and $(0,2)$ . The authors chose eight pixels as a trade-off between computation, robustness to motion blur and information, and note that the SSD over such a neighbourhood carries information similar to adding first- and second-order irradiance derivative constancy terms to the central pixel. A pattern of more than one pixel is necessary: with one unknown per point and no regularizer, a single pixel would leave the parameters under-constrained when only two frames are optimized (Figure 4 caption).

Each pixel is weighted by

w_\mathbf{q} := \frac{c^2}{c^2 + \lVert \nabla I_i(\mathbf{q}) \rVert_2^2},

which down-weights high gradients. The authors interpret it as adding small, independent geometric noise to the projected position and marginalizing it (Eq. 7). In the code, $c = 50$ on the 0–255 intensity scale and the Huber threshold is 9 (setting_outlierTHSumComponent = 50*50, setting_huberTH = 9). The code multiplies each residual by the square root of the ratio, so the squared error carries the ratio itself, and it averages a factor computed from the target frame’s gradient with one stored for the host pixel when the point was created (Residuals.cpp, ImmaturePoint.cpp).

The total error sums over all frames $i$ , all points $\mathbf{p} \in \mathcal{P}_i$ hosted in them and all frames $j \in \mathrm{obs}(\mathbf{p})$ in which they are visible (Eq. 8). Each term depends on the point’s inverse depth, on the intrinsics, on two poses and on four brightness parameters. If exposure times are known, a prior pulls the brightness parameters to zero:

E_\text{prior} := \sum_{i} \big( \lambda_a a_i^2 + \lambda_b b_i^2 \big).

Without a photometric calibration, $t_i = 1$ and $\lambda_a = \lambda_b = 0$ , so the affine parameters must absorb the unknown exposure changes (Eq. 9). The code’s three modes implement these cases: mode=0 uses the calibration with strong priors (setting_affineOptModeA = 1e12, setting_affineOptModeB = 1e8); mode=1 disables the calibration and the priors; mode=2, for synthetic, photometrically perfect images, disables the calibration, fixes $a$ and $b$ , and lowers the gradient threshold for point selection from 7 to 3 (main_dso_pangolin.cpp).

Hessian structure

Engel et al. (2016, Figure 5) draw the error as a factor graph in which each term connects a point, its host keyframe and the keyframe that observes it. The only difference from the classical reprojection error is the dependence on the host frame’s pose. This adds off-diagonal entries to the pose–pose block of the Hessian, but it does not change the sparsity pattern after the Schur complement eliminates the points, so the system can be solved as in the indirect case. The Jacobians with respect to the two poses are linearly related by the adjoint of their relative pose, so the factor can be pulled out of the sums that build the Hessian. The Bundle Adjustment article derives the reduced camera system.

Gauss–Newton with first-estimate Jacobians

All variables are collected in $\zeta \in \mathrm{SE}(3)^n \times \mathbb{R}^m$ , written as $\zeta = \mathbf{x} \boxplus \zeta_0$ : an evaluation point $\zeta_0$ that fixes the tangent space, and accumulated increments $\mathbf{x}$ . Following the approach of Leutenegger et al. (2015), as the authors write, the window is optimized with Gauss–Newton, $H = J^\top W J$ , $\mathbf{b} = -J^\top W \mathbf{r}$ , and $\boldsymbol\delta = H^{-1}\mathbf{b}$ is added to $\mathbf{x}$ (Engel et al., 2016, Eqs. 10–14). Each residual’s Jacobian splits into the image gradient $J_I$ , a geometric part $J_\text{geo}$ (poses, inverse depth, intrinsics) and a photometric part $J_\text{photo}$ (the four brightness parameters), and two approximations are made:

First-estimate Jacobians. $J_\text{geo}$ and $J_\text{photo}$ are evaluated at $\mathbf{x} = 0$ , at the evaluation point. The energy has non-linear null spaces, absolute pose and scale, and adding linearizations taken around different points would eliminate them and slowly corrupt the system. The image gradient $J_I$ is much less smooth but does not affect the null spaces, so it is evaluated at the current estimate, from centred differences at integer positions interpolated bilinearly.
One geometric Jacobian per point. $J_\text{geo}$ is evaluated at the centre pixel and shared by the eight residuals of the pattern; the authors report no notable effect on accuracy on their datasets.

After each step, the evaluation point is moved to the current estimate for every variable not part of the marginalization prior, in practice all depths and the newest keyframe’s pose. When a keyframe is added, up to six iterations run, ending early when the step is small. The authors found Levenberg–Marquardt damping unnecessary, since the optimization never starts far from the minimum (Section 2.3). The released code has a damping parameter, but its default solver mode fixes it at $10^{-5}$ and accepts every step (setting_solverMode = SOLVER_FIX_LAMBDA | ..., setting_forceAceptStep = true, EnergyFunctional.cpp).

Marginalization

When the window grows too large, old variables are marginalized with the Schur complement. The part of the energy $E'$ that contains the residuals depending on the variables to be removed is approximated to second order around the current state,

E'(\mathbf{x} \boxplus \zeta_0) \approx 2\,\mathbf{x}^\top \underbrace{(\mathbf{b} - H \mathbf{x}_0)}_{\mathbf{b}'} + \mathbf{x}^\top H \mathbf{x} + \text{const},

and the variables $\beta$ to be removed are eliminated from the variables $\alpha$ to be kept by the Schur complement, giving $\hat{H}_{\alpha\alpha}$ and $\hat{\mathbf{b}}'_\alpha$ (the factor graph article covers marginalization in information form). The quadratic $2\mathbf{x}_\alpha^\top \hat{\mathbf{b}}'_\alpha + \mathbf{x}_\alpha^\top \hat{H}_{\alpha\alpha} \mathbf{x}_\alpha$ replaces the removed nonlinear terms in all later optimization, which is why the tangent space of every variable in it must stay fixed (Engel et al., 2016, Eqs. 15–19). To keep the Hessian sparse, residuals that would fill it in are dropped: when frame $i$ is marginalized, first all points it hosts and all points not observed in the last two keyframes are marginalized, and the remaining observations of active points in frame $i$ are dropped. The authors call this clearly suboptimal; in practice about half of all residuals are dropped for this reason (Section 3.1, Step 3). In the code, the Hessian and gradient are also multiplied by $0.5^2$ when added to the prior, a factor commented as accounting for inaccurate linearization points (setting_margWeightFac, EnergyFunctional.cpp); the paper does not mention it.

Which keyframes to marginalize

With $I_1$ the newest and $I_n$ the oldest active keyframe, the rules are (Section 3.1, Step 3):

The two newest keyframes, $I_1$ and $I_2$ , are always kept.
Keyframes with less than 5% of their points visible in $I_1$ are marginalized.
If more than $N_f = 7$ keyframes remain active, the one (other than $I_1$ and $I_2$ ) with the highest distance score is marginalized:

s(I_i) = \sqrt{d(i, 1)} \sum_{j \in [3, n] \setminus \{i\}} \big(d(i, j) + \epsilon\big)^{-1},

with $d(i, j)$ the Euclidean distance between keyframes and $\epsilon$ a small constant. The heuristic keeps keyframes spread in space, with more of them near the newest one. In the code, a keyframe is also marginalized when its brightness differs from the newest frame’s by a factor larger than $e^{0.7}$ , and at least five keyframes are kept (setting_maxLogAffFacInWindow = 0.7, setting_minFrames = 5, FullSystemMarginalize.cpp).

Point Selection and Management

DSO aims to keep a fixed number $N_p = 2000$ of active points, spread evenly across the image and the window (Engel et al., 2016, Section 3.2). The authors report that image data are highly redundant and that the benefit of using more points quickly flattens off, so they sub-sample heavily.

Candidate selection. Points should be well distributed and have a high gradient relative to their surroundings. The image is split into $32 \times 32$ blocks, each with a region-adaptive threshold $\bar g + g_\text{th}$ : the median absolute gradient $\bar g$ of the block plus a global constant $g_\text{th} = 7$ . The image is then split into $d \times d$ cells, and from each cell the pixel with the largest gradient is selected if it passes its block’s threshold. To also capture weak intensity variations, such as smooth illumination changes on white walls, the procedure is repeated twice with a decreased threshold and cells of $2d$ and $4d$ . The cell size $d$ adapts from frame to frame so that the procedure yields the desired number of candidates. The paper’s Figure 9 shows 2,000 candidates per frame, the first pass spread over gradient-rich areas and the later passes, much sparser, also covering weakly textured regions.

The released code differs in details (PixelSelector2.cpp, settings.cpp). The block threshold is the median plus 7 of a gradient histogram, averaged over the neighbouring blocks. Within a cell, it does not take the largest gradient magnitude but the largest projection of the gradient onto a random direction drawn for that cell, which, in our reading, spreads the gradient directions of the selected points (setting_selectDirectionDistribution = true). The second and third passes test gradients on the next two pyramid levels against thresholds on squared gradients that are lowered by a factor of 0.75 and of $0.75^3$ (setting_gradDownweightPerLevel). A cell of the larger size gets a point only when none of its smaller cells did. The default asks for 1,500 candidates per keyframe (setting_desiredImmatureDensity) and 2,000 active points; when more candidates than needed are found, the excess is removed with a fixed pseudo-random pattern seeded for determinism.

Candidate tracking. Candidates are not optimized at once. They are tracked in subsequent frames by a discrete search along the epipolar line that minimizes the photometric error, an approach the authors say is inspired by LSD-SLAM. The best match gives a depth and a variance, which bound the next search. Candidates whose minimum is not distinct enough are discarded permanently, which removes many false matches in repetitive regions (Section 3.2).

Activation. When points have been marginalized, candidates from all frames in the window are activated to replace them. All active points are projected into the newest keyframe, and the candidates that, also projected there, lie farthest from any existing point are activated, with a larger required distance for candidates from the second and third passes.

Outliers. Observations whose photometric error exceeds a threshold are removed. The threshold adapts to the median residual of the frame, so it is looser for bad frames with much motion blur and stricter for good ones.

Parameters

The defaults in the paper and in the released code (commit 7b0c99f) agree on most points:

Parameter	Paper (preprint)	Code default
Active keyframes $N_f$	7	5 to 7 (`setting_minFrames`, `setting_maxFrames`)
Active points $N_p$	2,000	2,000; 1,500 candidates per keyframe
Gradient threshold $g_\text{th}$	7	7 (3 in `mode=2`)
Residual pattern	8 pixels	`staticPattern[8]`
Huber threshold	not given	9 intensity units
Gradient weight constant $c$	not given	50
Gauss–Newton iterations per keyframe	up to 6	1 to 6
Keyframe threshold $T_\text{kf}$	1, about 7.5 keyframes per second on TUM monoVO (Figure 18; the text says 8)	weighted flow and brightness sum above 1
Re-tracking	RMSE more than twice the previous frame’s; up to 27 small rotations	RMSE more than 1.5 times the previous frame’s; five motion guesses and 26 rotations

The re-tracking row compares the paper’s description with trackNewCoarse in FullSystem.cpp. There, the candidates are a constant-motion prediction, double, half and zero motion, zero motion from the keyframe, and 26 small rotations around the constant-motion prediction, and the search stops at the first that brings the RMSE below 1.5 times the previous frame’s (setting_reTrackThreshold = 1.5). The keyframe test in the code also creates a keyframe when the tracking RMSE has doubled since the first frame tracked against the current keyframe.

The parameter studies on TUM monoVO (Engel et al., 2016, Section 4.2) indicate how sensitive these choices are. Using more than about 500 points or 7 active keyframes had only marginal effect on the alignment error, and runtime grows quadratically with $N_f$ ; the defaults of 2,000 points and 7 keyframes were chosen mainly to give denser reconstructions. Always marginalizing the oldest keyframe (fixed lag) instead of using the distance score performed significantly worse. The gradient threshold had relatively low impact, with a sweet spot around 7. Taking fewer than 4 keyframes per second reduced robustness, mainly through strong occlusions such as walking through doors, and taking more than 15 reduced accuracy, because keyframes were marginalized earlier around less accurate linearization points. A $3 \times 3$ pattern seemed slightly worse, and patterns larger than eight pixels brought little; the authors add that the result may depend on the camera’s point spread function.

Complexity

The cost per keyframe is dominated by the window optimization over at most $N_f$ keyframes and $N_p$ points, each point contributing eight residuals per observing keyframe. The authors report that their CPU implementation runs in real time on a laptop, and that with reduced settings (800 points, 6 active keyframes, $424 \times 320$ pixels, at most 4 Gauss–Newton iterations per keyframe) it ran at five times real-time speed while still giving very good accuracy and robustness on all three of their datasets (Engel et al., 2016, Sections 1.2 and 4.1). Code presets 2 and 3 implement these fast settings (main_dso_pangolin.cpp).

For the evaluation, DSO ran in a sequentialized, single-threaded mode about four times slower than real time, with the same parameters as in real-time runs. In the real-time, multi-threaded mode, tracking new frames and creating keyframes run in parallel, so a new frame is sometimes tracked against the second-newest keyframe; the authors note that this occasionally makes initial alignment fail, particularly during strong exposure changes (Section 4.1).

Implementation

The official implementation is written in C++ and requires Eigen and SuiteSparse. OpenCV, for reading and displaying images, and Pangolin, for the GUI, are optional but recommended; without OpenCV the dataset binary is not built. It is licensed under GPLv3, and the README notes that a commercial version exists (DSO README, as of November 2026). The calibration code for the TUM monoVO dataset, including the response and vignette calibration and the MATLAB evaluation of the alignment error, is a separate repository under a BSD licence (mono_dataset_code). The README gives practical advice that follows from the paper’s findings: calibrate the geometry carefully, because geometric distortions of 1.5 pixels already reduce accuracy by a factor of ten; do not use a rolling-shutter camera; use a photometric calibration; and avoid rotating without translating.

Example: photometric calibration in sparse alignment

The example tracks one frame against a reference with DSO’s ingredients and asks which parts of the photometric calibration matter when the images carry vignetting or a nonlinear response. It is not DSO: there is no window, no marginalization and no depth estimation. The Direct Methods article’s example covers global exposure changes and clipping in dense alignment; this one uses neither.

A three-plane scene (a back wall 3 m away, a floor and a side wall, textured with images bundled with scikit-image) is rendered as irradiance $B \in [0.05, 0.45]$ into a $320 \times 240$ camera with a focal length of 250 pixels. Each image is recorded as $I = G(t\,V B)$ plus Gaussian noise of one intensity level, rounded to 8 bits. The vignetting is the natural $\cos^4$ fall-off of the off-axis angle, which reaches 0.37 in the corners. The responses are our choices, not measured curves: linear, a gamma curve $E^{1/2.2}$ , and a saturating curve $(1 - e^{-3E})/(1 - e^{-3})$ . A pure power law turns an exposure change into an exact gain, $G(2E) = 2^{1/2.2} G(E)$ , which an affine model can absorb, so the shape of the curve is an open choice and both nonlinear curves are printed. The second frame is taken after a motion of 4.40° and 226 mm, with the same exposure time or twice the reference’s, which saturates no pixel.

Points are selected on the raw reference image, as DSO does, by a simplified version of its rule: block thresholds of the median gradient plus 7 in $32 \times 32$ blocks, and three passes over $8 \times 8$ , $16 \times 16$ and $32 \times 32$ cells, with the thresholds multiplied by 1, 0.75 and $0.75^2$ . The simplifications against the code are fixed cell sizes, no smoothing of thresholds between blocks, no random gradient direction, full-resolution gradients in every pass and thresholds on gradient magnitudes rather than squared gradients. Each point uses DSO’s eight-pixel pattern, the gradient weight with $c = 50$ , a Huber threshold of 9 and one geometric Jacobian per point. The points’ true depths are given, so the example isolates the photometric model. The residual is $r = I_t[\mathbf{p}'] - b - k\,e^{a} I_r[\mathbf{p}]$ , with the pose convention $\mathbf{X}_t = R\,\mathbf{X}_r + \mathbf{t}$ and increments applied on the left, as in DSO. Gauss–Newton runs on a three-level pyramid built by $2 \times 2$ averaging. The five models follow the ablation of Engel et al. (2016, Figure 15):

full calibration: images corrected to $G^{-1}(I)/V$ , the exposure ratio $k$ known, $a = b = 0$ fixed;
unknown t: corrected images, $k = 1$ , $a$ and $b$ estimated;
unknown t, V: only $G^{-1}$ applied, $a$ and $b$ estimated;
unknown t, V, G: raw images, $a$ and $b$ estimated;
brightness constancy: raw images, $a = b = 0$ .

Every model starts from the same motion prediction, half the true rotation and translation. This stands in for DSO’s constant-motion initialization (Section 3.1) and is our choice. The plan was to start from the identity. A pilot on images without photometric distortion found all five models failing from there, and two diagnostic runs on the same images, one comparing residuals at the identity and at the true pose and one starting the alignment from fractions of the true motion, showed that the aligner itself worked. The half-motion start was then fixed by editorial decision and confirmed by a second pilot on the same undistorted images, before any distorted condition was run. After the first full run, only the table header was changed, because it printed the condition names truncated and run together; no number changed.

import sys
import numpy as np
from scipy import ndimage
from scipy.spatial.transform import Rotation
from skimage import color, data

H, W, f = 240, 320, 250.0
cx, cy = (W - 1) / 2, (H - 1) / 2

gray = lambda im: color.rgb2gray(im) if im.ndim == 3 else im / 255.0
# Planes n . X = d in the reference camera frame (x right, y down, z forward), with a texture
# and the two coordinates used as texture axes (column, row), tiled at 120 texels per metre.
planes = [([0, 0, 1.0], 3.0, gray(data.chelsea()), 0, 1),    # back wall
          ([0, 1.0, 0], 0.8, gray(data.brick()), 0, 2),      # floor
          ([1.0, 0, 0], -1.2, gray(data.coffee()), 2, 1)]    # left wall

def render(R, t, ss=3):
    """Irradiance and depth seen by a camera with X_cam = R X_ref + t (ss x ss supersampling)."""
    v, u = np.mgrid[0:H * ss, 0:W * ss] / ss + (1 / ss - 1) / 2
    rays = np.stack([(u - cx) / f, (v - cy) / f, np.ones_like(u)], -1) @ R
    c = -R.T @ t
    depth, tex = np.full(u.shape, np.inf), np.zeros(u.shape)
    for n, d, im, a0, a1 in planes:
        with np.errstate(divide="ignore", invalid="ignore"):
            s = (d - c @ n) / (rays @ n)
        s[~np.isfinite(s) | (s <= 0)] = np.inf
        X = c + s[..., None] * rays
        val = ndimage.map_coordinates(im, [X[..., a1] * 120, X[..., a0] * 120], order=1, mode="grid-wrap")
        hit = s < depth
        tex[hit], depth[hit] = val[hit], s[hit]
    B = 0.05 + 0.4 * tex.reshape(H, ss, W, ss).mean((1, 3))
    return B, depth.reshape(H, ss, W, ss).mean((1, 3))

R_true = Rotation.from_euler("xyz", [1.5, 4.0, -1.0], degrees=True).as_matrix()
t_true = np.array([0.20, 0.03, 0.10])
# Motion prediction (as from DSO's constant-motion model): half the true motion, the same for every model.
R_pred = Rotation.from_rotvec(0.5 * Rotation.from_matrix(R_true).as_rotvec()).as_matrix()
t_pred = 0.5 * t_true

v_, u_ = np.mgrid[0:H, 0:W]
V_cos4 = (1 + ((u_ - cx) ** 2 + (v_ - cy) ** 2) / f ** 2) ** -2          # cos^4 of the off-axis angle
responses = {"linear": lambda E: 255 * E,
             "gamma": lambda E: 255 * E ** (1 / 2.2),
             "saturating": lambda E: 255 * (1 - np.exp(-3 * E)) / (1 - np.exp(-3))}
LUT_E = np.linspace(0, 1, 100001)

def inverse_response(name):
    """256-entry inverse response, mapping each 8-bit value to irradiance x 255 (as DSO's pcalib file)."""
    Gv = responses[name](LUT_E)
    return 255 * np.interp(np.arange(256), Gv, LUT_E)

def capture(B, V, resp, t_exp, rng):
    I = responses[resp](np.clip(t_exp * V * B, 0, 1)) + rng.standard_normal(B.shape)
    return np.clip(np.round(I), 0, 255).astype(np.uint8)

PATTERN = np.array([[0, -2], [-1, -1], [1, -1], [-2, 0], [0, 0], [2, 0], [-1, 1], [0, 2]])  # DSO staticPattern[8]

def select_points(I, d=8, g_th=7.0):
    """DSO-style candidate selection on the raw image (simplified, see text)."""
    I = I.astype(float)
    gy, gx = np.gradient(I)
    g = np.hypot(gx, gy)
    th = np.zeros_like(g)
    for y in range(0, H, 32):
        for x in range(0, W, 32):
            th[y:y + 32, x:x + 32] = np.median(g[y:y + 32, x:x + 32]) + g_th
    margin = 4
    valid = np.zeros_like(g, bool); valid[margin:H - margin, margin:W - margin] = True
    chosen = np.zeros_like(g, bool)
    counts = []
    for cell, fac in [(d, 1.0), (2 * d, 0.75), (4 * d, 0.75 ** 2)]:
        n = 0
        for y in range(0, H, cell):
            for x in range(0, W, cell):
                if chosen[y:y + cell, x:x + cell].any():
                    continue
                gg = np.where(valid[y:y + cell, x:x + cell] & (g[y:y + cell, x:x + cell] > fac * th[y:y + cell, x:x + cell]),
                              g[y:y + cell, x:x + cell], -1)
                k = np.argmax(gg)
                if gg.flat[k] > 0:
                    chosen[y + k // gg.shape[1], x + k % gg.shape[1]] = True
                    n += 1
        counts.append(n)
    vv, uu = np.nonzero(chosen)
    return np.column_stack([uu, vv]).astype(float), counts

def pyramid(I, levels=3):
    out = [I]
    for _ in range(levels - 1):
        I = out[-1]
        h, w = I.shape[0] // 2 * 2, I.shape[1] // 2 * 2
        out.append(I[:h, :w].reshape(h // 2, 2, w // 2, 2).mean((1, 3)))
    return out

def twist_rows(Ju, Jv, X):
    x, y, z = X.T
    Jz = -(Ju * x + Jv * y) / z
    return np.stack([Ju, Jv, Jz, -Jv * z + Jz * y, Ju * z - Jz * x, -Ju * y + Jv * x], 1)

def align(Ir, It, pts, Z, k, affine, R0, t0, levels=3, iters=30, huber=9.0, c=50.0):
    """Sparse direct alignment: r = I_t[p'] - b - k e^a I_r[p] over DSO's 8-pixel pattern."""
    Pr, Pt = pyramid(Ir, levels), pyramid(It, levels)
    R, t, a, b = R0.copy(), t0.copy(), 0.0, 0.0
    for l in reversed(range(levels)):
        s = 2.0 ** l
        fl, cxl, cyl = f / s, (cx + 0.5) / s - 0.5, (cy + 0.5) / s - 0.5
        ul, vl = (pts[:, 0] + 0.5) / s - 0.5, (pts[:, 1] + 0.5) / s - 0.5
        X = np.column_stack([(ul - cxl) / fl * Z, (vl - cyl) / fl * Z, Z])
        Ir_l, It_l = Pr[l], Pt[l]
        gy_r, gx_r = np.gradient(Ir_l)
        gy_t, gx_t = np.gradient(It_l)
        qu = ul[:, None] + PATTERN[None, :, 0]; qv = vl[:, None] + PATTERN[None, :, 1]
        ref = ndimage.map_coordinates(Ir_l, [qv.ravel(), qu.ravel()], order=1).reshape(qu.shape)
        gr2 = (ndimage.map_coordinates(gx_r, [qv.ravel(), qu.ravel()], order=1) ** 2
               + ndimage.map_coordinates(gy_r, [qv.ravel(), qu.ravel()], order=1) ** 2).reshape(qu.shape)
        wg = c ** 2 / (c ** 2 + gr2)
        for _ in range(iters):
            Xc = X @ R.T + t
            pu, pv = fl * Xc[:, 0] / Xc[:, 2] + cxl, fl * Xc[:, 1] / Xc[:, 2] + cyl
            Qu, Qv = pu[:, None] + PATTERN[None, :, 0], pv[:, None] + PATTERN[None, :, 1]
            ok = ((Qu > 0) & (Qu < It_l.shape[1] - 1) & (Qv > 0) & (Qv < It_l.shape[0] - 1)).all(1)
            at = [Qv[ok].ravel(), Qu[ok].ravel()]
            Iw, Gx, Gy = (ndimage.map_coordinates(im, at, order=1) for im in (It_l, gx_t, gy_t))
            rf = ref[ok].ravel()
            r = Iw - b - k * np.exp(a) * rf
            Xo = np.repeat(Xc[ok], 8, 0)                      # centre-pixel geometry for the whole pattern
            Ju, Jv = Gx * fl / Xo[:, 2], Gy * fl / Xo[:, 2]
            J = twist_rows(Ju, Jv, Xo)
            if affine:
                J = np.hstack([J, -k * np.exp(a) * rf[:, None], -np.ones((len(r), 1))])
            w = wg[ok].ravel() * np.minimum(1.0, huber / np.maximum(np.abs(r), 1e-12))
            dx = np.linalg.solve((J * w[:, None]).T @ J, -(J * w[:, None]).T @ r)
            dR = Rotation.from_rotvec(dx[3:6]).as_matrix()
            R, t = dR @ R, dR @ t + dx[:3]
            if affine:
                a, b = a + dx[6], b + dx[7]
            if np.linalg.norm(dx) < 1e-6:
                break
    return R, t, np.exp(a), b

def pose_error(R, t):
    return np.degrees(Rotation.from_matrix(R @ R_true.T).magnitude()), 1000 * np.linalg.norm(t - t_true)

B_ref, Z_ref = render(np.eye(3), np.zeros(3))
B_tgt, _ = render(R_true, t_true)
print(f"true motion {np.degrees(Rotation.from_matrix(R_true).magnitude()):.2f} deg, {1000*np.linalg.norm(t_true):.0f} mm;"
      f" depth {Z_ref.min():.2f}-{Z_ref.max():.2f} m; corner vignetting {V_cos4.min():.2f}")

def run_condition(resp, vign, ratio):
    rng = np.random.default_rng(0)
    V = V_cos4 if vign else np.ones((H, W))
    I_r = capture(B_ref, V, resp, 1.0, rng)
    I_t = capture(B_tgt, V, resp, ratio, rng)
    pts, counts = select_points(I_r)
    Z = ndimage.map_coordinates(Z_ref, [pts[:, 1], pts[:, 0]], order=1)
    Ginv = inverse_response(resp)
    raw = lambda I: I.astype(float)
    g_only = lambda I: Ginv[I]
    full = lambda I: Ginv[I] / V
    models = [("full calibration", full, ratio, False),
              ("unknown t", full, 1.0, True),
              ("unknown t, V", g_only, 1.0, True),
              ("unknown t, V, G", raw, 1.0, True),
              ("brightness constancy", raw, 1.0, False)]
    res = []
    for name, corr, k, aff in models:
        R, t, g, o = align(corr(I_r), corr(I_t), pts, Z, k, aff, R_pred, t_pred)
        res.append((name,) + pose_error(R, t) + (g,))
    return counts, res

conditions = [("none", "linear", False, 1.0, "none"), ("vignetting", "linear", True, 1.0, "V"),
              ("gamma, t x2", "gamma", False, 2.0, "gamma"), ("saturating, t x2", "saturating", False, 2.0, "sat"),
              ("gamma+V, t x2", "gamma", True, 2.0, "gamma+V"), ("saturating+V, t x2", "saturating", True, 2.0, "sat+V")]
results = {}
for label, resp, vign, ratio, short in conditions:
    counts, res = run_condition(resp, vign, ratio)
    results[label] = res
    print(f"{label:19s} points per pass {counts}, total {sum(counts)}")
names = [r[0] for r in results["none"]]
for title, col, fmt in [("translation error (mm)", 2, "{:9.2f}"), ("rotation error (deg)", 1, "{:9.3f}"),
                        ("estimated gain e^a", 3, "{:9.3f}")]:
    print(f"\n{title}")
    print(" " * 22 + "".join(f"{c[4]:>9s}" for c in conditions))
    for m, name in enumerate(names):
        print(f"{name:22s}" + "".join(fmt.format(results[c[0]][m][col]) for c in conditions))

Output (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2; the script runs in about 5 seconds):

true motion 4.40 deg, 226 mm; depth 1.67-3.00 m; corner vignetting 0.37
none                points per pass [443, 42, 3], total 488
vignetting          points per pass [325, 53, 2], total 380
gamma, t x2         points per pass [471, 38, 2], total 511
saturating, t x2    points per pass [618, 27, 2], total 647
gamma+V, t x2       points per pass [406, 47, 2], total 455
saturating+V, t x2  points per pass [562, 33, 3], total 598

translation error (mm)
                           none        V    gamma      sat  gamma+V    sat+V
full calibration           0.47     0.72     0.54     0.47     0.66     0.55
unknown t                  0.59     0.77     0.67     0.36     0.28     0.63
unknown t, V               0.59     4.48     0.67     0.36     4.45    10.88
unknown t, V, G            0.59     4.48     0.66     1.60     7.15    16.12
brightness constancy       0.47     7.18   945.55  1184.44  2137.80  2050.21

rotation error (deg)
                           none        V    gamma      sat  gamma+V    sat+V
full calibration          0.006    0.006    0.008    0.004    0.006    0.007
unknown t                 0.010    0.012    0.012    0.005    0.005    0.010
unknown t, V              0.010    0.049    0.012    0.005    0.109    0.239
unknown t, V, G           0.010    0.049    0.011    0.025    0.171    0.354
brightness constancy      0.006    0.078   15.945   16.712   33.480   23.493

estimated gain e^a
                           none        V    gamma      sat  gamma+V    sat+V
full calibration          1.000    1.000    1.000    1.000    1.000    1.000
unknown t                 0.941    0.936    1.893    1.905    1.891    1.897
unknown t, V              0.941    0.904    1.893    1.905    1.805    1.814
unknown t, V, G           0.941    0.904    1.302    1.020    1.209    1.150
brightness constancy      1.000    1.000    1.000    1.000    1.000    1.000

Without photometric distortion, every model recovers the pose to about half a millimetre and 0.01°, so the differences in the other columns come from the photometric model, not from the optimizer. The free affine gain settles below 1 although nothing changed between the frames, with little effect on the pose; in our reading, this is the drift of an unconstrained gain estimated from two noisy images that the Direct Methods article describes.

Vignetting alone is not a global gain. A point moves across the image between the frames, so the ratio of the attenuation factors at its two positions varies from point to point, and a single gain and offset can absorb only part of it. The models that do not divide by $V$ are biased to about 4.5 mm, six times the error with the vignetting removed, and brightness constancy reaches 7.2 mm. With a linear response, applying $G^{-1}$ changes nothing, so the “unknown t, V” and “unknown t, V, G” rows coincide in this column.

With a doubled exposure, brightness constancy fails completely, with errors of about a metre or more and 16° to 33°. Under the gamma response, the affine model on raw images recovers the pose, as the power law implies; its gain differs from $2^{1/2.2} = 1.37$ , in our reading because the offset, which is not printed, is also free. Under the saturating response, the transfer between the two raw images is no longer affine, and the same model’s error rises to 1.6 mm, between three and four and a half times the calibrated models’. With vignetting added to either response, only the models that remove both $G$ and $V$ stay below 0.7 mm; the others are biased by several millimetres or more.

Knowing the exposure time adds little once $G$ and $V$ are removed: the “unknown t” row is as accurate as full calibration, and its estimated gain, close to 1.9, stands in for the true ratio of 2. The paper’s ablation, described under Results, also found known exposure times to have little effect and brightness constancy clearly worst, but the set-ups differ: it ran the full system on 50 real sequences with their measured calibration and scored drift over whole trajectories, while the example aligns one synthetic frame pair with chosen curves. The other findings do not carry over. In the paper, removing the vignette and response calibration decreased accuracy only slightly, and removing only the vignette was slightly worse than removing both; here vignetting causes the largest errors after brightness constancy, and under the nonlinear responses removing both $G$ and $V$ is worse than removing only $V$ . In our reading, the strong $\cos^4$ fall-off and points reaching the image corners make vignetting matter more here. The number of selected points changes with the condition, from 380 to 647, because selection runs on the raw images, whose gradients depend on $G$ and $V$ . The example uses one noise seed and one motion, and gives each model the true depths; in DSO the depths, poses and brightness parameters are estimated together, so, in our reading, an unmodelled photometric effect would also bias the depths.

Properties and Behavior

Evaluation set-up

Engel et al. (2016, Section 4) evaluate on three datasets: TUM monoVO (50 photometrically calibrated sequences, 105 minutes), EuRoC MAV (11 sequences, no photometric calibration, so $\lambda_a = \lambda_b = 0$ and no correction) and the synthetic ICL-NUIM sequences (constant virtual exposure, $t = 1$ ). TUM monoVO is scored by the alignment error $e_\text{align}$ of its own paper, the others by the absolute trajectory error after Sim(3) alignment. On EuRoC the shaky start of each sequence, meant to initialize IMU biases, was cropped. Every sequence was run forwards and backwards, five times in each direction, giving 500 runs on TUM monoVO, 220 on EuRoC (left and right camera separately) and 80 on ICL-NUIM. Results are shown as cumulative error plots (see Direct Methods).

The TUM monoVO benchmark

The dataset was recorded with uEye UI-3241LE-M-GL cameras, monochrome global-shutter CMOS cameras capable of $1280 \times 1024$ pixels at up to 60 frames per second, with two fisheye lenses of $148° \times 122°$ and $98° \times 79°$ field of view (Engel, Usenko and Cremers, 2016, Section 2.1). Sequences were recorded at 20 to 50 frames per second with JPEG compression, some with hardware gamma correction and some without, and come with a camera model, exposure times as reported by the sensor, a response function and a vignetting map. In one sequence of the DSO paper’s Figure 3, with indoor and outdoor parts, the exposure varies by a factor of more than 500, from 0.018 to 10.5 ms. The authors contrast the dataset with those then in use: KITTI is limited to forward motion in street scenes and comes pre-rectified, TUM RGB-D contains strong motion blur, rolling-shutter artifacts and rotation-only motions that monocular odometry cannot track well, and the rendered ICL-NUIM sequences need no photometric calibration because their virtual exposure is constant; in contrast to existing datasets, they write, all monoVO sequences are photometrically calibrated (abstract and Section 1.1).

Instead of full ground truth, each sequence starts and ends with 10 to 20 seconds of slow, loopy motion over the same scene. Those segments are tracked and aligned precisely with LSD-SLAM, or with motion capture where available, and scale is normalized so that the trajectory is about 100 units long. A method’s trajectory is aligned in Sim(3) separately to the start and the end segment, giving $T_s$ and $T_e$ . The accumulated drift $T_e T_s^{-1}$ yields scale, rotation and translation drift, and the alignment error is the RMSE between the trajectory aligned by the two transformations (Section 3, Eq. 14):

e_\text{align} := \sqrt{\frac{1}{n} \sum_{i=1}^{n} \big\lVert T_s \mathbf{p}_i - T_e \mathbf{p}_i \big\rVert_2^2 } .

The authors chose it because it weights scale, rotation and translation drift by their effect on position, and depends less on when the drift occurs than the translation drift or the RMSE over the two segments. The protocol is valid only for methods that do not close loops, so loop closure must be disabled in SLAM systems. The authors argue that drift before a loop closure is still a good indicator of accuracy, while loop detection, relocalization and map maintenance are not evaluated. The dataset paper’s own analyses, on pinhole-rectified images, found that a smaller field of view significantly decreased accuracy and robustness for both DSO and ORB-SLAM, that DSO was only marginally affected by resolution while ORB-SLAM was strongly affected, which they attribute to the sub-pixel accuracy of direct approaches, and that ORB-SLAM performed significantly better on sequences played backwards while DSO was largely unaffected, which the authors call an example of dataset bias (Section 4, Figures 9–11).

Results

The Direct Methods article summarizes the comparison with ORB-SLAM on the three datasets and the geometric-noise study (Section 4.3).

The photometric calibration ablation (Section 4.2, Figure 15) disables the components one at a time on TUM monoVO: first the exposure times ( $t_i = 1$ , $\lambda_a = \lambda_b = 0$ ), then also the vignetting ( $V = 1$ ), then also the response ( $G^{-1}$ the identity), and finally the affine brightness correction (brightness constancy). Known exposure times seemed to have little effect on accuracy; removing the vignette and response calibration slightly decreased accuracy and robustness; removing only the vignette calibration performed slightly worse than removing both; and brightness constancy performed clearly worst, since it does not account for automatic exposure changes at all.

Restricting candidates to FAST corners significantly decreased performance, which the authors take as evidence that using edges and weakly textured regions brings a real benefit. With low FAST thresholds, many false corners appear along edges; DSO can still use them, while for an indirect method they would be outliers (Figure 17).

Behaviour

Density on demand. The reconstruction density follows $N_p$ : in textured scenes the points become sparse, while in scenes with few informative regions they approach the density of LSD-SLAM’s semi-dense maps; the authors report point clouds similar to LSD-SLAM’s but much more accurate in high, non-real-time settings (Engel et al., 2016, Sections 1.2 and 4.4).
Consistency of the observations. Overlapping patterns can use some pixel values more than once and others not at all, particularly in scenes with little texture. The authors argue the effect is negligible and could be avoided by removing or down-weighting observations that share pixels (Section 2.2).

Limitations

Pure odometry. DSO marginalizes everything that leaves the field of view, so it neither closes loops nor relocalizes, and drift accumulates. On EuRoC, when ORB-SLAM was prevented from matching points not seen in the last 10 seconds, removing the small loops it closes implicitly, its accuracy became similar to DSO’s, though it was less robust (Engel et al., 2016, Section 4.1).
Monocular scale. Scale is unobservable and drifts. Wang, Schwörer and Cremers (2017) show monocular DSO trajectories on KITTI with substantial scale drift that a single best-fit scale cannot remove.
Geometric calibration and rolling shutter. The direct model does not model geometric noise, and the authors’ noise study (see Direct Methods) suggests it is more susceptible to inaccurate intrinsics than indirect methods. They conclude that indirect methods suit smartphones and webcams, whose cameras were designed for human viewing, and direct methods suit dedicated machine-vision cameras (Section 4.3).
Photometric calibration. Without it, DSO falls back on the affine model. Wang et al. (2017) note that the performance of direct methods like DSO substantially degrades without photometric calibration, which many datasets do not provide, and their stereo system still failed on a Cityscapes scene combining a severe auto-exposure change with strong rotation.
Pinhole model only. Radial distortion is removed beforehand, which for wide-angle lenses reduces the field of view; the authors note the approach can be extended to other invertible camera models at higher cost (Section 2.1.1).
Non-convexity. The energy has the same structure as an indirect one, so, the authors write, it can be integrated with frameworks such as double-window bundle adjustment or incremental smoothing and mapping. The main challenge they name is the much greater non-convexity that comes from having the image in the error, which they expect to restrict the model to video (Section 5).

Variants

Stereo DSO

Wang, Schwörer and Cremers (2017) extend DSO to a stereo camera. Each point now has temporal multi-view residuals, as in DSO, and static stereo residuals against the right image of its host keyframe,

r^s_k = I^R_i[\mathbf{p}'(T_{ji}, d, c)] - b^R_i - \frac{e^{a^R_i}}{e^{a^L_i}} \big(I_i[\mathbf{p}] - b^L_i\big),

in which the left-to-right transformation $T_{ji}$ is fixed, so only the depth and intrinsics are geometric parameters. Each keyframe has separate affine parameters for its left and right images. A coupling factor $\lambda$ weights the static stereo energy against the temporal one; on KITTI sequence 06, $\lambda = 1$ or 2 reduced translational and rotational errors significantly, while values above 3 made the method more sensitive to wrong static matches. Static stereo gives metric scale, replaces DSO’s monocular initializer with a semi-dense depth map from NCC matching of $3 \times 5$ neighbourhoods along the horizontal epipolar line, initializes candidate depths, and, according to the authors, reduces sensitivity to large optical flow and rolling shutter. Point selection uses blocks proportional to the image size rather than square blocks of fixed size, which the authors found helpful for the wide KITTI images.

On the KITTI training sequences 00–10 (Table 1), with relative errors averaged over 100 m to 800 m intervals, Stereo DSO had a mean translational RMSE of 0.84% and rotational RMSE of 0.20° per 100 m, against 0.81% and 0.26° for ORB-SLAM2 run by the authors with loop closure and global bundle adjustment turned off, and 1.14% and 0.40° for Stereo LSD-SLAM’s odometry as reported in its paper. Rotational errors were better than ORB-SLAM2’s on every sequence and translational errors mixed, which the authors suggest may come from the dataset’s low frame rate. On the test sequences 11–21, the authors report Stereo DSO, as pure odometry, below the errors of ORB-SLAM2 and Stereo LSD-SLAM run as full SLAM systems with loop closure, at all path lengths and speeds (Figure 5).

LDSO

LDSO (Gao, Wang, Demmel and Cremers, 2018) turns DSO into a monocular SLAM system with loop closure. Only its abstract was read for this article; it describes the design as follows. The point selection favours corners, which keeps some points repeatable, so that loop candidates can be detected with a conventional bag-of-words approach. Candidates are verified geometrically, and a Sim(3) relative pose is estimated by jointly minimizing 2D and 3D geometric errors. These constraints are fused with a co-visibility graph of relative poses taken from DSO’s window optimization, and the resulting pose graph is optimized, without global bundle adjustment.

The official code shows how (tum-vision/LDSO, commit 9b0c48c). The default selection (setting_pointSelection = 1) divides the image into a grid, keeps the pixels of each cell whose squared gradient exceeds half the cell’s maximum (and at least 5), ranks them by Shi–Tomasi score and keeps the best few per cell. Of these, it marks as corners those scoring above 1% of the best score in the image and not within 5 pixels of a better one, and computes ORB-style descriptors for them. Setting 0 restores DSO’s selection, and the README notes that loop closing needs setting 1. Loop candidates are queried in a DBoW3 database, verified with OpenCV’s PnP with RANSAC, refined as a Sim(3) pose with g2o, and added as edges of a pose graph over all keyframes. The README stresses that LDSO still cannot recover the real scale; loop closure only makes it more consistent over long trajectories.

Other extensions

Visual-inertial. VI-DSO (von Stumberg, Usenko and Cremers, 2018) adds preintegrated IMU measurements between keyframes to DSO’s photometric energy and includes scale and gravity direction as variables. Because the scale is often not immediately observable from IMU data, the system starts with an arbitrary scale, and a strategy the authors call dynamic marginalization keeps partial marginalization consistent while the scale estimate is still far from its optimum.
Rolling shutter. Schubert et al. (2018) extend DSO with a rolling-shutter projection model, estimating a velocity at each keyframe under a constant-velocity prior, and report better results on rolling-shutter sequences than global-shutter VO.
Online photometric calibration. Bergmann, Wang and Cremers (2018) estimate exposure times, the response function and vignetting online from gain-robust KLT feature tracks, so that auto-exposure video from an uncalibrated camera can be processed. They report visual odometry precision on par with photometrically calibrated video.

The TUM project page lists further extensions, among them deep virtual stereo odometry and omnidirectional DSO (as of November 2026).