File: differencing.org

package info (click to toggle)
mrcal 2.5-3
  • links: PTS, VCS
  • area: main
  • in suites: forky
  • size: 8,992 kB
  • sloc: python: 40,651; ansic: 15,632; cpp: 1,754; perl: 303; makefile: 160; sh: 99; lisp: 84
file content (632 lines) | stat: -rw-r--r-- 30,186 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
#+TITLE: Model differencing
#+OPTIONS: toc:t

* Intrinsics differences
mrcal provides the [[file:mrcal-show-projection-diff.html][=mrcal-show-projection-diff=]] tool to compute and display
projection differences between several models (implemented using
[[file:mrcal-python-api-reference.html#-show_projection_diff][=mrcal.show_projection_diff()=]] and [[file:mrcal-python-api-reference.html#-projection_diff][=mrcal.projection_diff()=]]). This has numerous
applications. For instance:

- evaluating the manufacturing variation of different lenses
- quantifying intrinsics drift due to mechanical or thermal stresses
- testing different solution methods
- underlying a cross-validation scheme

** What is meant by a "difference"?
We want to compare different representations of the same lens, so we're /not/
interested in extrinsics. At a very high level, to evaluate the projection
difference at a pixel coordinate $\vec q_0$ in camera 0 we need to:

1. Unproject $\vec q_0$ to a fixed point $\vec p$ using lens 0
2. Project $\vec p$ back to pixel coords $\vec q_1$ using lens 1
3. Report the reprojection difference $\vec q_1 - \vec q_0$

[[file:figures/diff-notransform.svg]]

This simple definition is conceptually sound, but isn't applicable in practice.
In the [[file:tour-differencing.org][tour of mrcal]] we calibrated the same lens using the same data, but with
two different lens models. The models are describing the same lens, so we would
expect a low difference. However, the simple algorithm above produces a
difference that is non-sensically high. As a heat map:

#+begin_src sh
mrcal-show-projection-diff \
  --intrinsics-only        \
  --cbmax 15               \
  --unset key              \
  opencv8.cameramodel      \
  splined.cameramodel
#+end_src
#+begin_src sh :exports none :eval no-export
# THIS IS GENERATED IN tour-differencing.org
#+end_src

[[file:external/figures/diff/diff-radius0-heatmap-splined-opencv8.png]]

And as a vector field:

#+begin_src sh
mrcal-show-projection-diff \
  --intrinsics-only        \
  --vectorfield            \
  --vectorscale 5          \
  --gridn 30 20            \
  --cbmax 15               \
  --unset key              \
  opencv8.cameramodel      \
  splined.cameramodel
#+end_src
#+begin_src sh :exports none :eval no-export
# THIS IS GENERATED IN tour-differencing.org
#+end_src

[[file:external/figures/diff/diff-radius0-vectorfield-splined-opencv8.svg]]

The reported differences are in pixels.

The issue is similar to the one encountered by the [[file:uncertainty.org::#propagating-through-projection][projection uncertainty]]
routine: each calibration produces noisy estimates of all the intrinsics and all
the coordinate transformations:

[[file:figures/uncertainty.svg]]

The above plots projected the same $\vec p$ in the camera coordinate system, but
that coordinate system shifts with each model computation. So in the /fixed/
coordinate system attached to the camera housing, we weren't in fact projecting
the same point.

There exists some transformation between the camera coordinate system from the
solution and the coordinate system defined by the physical camera housing. It is
important to note that *this implied transformation is built-in to the
intrinsics*. Even if we're not explicitly optimizing the camera pose, this
implied transformation is still something that exists and moves around in
response to noise.

The above vector field suggests that we need to pitch one of the cameras. We can
automate this by adding a critical missing step to the procedure above between
steps 1 and 2:

- Transform $\vec p$ from the coordinate system of one camera to the coordinate
  system of the other camera

[[file:figures/diff-yestransform.svg]]

We don't know anything about the physical coordinate system of either camera, so
we do the best we can: we compute a fit. The "right" transformation will
transform $\vec p$ in such a way that the reported mismatches in $\vec q$ will
be small. Previously we passed =--intrinsics-only= to bypass this fit. Let's
omit that option to get the the diff that we expect:

#+begin_src sh
mrcal-show-projection-diff \
  --unset key              \
  opencv8.cameramodel      \
  splined.cameramodel
#+end_src
#+begin_src sh :exports none :eval no-export
# THIS IS GENERATED IN tour-differencing.org
#+end_src

[[file:external/figures/diff/diff-splined-opencv8.png]]

** Implied transformation details
:PROPERTIES:
:CUSTOM_ID: implied-transformation
:END:

As with [[file:uncertainty.org::#effect-of-range][projection uncertainty]], the difference computations are not invariant to
range. So we always compute "the projection difference when looking out to $r$
meters" for some possibly-infinite $r$. The procedure we implement is:

1. Regularly sample the two imagers, to get two corresponding sets of pixel
   coordinates $\left\{\vec q_{0_i}\right\}$ and $\left\{\vec q_{1_i}\right\}$

2. Unproject the camera-0 pixel coordinates to a set of points $\left\{\vec
   p_{0_i}\right\}$ in the camera-0 coordinate system. The range is given with
   =--distance=. [[file:mrcal-python-api-reference.html#-sample_imager_unproject][=mrcal.sample_imager_unproject()=]] function does this exactly
3. Unproject the camera-1 pixel coordiantes to a set of normalized unit
   observation vectors $\left\{\vec v_{1_i}\right\}$ in the camera-1 coordinate
   system. These correspond to $\left\{\vec p_{0_i}\right\}$, except they're
   unit vectors instead of points in space.
4. Compute the implied transformation $\left(R,t\right)$ as the one to maximize
   \[ \sum_i w_i \vec v_{1_i}^T \frac{R \vec p_{0_i} + t}{\left|R \vec p_{0_i} +
   t\right|} \] where $\left\{w_i\right\}$ is a set of [[#fit-weighting][weights]]. As with main
   calibration optimization, this one is unconstrained, using the
   [[file:conventions.org::#pose-representation][=rt= transformation
   representation]]. The inner product above is $\cos \theta$ where $\theta$ is
   the angle between the two observation vectors in the camera-1 coord system.

When looking out to infinity the $t$ becomes insignificant, and we do a
rotation-only optimization.

This is the logic behind [[file:mrcal-python-api-reference.html#-implied_Rt10__from_unprojections][=mrcal.implied_Rt10__from_unprojections()=]] and
[[file:mrcal-python-api-reference.html#-projection_diff][=mrcal.projection_diff()=]].

*** Which transform is "right"?
We just described a differencing method that computes an implied transformation
given a range $r$. This will produce a different result for each $r$, but in
reality, there's only a single /true/ transformation, and the solutions at
different ranges are different estimates of it.

If we just need to compare two different representations of the same lens, then
we don't care about the implied transformation itself. The correct thing to do
here would be to set $r$ to the intended working distance of the system: for me
this is generally infinity. Looking at a single $r$, these
implied-transformation fits will always overfit a little bit, but from
experience, this doesn't affect the resulting model-difference output.

However, in some applications ([[#extrinsics-diff][extrinsics differencing]] for instance), we /do/
want the true implied transform. And when we look at this computed transform
more deeply, we see that with some $r$ we get a result that is clearly wrong.
For instance:

#+begin_example
$ mrcal-show-projection-diff \
    --no-uncertainties       \
    --unset key              \
    --distance 1000          \
    opencv8.cameramodel      \
    splined.cameramodel

Transformation cam1 <-- cam0:  rotation: 0.316 degrees, translation: [-2.67878241 -0.832383    1.8063095 ] m
## WARNING: fitted camera moved by 3.336m. This is probably aphysically high, and something is wrong. Pass both a high and a low --distance?
#+end_example

More on =--no-uncertainties= later; it's here to speed things up, and isn't
important. We asked for a diff at 1000m out, and the solver said that the
optimal transform moves the camera coordinate system back by 1.8m and to the
right by 2.7m. This is looking at the same data as before: comparing two solves
from the [[file:tour-differencing.org][tour of mrcal]]. Nothing moved. The camera coordinate system could have
shifted inside the camera housing a tiny bit, but the solved shifts are huge,
and clearly aren't inside the housing anymore. This needs a deeper investigation
on how to do it "properly", but for practical use I have a working solution:
solve using points at two ranges at the same time, a near range and a far range:

#+begin_example
$ mrcal-show-projection-diff \
    --no-uncertainties       \
    --unset key              \
    --distance 1,1000        \
    opencv8.cameramodel      \
    splined.cameramodel

Transformation cam1 <-- cam0:  rotation: 0.316 degrees, translation: [2.90213542e-05 2.59871714e-06 1.73397909e-03] m
#+end_example

Much better. This is analogous to the multiple ranges-to-the camera we need when
we compute a [[file:recipes.org::#surveyed-calibration][surveyed calibration]] and the multiple ranges we get when we [[file:tour-choreography.org::#tilt][tilt
the chessboard during calibration]].

*** Selection of fitting data
:PROPERTIES:
:CUSTOM_ID: fitting-data-selection
:END:

When we use a fit to compute the implied transformation, we minimize the
reprojection error. We hope that the main contributions to this error come from
geometric misalignment. If this was the case, minimizing this reprojection error
would produce the correct implied $\left(R,t\right)$. Applying this
transformation would correct the misalignment, leaving a very small residual. It
is possible to have other error contributions though, which would break this
method of finding the implied $\left(R,t\right)$. This is common when intrinsics
differences are significant; let's demonstrate that.

We just saw a difference result from the [[file:tour-differencing.org][tour of mrcal]], showing that the two
models are similar but not identical. Once again:

#+begin_src sh
mrcal-show-projection-diff \
  --unset key              \
  opencv8.cameramodel      \
  splined.cameramodel
#+end_src
#+begin_src sh :exports none :eval no-export
# THIS IS GENERATED IN tour-differencing.org
#+end_src

[[file:external/figures/diff/diff-splined-opencv8.png]]

This uses the default data-selection behavior of [[file:mrcal-show-projection-diff.html][=mrcal-show-projection-diff=]]:

- [[#fit-weighting][uncertainties are used to weight the sampled points]]
- points are sampled from the /whole/ imager (without uncertainties the default
  behavior is to use a limited region in the center instead)

What if we used all the data in the imager, but weighed them all equally?

#+begin_src sh
mrcal-show-projection-diff \
  --unset key              \
  --no-uncertainties       \
  --radius 1e6             \
  --cbmax 50               \
  opencv8.cameramodel      \
  splined.cameramodel
#+end_src
#+begin_src sh :exports none :eval no-export
D=~/projects/mrcal-doc-external/2022-11-05--dtla-overpass--samyang--alpha7/2-f22-infinity/
mrcal-show-projection-diff           \
  --unset key                                         \
  --no-uncertainties       \
  --radius 1e6             \
  --cbmax 50               \
  $D/{opencv8,splined}.cameramodel         \
  --hardcopy ~/projects/mrcal-doc-external/figures/diff/diff-splined-opencv8--all-data.png \
  --terminal 'pngcairo size 1024,768 transparent noenhanced crop font ",12"'
#+end_src

[[file:external/figures/diff/diff-splined-opencv8--all-data.png]]

Note the completely different color-bar scaling: this difference is huge! The
reason this doesn't work is that projections of the two models don't just differ
due to an implied transformation: there's also a significant difference in
intrinsics. This difference is small near the center, and huge in the corners,
so even the optimal $\left(R,t\right)$ fits badly. We can solve this problem by
either

- [[#fit-weighting][Giving the bad points low weights during the solve]] by using the default
  behavior, without =--no-uncertainties=. This is described in the next section
- Excluding the bad points entirely by passing a =--radius=

Let's demonstrate the effect of =--radius=. I re-ran the calibration from the
[[file:tour.org][tour of mrcal]] using [[file:lensmodels.org::#lensmodel-opencv][=LENSMODEL_OPENCV4=]]. This model fits the wide lens even
worse than the [[file:lensmodels.org::#lensmodel-opencv][=LENSMODEL_OPENCV8=]] fits we studied, and is not expected to work
at all. But it is very good for this demo. If we run this solve, the [[file:formulation.org::#outlier-rejection][outlier
rejection]] logic kicks in, and makes the solve converge despite the model not
fitting the data well. The resulting model is available [[file:external/2022-11-05--dtla-overpass--samyang--alpha7/2-f22-infinity/opencv4.cameramodel][here]]. Using an
insufficient model like this isn't a good way to run calibrations: the outlier
rejection will throw away the clearly-ill-fitting measurements, but marginal
measurements will make it through, which will slightly poison the result. I'm
only doing this for this demo. If you really do want to use an insufficient
model, and you really do care only about points near the center, you should use
the [[file:mrcal-cull-corners.html][=mrcal-cull-corners=]] tool to throw away the data you don't care about, prior
to calibrating.

Let's compute the diff between the narrow-only =LENSMODEL_OPENCV4= lens model
and the mostly-good-everywhere =LENSMODEL_OPENCV8= model, using an expanding
radius of points. We expect this to work well when using a small radius, and we
expect the difference to degrade as we use more and more data away from the
center.

#+begin_src sh
# This is a zsh loop
for r (200 600 800 1000) {
  mrcal-show-projection-diff \
    --no-uncertainties       \
    --distance 1             \
    --radius $r              \
    --cbmax 15               \
    --unset key              \
    --extratitle "radius $r" \
    opencv[48].cameramodel
}
#+end_src
#+begin_src sh :exports none :eval no-export
D=~/projects/mrcal-doc-external/2022-11-05--dtla-overpass--samyang--alpha7/2-f22-infinity/
for r (200 600 800 1000) {
  ~/projects/mrcal/mrcal-show-projection-diff                    \
    --no-uncertainties                                           \
    --distance 1                                                 \
    --radius $r                                                  \
    --cbmax 15                                                   \
    --unset key                                                  \
    --extratitle "radius $r"                                     \
    $D/opencv[48].cameramodel                                    \
    --hardcopy ~/projects/mrcal-doc-external/figures/diff/diff-radius$r-opencv4-opencv8.png \
    --terminal 'pngcairo size 1024,768 transparent noenhanced crop font ",12"'
}
#+end_src

[[file:external/figures/diff/diff-radius200-opencv4-opencv8.png]]

[[file:external/figures/diff/diff-radius600-opencv4-opencv8.png]]

[[file:external/figures/diff/diff-radius800-opencv4-opencv8.png]]

[[file:external/figures/diff/diff-radius1000-opencv4-opencv8.png]]

All of these radii are equally "right", and there's a trade-off in picking one.
The models agree very well in a small region at the center, but they can also
agree decently well in a much larger region, if we are willing to accept a
higher level of error in the middle.

For the purposes of computing the implied transform you don't need a lot of
data, so I generally use a small region in the center, where I have reasonable
confidence that the intrinsics match. For a camera like this one =--radius 500=
is usually plenty.

*** Fit weighting
:PROPERTIES:
:CUSTOM_ID: fit-weighting
:END:

Clearly the =LENSMODEL_OPENCV4= solve does agree with the =LENSMODEL_OPENCV8=
solve well, but /only/ in the center of the imager. The issue from a tooling
standpoint is that in order for the tool to tell us that, *we* needed to tell
the tool to only look at the center. That is unideal.

The issue we observed is that some regions of the imager have unreliable
behavior, which poisons the fit. But we know where the fit is reliable: in the
areas where the [[file:uncertainty.org][projection uncertainty]] is low. So we can weigh the fit by the
inverse of the projection uncertainty, and we will then favor the "good"
regions. Without requiring the user to specify the good-projection region.

This works, but with a big caveat. As described on the [[file:uncertainty.org][projection uncertainty]]
page, lean models report overly-optimistic uncertainties. Thus when used as
weights for the fit, areas that actually are unreliable will be weighted too
highly, and will still poison the fit. We see that here, when comparing the
=LENSMODEL_OPENCV4= and =LENSMODEL_OPENCV8= results. The above plots show that
the =LENSMODEL_OPENCV4= result is only reliable within a few 100s of pixels
around the center. However, =LENSMODEL_OPENCV4= is a very lean model, so its
uncertainty at 1m out (near the sweet spot, where the chessboards were) looks
/far/ better than that:

#+begin_src sh
mrcal-show-projection-uncertainty \
    --distance 1                  \
    --unset key                   \
    opencv4.cameramodel
#+end_src
#+begin_src sh :exports none :eval no-export
D=~/projects/mrcal-doc-external/2022-11-05--dtla-overpass--samyang--alpha7/2-f22-infinity/
~/projects/mrcal/mrcal-show-projection-uncertainty               \
    --distance 1                                                 \
    --unset key                                                  \
    $D/opencv4.cameramodel                                       \
    --hardcopy ~/projects/mrcal-doc-external/figures/uncertainty/uncertainty-opencv4-1m.png \
    --terminal 'pngcairo size 1024,768 transparent noenhanced crop font ",12"'
#+end_src

[[file:external/figures/uncertainty/uncertainty-opencv4-1m.png]]

And the diff using that uncertainty as a weight /without/ specifying a radius
looks poor:

#+begin_src sh
mrcal-show-projection-diff \
  --distance 1             \
  --unset key              \
  opencv[48].cameramodel
#+end_src
#+begin_src sh :exports none :eval no-export
D=~/projects/mrcal-doc-external/2022-11-05--dtla-overpass--samyang--alpha7/2-f22-infinity/
~/projects/mrcal/mrcal-show-projection-diff                    \
    --distance 1                                               \
    --unset key                                                \
    $D/opencv[48].cameramodel                       \
  --hardcopy ~/projects/mrcal-doc-external/figures/diff/diff-weighted-opencv4-opencv8.png \
  --terminal 'pngcairo size 1024,768 transparent noenhanced crop font ",12"'
#+end_src

[[file:external/figures/diff/diff-weighted-opencv4-opencv8.png]]

Where this technique /does/ work well is when using [[file:splined-models.org][splined models]], which
produce realistic uncertainty estimates. To demonstrate, let's produce a
splined-model calibration that is only reliable in a particular region of the
imager. We do this by culling the [[file:tour.org][tour of mrcal]] calibration data to throw out
all points outside of a circle at the center, calibrate off /that/ data, and run
a diff on /those/ results:

#+begin_src sh
< corners.vnl \
  mrcal-cull-corners --imagersize 6000 3376 --cull-rad-off-center 1500 \
> /tmp/raw.vnl &&
  vnl-join --vnl-sort - -j filename /tmp/raw.vnl \
     <(< /tmp/raw.vnl vnl-filter -p filename --has level |
                       vnl-uniq -c |
                       vnl-filter 'count > 20' -p filename ) \
  > corners-rad1500.vnl

mrcal-calibrate-cameras                                                         \
  --corners-cache corners-rad1500.vnl                                           \
  --lensmodel LENSMODEL_SPLINED_STEREOGRAPHIC_order=3_Nx=30_Ny=18_fov_x_deg=150 \
  --focal 1900                                                                  \
  --object-spacing 58.8e-3                                                      \
  --object-width-n 14                                                           \
  '*.JPG'

mrcal-show-projection-uncertainty \
  --distance 1                    \
  --unset key                     \
  splined-rad1500.cameramodel

mrcal-show-projection-diff \
  --distance 1             \
  --unset key              \
  splined{,-rad1500}.cameramodel
#+end_src
#+begin_src sh :exports none :eval no-export
D=~/projects/mrcal-doc-external/2022-11-05--dtla-overpass--samyang--alpha7/2-f22-infinity/

< $D/corners.vnl \
  mrcal-cull-corners --imagersize 6000 3376 --cull-rad-off-center 1500 \
> /tmp/raw.vnl &&
  vnl-join --vnl-sort - -j filename /tmp/raw.vnl \
     <(< /tmp/raw.vnl vnl-filter -p filename --has level |
                       vnl-uniq -c |
                       vnl-filter 'count > 20' -p filename ) \
  > $D/corners-rad1500.vnl

~/projects/mrcal/mrcal-calibrate-cameras               \
  --corners-cache $D/corners-rad1500.vnl \
  --lensmodel LENSMODEL_SPLINED_STEREOGRAPHIC_order=3_Nx=30_Ny=18_fov_x_deg=150 \
  --imagersize 6000 3376              \
  --focal 1900                        \
  --object-spacing 58.8e-3            \
  --object-width-n 14                 \
  --out /tmp                          \
  '*.JPG'

mv /tmp/camera-0.cameramodel $D/splined-rad1500.cameramodel

~/projects/mrcal/mrcal-show-projection-uncertainty                     \
  --distance 1                                                         \
  --unset key                                                          \
  $D/splined-rad1500.cameramodel                                       \
  --hardcopy ~/projects/mrcal-doc-external/figures/uncertainty/uncertainty-splined-rad1500-1m.png \
  --terminal 'pngcairo size 1024,768 transparent noenhanced crop font ",12"'

~/projects/mrcal/mrcal-show-projection-diff                            \
  --distance 1                                                         \
  --unset key                                                          \
  $D/splined{,-rad1500}.cameramodel                                    \
  --hardcopy ~/projects/mrcal-doc-external/figures/diff/diff-weighted-splined-rad1500.png \
  --terminal 'pngcairo size 1024,768 transparent noenhanced crop font ",12"'
#+end_src

The cut-down corners are [[file:external/2022-11-05--dtla-overpass--samyang--alpha7/2-f22-infinity/corners-rad1500.vnl][here]] and the resulting model is [[file:external/2022-11-05--dtla-overpass--samyang--alpha7/2-f22-infinity/splined-rad1500.cameramodel][here]]. The uncertainty
of this model looks like this:

[[file:external/figures/uncertainty/uncertainty-splined-rad1500-1m.png]]

and the diff like this:

[[file:external/figures/diff/diff-weighted-splined-rad1500.png]]

This is yet another reason to use only splined models for real-world lens
modeling.

We just demonstrated that projection uncertainties provide a working method for
selecting the points used to fit the implied transformation. From a practical
standpoint, there's a downside: computing uncertainties makes the differencing
method much slower. So while the uncertainty-aware method is available in
[[file:mrcal-show-projection-diff.html][=mrcal-show-projection-diff=]], most of the time I run something like
=mrcal-show-projection-diff --no-uncertainties --radius 500=. This falls back on
simply using the points in a 500-pixel-radius circle at the center. This
computes much faster, and with good-coverage-low-error calibrations usually
produces a very similar result. If something looks odd, I will run the
uncertainty-aware diff to debug. But my usual default is =--no-uncertainties=.

* Extrinsics differences
:PROPERTIES:
:CUSTOM_ID: extrinsics-diff
:END:

The above technique can be used to quantify intrinsics differences between two
representations of the same lens. This is useful, for instance, to detect any
intrinsics calibration drift over time. This doesn't tell us anything about
extrinsics drift, however. It's possible to have a multi-camera system composed
of very stable lenses mounted on a not-rigid-enough mount. Any mechanical stress
wouldn't affect the intrinsics, but the extrinsics /would/ shift. And we want a
method to quantify this shift .

** Extrinsics only
I haven't needed to do this very often, so the technique I'm using isn't mature
yet. The extrinsics diff computation is implemented in the
[[https://github.com/dkogan/mrcal/blob/master/analyses/extrinsics-stability.py][=extrinsics-stability.py= tool]]. This isn't stable yet, and only exists in the
mrcal sources for now. But I used it several times, and it appears to work well.
Here I'm describing an extrinsics-only method: I'm computing and applying an
implied transform to cancel out any intrincs shift, leaving the effect of
extrinsics only.

In this description I will consider 2-camera systems, but the approach is
general for any number of cameras. Let's say we have two calibrations (0 and 1)
of a stereo pair ($\mathrm{left}$ and $\mathrm{right}$ cameras). Between the
calibrations the system was stressed (shaked, flipped, heated, etc), and we want
to know if the camera geometry shifted as a result. The obvious technique is to
compare the transformation $T_{0\mathrm{left},0\mathrm{right}}$ and
$T_{1\mathrm{left},1\mathrm{right}}$. Each of these is available directly in the
two calibrations, and we can compute the difference
$T_{0\mathrm{right},1\mathrm{right}} = T_{0\mathrm{right},0\mathrm{left}}
T_{1\mathrm{left},1\mathrm{right}}$. This is the transform between the right
cameras in the two calibrations if we line up the two left cameras.

This approach sounds good, but it is incomplete because it ignores the
transformation implied by the different intrinsics, as described above. So
lining up the two left cameras does not line up their projection behavior. And
comparing the poses of the right cameras does not compare their projection
behavior.

We just talked about how to compute the implied transforms
${T^\mathrm{implied}_{0\mathrm{left},1\mathrm{left}}}$ and
$T^\mathrm{implied}_{0\mathrm{right},1\mathrm{right}}$. So the final
transformation describing the shift of the right camera is

\[
T^\mathrm{implied}_{1\mathrm{right},0\mathrm{right}}
T_{0\mathrm{right},0\mathrm{left}}
T^\mathrm{implied}_{0\mathrm{left},1\mathrm{left}}
T_{1\mathrm{left},1\mathrm{right}}
\]

The [[https://github.com/dkogan/mrcal/blob/master/analyses/extrinsics-stability.py][=extrinsics-stability.py= tool]] implements this logic. To compute the implied
transformations we want the "true" transform, not a transform at any particular
range, so we use a near and a far range.

I tried this out in practice on a physical camera pair. I had calibrations
before and after a lot of mechanical jostling happened. And for each calibration
I had an odd and even set for [[file:tour-cross-validation.org][cross-validation]], which reported very low
intrinsics differences: the intrinsics were "correct". I evaluated the
extrinsics stability, looking at the odd calibrations:

#+begin_example
$ ~/projects/mrcal/analyses/extrinsics-stability.py \
  BEFORE/camera-0-odd-SPLINED.cameramodel           \
  BEFORE/camera-1-odd-SPLINED.cameramodel           \
  AFTER/camera-0-odd-SPLINED.cameramodel            \
  AFTER/camera-1-odd-SPLINED.cameramodel

translation: 0.04mm in the direction [0.13 0.06 0.99]
rotation:    0.01deg around the axis [ 0.93  0.32 -0.18]
#+end_example

So it claims that the right camera shifted by 0.04mm and yawed by 0.01deg. This
sounds low-enough to be noise, but what /is/ the noise level? The tool also
reports the camera resolution for comparison against the reported rotation:

#+begin_example
Camera 0 has a resolution of 0.056 degrees per pixel at the center
Camera 1 has a resolution of 0.057 degrees per pixel at the center
Camera 2 has a resolution of 0.056 degrees per pixel at the center
Camera 3 has a resolution of 0.057 degrees per pixel at the center
#+end_example

So this rotation is far smaller than one pixel. What if we look at the
allegedly-identical "even" calibration:

#+begin_example
$ ~/projects/mrcal/analyses/extrinsics-stability.py \
  BEFORE/camera-0-even-SPLINED.cameramodel          \
  BEFORE/camera-1-even-SPLINED.cameramodel          \
  AFTER/camera-0-even-SPLINED.cameramodel           \
  AFTER/camera-1-even-SPLINED.cameramodel

translation: 0.07mm in the direction [ 0.63 -0.73  0.25]
rotation:    0.00deg around the axis [ 0.55 -0.32  0.77]
#+end_example

That's similarly close to 0. What if instead of comparing before/after we
compare odd/even before and then odd/even after? odd/even happened at the same
time, so there was no actual shift, and the result will be at the noise floor.

#+begin_example
$ ~/projects/mrcal/analyses/extrinsics-stability.py \
  BEFORE/camera-0-odd-SPLINED.cameramodel           \
  BEFORE/camera-1-odd-SPLINED.cameramodel           \
  BEFORE/camera-0-even-SPLINED.cameramodel          \
  BEFORE/camera-1-even-SPLINED.cameramodel

translation: 0.07mm in the direction [0.09 0.76 0.64]
rotation:    0.00deg around the axis [ 0.77 -0.1  -0.63]


$ ~/projects/mrcal/analyses/extrinsics-stability.py \
  AFTER/camera-0-odd-SPLINED.cameramodel            \
  AFTER/camera-1-odd-SPLINED.cameramodel            \
  AFTER/camera-0-even-SPLINED.cameramodel           \
  AFTER/camera-1-even-SPLINED.cameramodel

translation: 0.05mm in the direction [9.69e-01 8.13e-04 2.46e-01]
rotation:    0.00deg around the axis [-0.43 -0.78  0.44]
#+end_example

So with these calibrations we have strong evidence that no extrinsics drift has
occurred. More testing and development of this method are planned, and this tool
will be further documented and released.

** Joint intrinsics+extrinsics
There's also an experimental [[file:mrcal-show-stereo-pair-diff.html][=mrcal-show-stereo-pair-diff=]] tool that shows the
pixel shift (not geometry shift) of the intrinsics and extrinsics together. This
works, but more is needed to figure out when a intrinsics-only, extrinsics-only
or joint diff is appropriate. More tooling is needed to support all the cases
and to report the joint and extrinsics-only uncertainties