cuboid chain rotation - reversed case

[alexas]

I have a chain of cuboids. Every cube's new position is based on the previous cube. The cubes are rotated relative to ABSOLUTE space using axis-angle notion.

Details in this thread, where the rotation was solved with the help of @lastchance : https://www.cplusplus.com/forum/beginner/279658/2/

I'm currently stuck on the reverse rotation, which means I want to make the first cube the final one, and vice versa.

- I am using std::reverse to reverse the rotation vector;

If the chain is relatively simple, I can multiply the angle by -1 and get what I need. However, if I have complex shape, I cannot do that (somehow it destroys the original shape), and simply reversing the rotation gives me like mirrored image, see picture: https://ibb.co/Hg1cHv3

How could I properly reverse the order of rotations/cuboids to preserve the exactly same structure?

[lastchance]

It depends how your chain of cubes is constructed.

In my code here:
https://www.cplusplus.com/forum/beginner/279658/2/#msg1208888
the chain was constructed by offsetting successive cubes by (side length times) the rotation of the z unit vector. The sides of the resulting cubes were (side length times) the rotation of the x,y,z unit vectors (look in int main()).

In the original scheme the offset in the chain goes from the bottom to top of each cube. If you reverse the order of storage then you will need to go from top to bottom. Hence, for each cube, you will need to reverse the z unit vector first (e.g. by rotating PI radians about the y axis) and THEN apply the rotation.

The following routine will do it in my code above. Note that, to do the reversing, the M matrix must be applied to the unit vector first, and hence it must be the SECOND term in the matrix product.

void Polymer::backward()
{
   Mat M = rotationMatrix( Pt{ 0.0, 1.0, 0.0}, PI );   // This will rotate unit cube first, before main rotation
   reverse( cubes.begin(), cubes.end() );              // Cubes now listed in reverse order
   for ( auto &C : cubes )
   {
      Mat P = rotationMatrix( C.first, C.second );     // Convert to rotation matrix form
      Mat R = P * M;                                   // Order matters!!!
      axisAngle( R, C.first, C.second );               // Convert back to axis/angle form
   }
}

Note also:-
- If your chain is constructed in a different way then you will have to do whatever is necessary to your own code;
- This fixes the orientation of the whole structure; as it stands, the end point of the reversed structure will coincide with the start point of the original and you will have to translate the new start point if you want to do anything other than that. In my code, just
pb2.start = pb2.start - Pt{ 5.0, 0.0, 0.0 };
would do.
- This will only work satisfactorily for objects (like cubes) where it doesn't matter which side is the rotation of x and which the rotation of y; it would be more complicated for cross-sections without two planes of symmetry as you would probably end up having to flip from a right-handed to a left-handed coordinate system.

[alexas]

Thank you very much, you are helping immensely.
One (hopefully last) more question just to get deeper in this,

If I have one cube with position x1, y1, z1 and rotation {rx1,ry1,rz1, theta1}, and another cube with position x2, y2, z2 and rotation {rx2,ry2,rz2, theta2}, what rotation should be applied to the second cube, that its rotation would become {rx1,ry1,rz1, theta1+ 0.5}, for example?

[lastchance]

alexas wrote:

another cube with position x2, y2, z2 and rotation {rx2,ry2,rz2, theta2}, what rotation should be applied to the second cube, that its rotation would become {rx1,ry1,rz1, theta1+ 0.5

In matrix terms, if you wanted to change a rotation of R₁ into one of R₂ then you would apply the matrix
R₂.R₁^-1
to the current matrix.

The inverse of a rotation matrix is easy to compute because it is an orthogonal matrix: its transpose is its inverse.

[alexas]

Thank you again - it works indeed. I am learning a lot about the rotation matrixes!

[alexas]

Small update regarding reversing. I constructed the following function and it works perfectly with the complex case - but apparently it does not work with simple cases. Could you maybe see the mistake? Or maybe I did not get smth.

The function:

void
Polymer::reverse()
{
    Mat M = rotationMatrix( Pt{ 0.0, 1.0, 0.0}, getPiValue());
    std::reverse(_monomers.begin(), _monomers.end());
    std::reverse(_rotations.begin(), _rotations.end());
    for (int32_t i = 0; i < number_of_monomers(); i++)
    {
        float rx, ry, rz, theta;
        _rotations[i].get_rotation(rx, ry, rz, theta);
        Pt monomer_axis = {rx, ry, rz};
        Mat P = rotationMatrix(monomer_axis, theta);
        Mat R = P * M;
        axisAngle(R, monomer_axis, theta);
        ExponentialMap mr(monomer_axis.x, monomer_axis.y, monomer_axis.z, theta);
        set_rotation(mr, i);
        fill_coords(i);
    }
}

works fine with the first rotation chain, but not with the second:

//works well
{ { { -0.79948  ,  0.504813 , 0.325568 }, 0.0888429 },
                  { { -0.941624 ,  0.0857279, 0.325568 }, 6.06663   },
                  { { -0.940697 , -0.0953601, 0.325568 }, 0.454249  },
                  { {  0.945361 , -0.017274 , 0.325568 }, 0.932644  },
                  { {  0.936006 , -0.133785 , 0.325568 }, 0.585938  },
                  { {  0.0287705,  0.945081 , 0.325568 }, 6.16575   },
                  { { -0.190289 , -0.926173 , 0.325568 }, 1.48524   },
                  { { -0.922358 , -0.207991 , 0.325568 }, 2.83318   },
                  { { -0.0299494,  0.945044 , 0.325568 }, 2.60525   } } );

//does not work
                {{0.0f, 0.0f, 1.0f, 0.0f},
                {0.2f, 0.0f, 0.7f, 1.2f},
                {0.1f, 0.0f, 0.8f, 1.4f},
                {0.3f, 0.0f, 0.9f, 1.6f},
                {0.5f, 0.0f, 0.5f, 1.7f}}
//outcome of the reversed latter:
1	0	0	3.14159
1	0	0	3.14159
1	0	0	3.14159
1	0	0	3.14159
1	0	0	3.14159

[lastchance]

Mmmm, you've managed to hit the only edge case where the conversion from rotation matrix to axis-angle form fails - when sin(angle) = 0.

The case when angle=0 doesn't matter: you could use any axis you like if the rotation angle is zero. However, you've managed to highlight the case when angle = PI.

You can try the following revised code for the conversion back to axis-angle form. As you see, it has increased in length just to deal with a rare edge case. (Your second example then runs in my code).

void axisAngle( const Mat &R, Pt &axis, double &angle )
{
   angle = acos( 0.5 * ( R[0][0] + R[1][1] + R[2][2] - 1 ) );
   double s = sin( angle );
   if ( abs( s ) < 1.0e-6 )       // Edge case sin(theta)=0: painful!
   {
      // Find largest diagonal element and use that column
      int i = 0;
      if ( R[1][1] > R[i][i] ) i = 1;
      if ( R[2][2] > R[i][i] ) i = 2;
      axis = Pt( R[0][i], R[1][i], R[2][i] );
      if ( i == 0 ) axis.x += 1;
      if ( i == 1 ) axis.y += 1;
      if ( i == 2 ) axis.z += 1;
      // Note: axis must be normalised below; it currently has length 2n_i
   }
   else
   {
      // Antisymmetric part should determine axis without problems
      axis = ( 0.5 / s ) * Pt( R[2][1] - R[1][2], R[0][2] - R[2][0], R[1][0] - R[0][1] );
   }
   axis = axis / len( axis );     // Normalise to length 1
}

For the record the problem is as follows. The rotation matrix is

   R[0][0] = cosA + n.x * n.x * (1-cosA)
   R[0][1] =      + n.x * n.y * (1-cosA) - n.z * sinA
   R[0][2] =      + n.x * n.z * (1-cosA) + n.y * sinA
   R[1][0] =      + n.y * n.x * (1-cosA) + n.z * sinA
   R[1][1] = cosA + n.y * n.y * (1-cosA)
   R[1][2] =      + n.y * n.z * (1-cosA) - n.x * sinA
   R[2][0] =      + n.z * n.x * (1-cosA) - n.y * sinA
   R[2][1] =      + n.z * n.y * (1-cosA) + n.x * sinA
   R[2][2] = cosA + n.z * n.z * (1-cosA)

In most cases you can use the antisymmetric part to extract the axis (nx,ny,nz), with subsequent normalisation removing the factor (2.sinA):

       R[2][1] - R[1][2] =  (2.sinA) nx
       R[0][2] - R[2][0] =  (2.sinA) ny
       R[1][0] - R[0][1] =  (2.sinA) nz

The only place where that doesn't work is if sinA = 0, in which case you have to try to extract the components of the axis vector from one of the columns (with cosA being either 1 or -1). The safest way (since up to two of nx, ny or nz could be zero) is to use the one with the largest diagonal element.

[alexas]

Thank you very much.