8

I've seen several similar questions, and have some ideas of what I might try, but I don't remember seeing anything about spread.

So: I am working on a measurement system, ultimately computer vision based.

I take N captures, and process them using a library which outputs pose estimations in the form of 4x4 affine transformation matrices of translation and rotation.

There's some noise in these pose estimations. The standard deviation in Euler angles for each axis of rotation is less than 2.5 degrees, so all orientations are pretty close to each other (for a case where all Euler angles are close to 0 or 180). Standard errors of less than 0.25 degrees are important to me. But I have already run into the problems endemic to Euler angles.

I want to average all these pretty-close-together pose estimates to get a single final pose estimate. And I also want to find some measure of spread so that I can estimate accuracy.

I'm aware that "average" isn't actually well defined for rotations.

(For the record, my code is in Numpy-heavy Python.)

I also may want to weight this average, since some captures (and some axes) are known to be more accurate than others.

My impression is that I can just take the mean and standard deviation of the translation vector, and that for the rotation I can convert to quaternions, take the mean, and re-normalize with OK accuracy since these quaternions are pretty close together.

I've also heard mentions of least-squares across all the quaternions, but most of my research into how this would be implemented has been a dismal failure.

Is this workable? Is there a reasonably well-defined measure of spread in this context?

Peter O.
  • 32,158
  • 14
  • 82
  • 96
ikrase
  • 335
  • 3
  • 13
  • I'd rather do the averaging on the matrix coefficients than on the Euler angles. –  Jul 12 '15 at 17:17
  • [Here](https://stackoverflow.com/q/12374087/1264582) is a related question about averaging quaternions with several answers containing code snippets in C++, Matlab, Python. – ketza Jan 25 '22 at 17:10

2 Answers2

3

Without more info about your geometry setup is hard to answer. Anyway for rotations I would:

  1. create 3 unit vectors

    x=(1,0,0),y=(0,1,0),z=(0,0,1)

    and apply the rotation on them and call the output

    x(i),y(i),z(i)

    it is just applying the matrix(i) with position at (0,0,0)

  2. do this for all measurements you have

  3. now average all vectors

    X=avg(x(1),x(2),...x(n))
Y=avg(y(1),y(2),...y(n))
Z=avg(z(1),z(2),...z(n))

  4. correct the vector values

    so make each of the X,Y,Z unit vectors again and take the axis which is more closest to the rotation axis as main axis. It will stay as is and recompute the remaining two axises as cross product of main axis and the other vector to ensure orthogonality. Beware of the multiplication order (wrong order of operands will negate the output)

  5. construct averaged transform matrix

    see transform matrix anatomy as origin you can use averaged origin of the measurement matrices

Spektre
  • 49,595
  • 11
  • 110
  • 380
  • I'm inclined to like this method, since it would seem closely linked to the physical rotation of an object with several features, thus averaging the positions of the features (and would also allow yielding standard deviation / standard error / variance / IQR with a well defined representation) – ikrase Jul 13 '15 at 16:57
  • What information about the geometry setup do you need? – ikrase Jul 13 '15 at 17:16
  • 1
    @ikrase maybe a sketch of your skeleton or whatever you are measuring pose of, and few example inputs. btw you can get the rotated X,Y,Z axises from your transform matrix directly (or from its inverse matrix depends on your matrix pipeline and representation) – Spektre Jul 13 '15 at 18:06
  • Actually, I'm a little confused by what you are doing in Step 4. You start with the axis that's closest to rotation axis (for me this will normally be Y), and then just take that one, then you take the second axis temporarily, calculate the third axis as 1st X 2nd, and then get the true 2nd as 3rd x 1st (or something, didn't check the commutativity there)? Am I getting that right? – ikrase Jul 13 '15 at 22:04
  • @ikrase the step 4 is just making your computed axises vectors parallel to each other. If you take any 2 non parallel vectors `a,b` and make cross product of them then you get parallel vector to booth of them. For any 2 non parallel vectors there are just 2 such vectors (in 3D) one is negation of the other `c=a x b` and `-c= a x b` So when you exploit this you can easily correct 2 axises from the main one to ensure orthogonality. So you are right you just need to set the order of multiplication correctly otherwise you will have some mirrored axises. So if you do just swap the cross operands – Spektre Jul 14 '15 at 06:52
  • 1
    @ikrase The cross product also changes the size of output vector but if both are unit then the resultant vector is also unit in size – Spektre Jul 14 '15 at 06:53
  • Is there a reason why not to just use SVD or other method of orthogonalizing the matrix of vectors directly? – ikrase Jul 14 '15 at 16:18
  • 1
    @ikrase I am used to pure C++ where I have not such things at my disposal and this is pretty direct way anyway... of coarse you can use any other method you know or can ... – Spektre Jul 14 '15 at 20:17
  • @Spektre I am not sure what you mean by `it is just applying the matrix(i) with position at (0,0,0)`. Do you mean take each of the rotation matrices and multiply each of $x(i), y(i), z(i)$ and get the corresponding rotated axes vectors? – Luca Nov 23 '18 at 16:54
  • @Luca if you multiply transform matrix offseted to `(0,0,0)` by a unit basis vector like `(1,0,0)` the result will be vector parallel to the matrix coordinate system x axis ... Take a look at: [Understanding 4x4 homogenous transform matrices](https://stackoverflow.com/a/28084380/2521214)... Its the same like extracting the basis vectors from the matrix directly ... – Spektre Nov 23 '18 at 17:10
1

Moakher wrote a paper that explains there are basically two ways to take an average of Rotation matrices. The first is a weighted average followed by a projection back to SO(3) using the SVD. The second is the Riemannian center of mass. That one is a closer notion to the geometric mean, and its more complicated to compute.

Charles F
  • 539
  • 6
  • 11