X Y Z angles of a D20 (icosahedron) for gravity projection

Hi all,

EDIT#1: Added the .m script as an attachement

-------------------------- EDIT#2: Please read the comments --------------------------

Presentation (you can skip)

Thanks for your interest. Sorry for the title, it is a bit of a mouthful. Let me explain...

I have a large 3D printed D20 also named icosahedron. Basically, it has 20 faces. Please have a look at the attached picture :) .

I can attach a 3D angular sensor to any of the faces. The sensor would return me the 3D projection of the gravity like [0 0 -9.81] or [0 -4.905 -4.905].

What I would like to do

In Matlab, I would like to simulate this. To be more specific, given any of the 20 faces as a reference, I would like to get the 20 projections of the gravity (because of the 20 faces).

Let me explain this a bit more :)

Each time I "roll" the dice, I basically force one of the triangular faces to be "flat" (parallel to the XY plane). The reference face (to which my sensor is attached) is at a different angle. Therefore the gravity is not [0 0 -9.81] but a projection.

These projections are what I would like to calculate.

Therefore I think I expect to get a matrix like this:

projection#1 ---> [ [X Y Z];
projection#2 --->   [X Y Z];
...
projection#20 --->  [X Y Z] ]

What I can do

1. I have managed to get a list of all the D20 vertices.
2. With that, I can construct the 20 triangular faces of the icosahedron.
3. And I can get the 3D normal vectors of each to each of the 20 faces.

Where I am stuck

But the part where I am stuck is calculating the 3D angles (can I name them roll/pitch/yaw?) given 2 normals to the planes:

• the one which is flat
• the one on which the sensor is attached

Do I need a quaternion? a rotation matrix? What would be the best way to calculate the 20 orientations?

I tried something like the following but it only retuns a single angle:

 Angles_roll_pitch_yaw(cnt_face,:) = atan2(norm(...
         cross(refrence_3D_vector,Normals_Vector(cnt_face,:))), ...
         dot(refrence_3D_vector,Normals_Vector(cnt_face,:)));

Thank you very much for your time.

Sorry for the long post.

EDIT: I really can't understand if I need and what the function "vrrotvec" returns

Answer 1

[James Tursa]

If you have the 3D normal vectors, then all you need to do is take the dot product between them and the "flat" face normal vector. This will give you the cos(angle) between them which you can use to multiply by the normal vectors to get your gravity projection. E.g.,

flat_normal = the 3D normal vector of the flat face

other_normal = the 3D normal vector of another arbitrary face

g = gravity magnitude

flat_gravity = g * flat_normal

other_gravity = g * dot(flat_normal,other_normal) * other_normal

Is this what you are trying to do?

Side Note: You may need to use -g instead of g to get things to match your sensors output. Strictly speaking accelerometers measure reaction force, not gravity, and reaction force is (sort of) in the opposite direction from gravity. The caveat is because if you are on the ground you would strictly need to account for rotation of the Earth to get a more accurate reaction force.

[NE]

Thanks for your answer.

|Is this what you are trying to do?

Yes, despite my bad explanation you got it! Thanks!

Thanks for your side note. The "-" sign is noted for the raction force.

Also if I understood correctly

g = gravity magnitude

is a scalar?

And let say I have the data as:

other_normal = [other_normal_x other_normal_y other_normal_z]; 
flat_normal  = [flat_normal_x  flat_normal_y  flat_normal_z ]; 

So my 3D reaction force for the flat face is a scalar that is multiplied element-wise to a 3D vector (normal vector)?

projected_gravity_flat_face = g * [flat_normal_x flat_normal_y flat_normal_z];

and the one of the other faces is:

other_gravity = g * dot(flat_normal,other_normal) * other_normal

I'll try this method asap and l'll mark the answer accepted.

[James Tursa]

Yes g is a scalar.

[NE]

Thanks

[NE]

So let's make it simple and say I have 2 faces:

1. face#1 facing up (the flat one)
2. face#2 facing left (the other one)

Their respective normal vectors (3-by-1) would be:

flat_normal = [0 0 1]';
other_normal = [0 1 0]';

Let the gravity be:

g = -9.81
g_vector = [0 0 g]'; % X Y Z

I think using ENU is simplier in that case.

The gravity sensed if the sensor would be on the flat face would be:

projected_gravity_flat_face = [0 0 -g];

The gravity sensed if the sensor would be on the other face would be (we would expect the sensor to "see" the gravity on its X axis (instead of its Z)):

projected_gravity_other_face = [-g 0 0];

So from the formula I have:

projected_gravity_other_face = g * dot(flat_normal,other_normal) * other_normal
projected_gravity_other_face = 9.81 * dot([0 0 1]',[0 1 0]') * [0 1 0]'
projected_gravity_other_face = 9.81 * 0 * [0 1 0]'
projected_gravity_other_face =  0 * [0 1 0]'
projected_gravity_other_face =  [0 0 0]'

Or directly from command line:

>> 9.81 * dot([0 0 1]',[0 1 0]') * [0 1 0]'
ans =
     0
     0
     0

Should't there be at least 1 axis that senses gravity?

Sorry if I made any mistakes.

[NE]

Hi James!

I think I found a way (but might not be the best, maybe Matlab has some function that already does that?).

Is this the way you would do it?

So what I need apparently is a rotation matrix that changes a vector from one frame to another. So I have to calculate 19 matrices (because 1 of them is the flat/refrence one).

There is a guy that did just that on a Github gist. Here is the source.

https://gist.github.com/kevinmoran/b45980723e53edeb8a5a43c49f134724

It looks like C but I put it to Matlab easily (below).

PS: if you think that is too easy you can look at the hardcore version here: https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d

% v1 is the reference (flat), v2 is the new (the other one)
% I replace those with my normal vectors
v1 = [1 0 0]';
v2 = [0 0 1]';

cross_p = cross( v1, v2 );

cosA = dot( v1, v2 );
k    = 1 / (1 + cosA);

rotmat_1_to_2 = [...
    [ (cross_p(indx_X) * cross_p(indx_X) * k) + cosA,                ...
      (cross_p(indx_Y) * cross_p(indx_X) * k) - cross_p(indx_Z),  	 ...
      (cross_p(indx_Z) * cross_p(indx_X) * k) + cross_p(indx_Y)];  	 ...
    [ (cross_p(indx_X) * cross_p(indx_Y) * k) + cross_p(indx_Z),  	 ...
      (cross_p(indx_Y) * cross_p(indx_Y) * k) + cosA,                ...
      (cross_p(indx_Z) * cross_p(indx_Y) * k) - cross_p(indx_X)];	 ...
    [ (cross_p(indx_X) * cross_p(indx_Z) * k) - cross_p(indx_Y),	 ...
      (cross_p(indx_Y) * cross_p(indx_Z) * k) + cross_p(indx_X),	 ...
      (cross_p(indx_Z) * cross_p(indx_Z) * k) + cosA]                ...
];

Let's check! The answer of each line should be 1 (logical).

% Check

%1) obtain v2 from v1 and rotmat_1_to_2
all(   (rotmat_1_to_2) * v1 == v2   )

%2) obtain v1 from v2 and inv(rotmat_1_to_2)
all(inv(rotmat_1_to_2) * v2 == v1   )
% or more efficient
all(   (rotmat_1_to_2) \ v2 == v1   )

[James Tursa]

OK, I may have misunderstood your original question. I thought you were asking how much gravity (or reaction force) would be sensed along the normal direction to each face, and that is what I provided. I.e., if you had a single sensor mounted on that surface that only sensed acceleration along the normal direction to that surface.

But now with your talk of rotation matrices it sounds like you have a 3-axis accelerometer package attached to each face, so you need 3D mapping to each face package. Is that what you actually have? If so, then I don't see how just the normal vector directions gives you enough information to recover the rotation matrix. I.e., this sounds more like an alignment between two coordinate frames that you are trying to calculate. Is that it?

[NE]

quick note: I find the lack of "quote" formatting quite disturbing on this site, but I'll use "||"

Sorry, I should have explained it better. I have 10 IMUs, so I can cover 1/2 of the faces at once.

To find all the possible accelerations (1 per face), I have 2 choices, I think.

• Choice1: I make an acquisition, then roll the dice only once, then make a second acquisition.
• Choice2: Or I can just make an acquisition, then unmount the 10 IMUs, mount them to the 10 other faces without moving the D20, and then make my second acquisition.

||Is that what you actually have?

This is what I have. I should have spent more time on the hardware description (this is my first post).

||If so, then I don't see how just the normal vector directions gives you enough information to recover the rotation matrix.

From the formula from the gist, you apparently need only 2x 3D vectors v1 and v2 (the normals in my case), and it gives me logical results for the gravity projection, so I am quite pleased with that simple approach.

||I.e., this sounds more like an alignment between two coordinate frames that you are trying to calculate. Is that it?

I wasn't sure at the beginning, but with the gist, I am sure that the alignment is at least one (of the possibly other) ways to do it.

So to describe the process:

1. Get the vertices coordinates
2. Get the vertice face map (what are the 3 vertices per face)
3. Calculate normal vector to each plane
4. Create rotation matrix (3-by-3) from a reference face or any other plane, with their normal vectors (see the gist method)
5. Use the rotation matrix to project the gravity vector (originally expressed in the reference face) on the new face.
6. Enjoy the algorithm

Note: this might not be the most elegant and memory-efficient way because I need to store a 3-by-3-by-20 (or 19) matrix for all the rotation matrices but ¯\_(ツ)_/¯

Final question: should I put my final Matlab script in the attachment for the other to see or not?

Thank you so much for your time and your perseverance.