rotation of the cuboid chain

[lastchance]

@Alexas,
You still have not finally cleared up how your cubes and molecule are stored.

In your post http://www.cplusplus.com/forum/beginner/279658/#msg1208673
you give the following data for a polymer:

Polymer pb({ma, ma, ma, ma, ma, ma, ma, ma, ma}, {"poly2"}, {}, 10.0f, 15.0f, 5.0f,
               {{-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}});

Please explain what the first container is:

{ma, ma, ma, ma, ma, ma, ma, ma, ma}

What is each of those 9 elements, and what is the value of ma?


What do the succeeding 5 fields

{"poly2"}, {}, 10.0f, 15.0f, 5.0f,

represent? (You have sort of said that the last one is a "start point", but does it represent the coordinates of the centre or some other facet of the cube? What is the orientation of the first cube?


What exactly do the following lines, like

{-0.79948, 0.504813, 0.325568, 0.0888429}

represent? If they represent an axis and a rotation angle in radians then rotation with respect to what? (Remember that the first cube has no "previous cube".)


The problem is that you are failing to answer our questions.

[alexas]

Sorry, I am doing my best to answer your questions, just sometimes maybe I don't get it right.
I am sincerely grateful for your patience.

- every cuboid has a position dependent on the previous cube;

We have origin coordinates, and the first cube. We rotate it, and we calculate the position of the top center of it, xyz. We position the second cube that its bottom center is equal to previously calculated xyz. We rotate it and calculate the position of the top center of the second cube, and so on.

- every cuboid's rotation is described by its own axis-angle notion, i.e. it has a vector and angle. Important: the point of rotation is the bottom center point of each cube.
What does "bottom" mean?
If rotation is zero, is the cube always axis-aligned, or does its orientation depend on the previous cube? 

- I mean the center point of the cuboid bottom face
- If rotation is zero, the cube is always axis-aligned.

Please explain what the first container is:

It describes the geometry of each cuboid, so it means 9 cuboids are represented with set length/width/height

What do the succeeding 5 fields {"poly2"}, {}, 10.0f, 15.0f, 5.0f,

the first two are irrelevant, the next three are the coordinates of the first cuboid, its bottom face center.

What is the orientation of the first cube? 
What exactly do the following lines, like
{-0.79948, 0.504813, 0.325568, 0.0888429}
if they represent an axis and a rotation angle in radians then rotation with respect to what? 

Yes, they describe the axis and the angle in the respect to the cuboids bottom center. Please find the picture below with examples:
https://ibb.co/2h4j0HP

[lastchance]

What is the value of ma (the size of each cube)?

Are you using your own data structure for your polymer, or some pre-existing one? If the latter, then please provide a link to it directly, as it may make more sense.

Which is the "bottom face" - in the negative x direction, y direction or z direction?

[lastchance]

Let's guess that the rotation axis and angle of each cube correspond to what you would do to a cube initially oriented with the coordinate axes.

For the overall rotation matrix, compute its matrix M. I gave a formula in your earlier thread:
http://www.cplusplus.com/forum/beginner/278645/#msg1203056

For each segment in your polymer:
- compute the rotation matrix corresponding to that segment (same formula as above); call that P
- the resultant rotation matrix is R=MP (that order matters);
- if you want the axis (n1,n2,n3) and angle (theta) back from the resultant R, then you can get them from the trace and the antisymmetric part of R:
cos(theta) = (R11+R22+R33 - 1) / 2
n1.sin(theta) = (R32-R23) / 2
n2.sin(theta) = (R13-R31) / 2
n3.sin(theta) = (R21-R12) / 2
(This has two equivalent solutions, depending on which way you wish to orient your axis).
You could also use quaternions if you wish - but I have no experience of that.

For the start point of the first cube just apply the overall rotation matrix M to its displacement from whichever point you are rotating about.

For the start point of successive cubes calculate successive displacements from the previous one using the
individual segment rotations R. (You presumably had to do this to draw the first figure).

[alexas]

@lastchance,
Thank you very much for your patience and time, I really appreciate it.
Your method will work indeed in this case, I believe, and we will try to implement it.

The only drawback I can see is that we will have so much calculations for the long chains (10000 and more) - as the individual rotations P will have to be multiplied for the consecutive cuboids, i.e. rotation for 5th cuboid will be R = M*P1*P2*P3*P4

Please correct me if I didn't get right this.

[lastchance]

The only drawback I can see is that we will have so much calculations for the long chains (10000 and more) - as the individual rotations P will have to be multiplied for the consecutive cuboids, i.e. rotation for 5th cuboid will be R = M*P1*P2*P3*P4

You haven't made clear what the "rotation" for the n'th cuboid means. One of your previous pictures suggests it would be an actual rotation from an original cube oriented with the coordinate axes. In which case
R = M * P5
for the 5th cube. There is no reason to cumulate products of rotation matrices (although it's only a single additional multiply on the previous cube if it were; you wouldn't start the product again from scratch.)

I'm assuming that you drew the original molecule, so your code for doing so should indicate what the "rotation" of each cube means. We are still waiting for a definitive explanation from you.

[alexas]

I think it would be best to show you how it actually rotates. So I prepared the chain from the five cuboids with the following rotations carefully selected in such way to demonstrate how it actually works, please look to the picture. It should clearly indicate how exactly the rotations works for the cuboids. If there is still questions, please write them as simple as possible for me to understand them well. Thank you,

https://ibb.co/vVDcxhT

                  {0, 0, 1, 0},
                { 1, 0, 0, 1.3},
                {1, 0, 1, 1.0},
                { 0, 0, 1, 0},
                {1, 0, 1, 1.0},

[lastchance]

Well, from the two in that image marked
{0, 0, 1, 0}
which corresponds to zero angle rotation (i.e. nothing happens) it would appear that your cubes are rotated relative to ABSOLUTE space, not the previous cube.

[alexas]

yes, I confused you. What I meant is that cube's new position is based on the previous cube. The rotation is absolute then.

Sorry, all this topic was so wrongly communicated from my side, the only excuse is that I learned much more about all this - I had a bad understanding of the problem myself.

In this case, can your proposed solution still work?

[lastchance]

Yes, it works.
This code creates an image file (stl.stl). Open directly or in
https://www.viewstl.com/

The original polymer (pb1) is copied to pb2, which is then rotated 90 degrees about the origin.

EDIT: note forced revision of part of code in https://www.cplusplus.com/forum/beginner/279885/#msg1210032

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cmath>
using namespace std;

using Mat = vector< vector<double> >;
const double PI = 4.0 * atan( 1.0 );

//====

class Stl;

//---

class Shape
{
public:
   virtual void addToPlot( Stl &stl ) = 0;
};

//====

struct Pt
{
   double x, y, z;
   Pt( double x = 0, double y = 0, double z = 0 ) : x( x ), y( y ), z( z ) {}
};

//---

Pt operator + ( const Pt &a, const Pt &b ){ return { a.x + b.x, a.y + b.y, a.z + b.z }; }
Pt operator - ( const Pt &a, const Pt &b ){ return { a.x - b.x, a.y - b.y, a.z - b.z }; }
Pt operator / ( const Pt &a, double d    ){ return { a.x / d, a.y / d, a.z / d }; }
Pt operator * ( double d   , const Pt &a ){ return { d * a.x, d * a.y, d * a.z }; }
Pt cross      ( const Pt &a, const Pt &b ){ return { a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x }; }
double dot    ( const Pt &a, const Pt &b ){ return a.x * b.x + a.y * b.y + a.z * b.z; }
double normsq ( const Pt &a              ){ return dot( a, a ); }
double len    ( const Pt &a              ){ return sqrt( normsq( a ) ); }
ostream & operator << ( ostream &out, const Pt a ){ return out << a.x << " " << a.y << " " << a.z << " "; }

Pt operator * ( const Mat &M, const Pt &v )
{
   return { M[0][0] * v.x + M[0][1] * v.y + M[0][2] * v.z,
            M[1][0] * v.x + M[1][1] * v.y + M[1][2] * v.z,
            M[2][0] * v.x + M[2][1] * v.y + M[2][2] * v.z };
}

Mat operator * ( const Mat &M, const Mat &N )
{
   Mat R( M.size(), vector<double>( N[0].size(), 0.0 ) );
   for ( int i = 0; i < M.size(); i++ )
   {
      for ( int j = 0; j < N[0].size(); j++ )
      {
         for ( int k = 0; k < N.size(); k++ ) R[i][j] += M[i][k] * N[k][j];
      }
   }
   return R;
}

//====

Mat rotationMatrix( const Pt &axis, double radians )
{
   Mat result( 3, vector<double>( 3, 0.0 ) );
   Pt n = axis / len( axis );
   double cosA = cos( radians ), sinA = sin( radians ), cosA1 = 1.0 - cosA;
   result[0][0] = cosA + n.x * n.x * cosA1;
   result[0][1] =      + n.x * n.y * cosA1 - n.z * sinA;
   result[0][2] =      + n.x * n.z * cosA1 + n.y * sinA;
   result[1][0] =      + n.y * n.x * cosA1 + n.z * sinA;
   result[1][1] = cosA + n.y * n.y * cosA1;
   result[1][2] =      + n.y * n.z * cosA1 - n.x * sinA;
   result[2][0] =      + n.z * n.x * cosA1 - n.y * sinA;
   result[2][1] =      + n.z * n.y * cosA1 + n.x * sinA;
   result[2][2] = cosA + n.z * n.z * cosA1;
   return result;
}

//---

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 ) axis = Pt( 1.0, 0.0, 0.0 );
   else                     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 );
}

//====

struct Triangle
{
   Pt v1, v2, v3;
   Triangle( const Pt &v1, const Pt &v2, const Pt &v3 ) : v1( v1 ), v2( v2 ), v3( v3 ) {}
};

//====

class Stl
{
   vector<Triangle> triangles;
public:
   void add( Shape *S ) { S->addToPlot( *this ); }
   void addTriangle( const Pt &v1, const Pt &v2, const Pt &v3 ) { triangles.push_back( Triangle( v1, v2, v3 ) ); }
   void addRectangle( const Pt &v1, const Pt &v2, const Pt &v3, const Pt &v4 ) { addTriangle( v1, v2, v3 );   addTriangle( v1, v3, v4 ); }
   void insert( const Stl &stl ) { triangles.insert( triangles.end(), stl.triangles.begin(), stl.triangles.end() ); }
   void translate( const Pt &v ) { for ( Triangle &t : triangles ) t = Triangle( t.v1 + v, t.v2 + v, t.v3 + v ); }
   void scale( double s ) { for ( Triangle &t : triangles ) t = Triangle( s * t.v1, s * t.v2, s * t.v3 ); }
   void scale( double s, const Pt &centre ) { translate( -1 * centre );   scale( s );   translate( centre ); }
   void rotate( const Mat &M ) { for ( Triangle &t : triangles ) t = Triangle( M * t.v1, M * t.v2, M * t.v3 ); }
   void rotate( const Mat &M, const Pt &centre ) { translate( -1 * centre );   rotate( M );   translate( centre ); }
   void draw( const string &filename );
};

//---

void Stl::draw( const string &filename )
{
   ofstream out( filename );
   out << "solid\n";

   for ( Triangle t : triangles )
   {
      Pt n = cross( t.v2 - t.v1, t.v3 - t.v1 );
      n = n / len( n );

      out << "   facet normal " << n << '\n';
      out << "      outer loop\n";
      out << "         vertex " << t.v1 << '\n';
      out << "         vertex " << t.v2 << '\n';
      out << "         vertex " << t.v3 << '\n';
      out << "      endloop\n";
      out << "   endfacet\n";
   }

   out << "endsolid\n";
}

//====

class Cuboid : public Shape
{
   Pt centre;
   Pt side1, side2, side3;

public:
   Cuboid( const Pt &c, double L ) : centre( c ), side1( Pt( L, 0, 0 ) ), side2( Pt( 0, L, 0 ) ), side3( Pt( 0, 0, L ) ) {}
   Cuboid( const Pt &c, const Pt &s1, const Pt &s2, const Pt &s3 ) : centre( c ), side1( s1 ), side2( s2 ), side3( s3 ) {}

   void addToPlot( Stl &stl );
};

//---

void Cuboid::addToPlot( Stl &stl )
{
   Pt v1 = centre - 0.5 * ( side1 + side2 + side3 );
   Pt v2 = v1 + side1, v3 = v2 + side2, v4 = v3 - side1;
   Pt v5 = v1 + side3, v6 = v2 + side3, v7 = v3 + side3, v8 = v4 + side3;
   stl.addRectangle( v1, v2, v6, v5 );
   stl.addRectangle( v2, v3, v7, v6 );
   stl.addRectangle( v3, v4, v8, v7 );
   stl.addRectangle( v4, v1, v5, v8 );
   stl.addRectangle( v1, v4, v3, v2 );
   stl.addRectangle( v5, v6, v7, v8 );
}

//====

class Polymer
{
public:
   double side;
   Pt start;
   vector< pair<Pt,double> > cubes;

   Polymer( double side, const Pt &start, const vector< pair<Pt,double> > &cubes ) : side( side ), start( start ), cubes( cubes ) {}
   void rotate( const Pt &centre, const Pt &axis, double angle );
};

//---

void Polymer::rotate( const Pt &centre, const Pt &axis, double angle )
{
   Mat M = rotationMatrix( axis, angle );
   start = centre + M * ( start - centre );
   for ( auto &C : cubes )
   {
      Mat P = rotationMatrix( C.first, C.second );
      Mat R = M * P;
      axisAngle( R, C.first, C.second );
   }
}

//====

int main()
{
   Polymer pb1( 1.0,
                { 10.0, 15.0, 5.0 },
                { { { -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   } } );
   Polymer pb2 = pb1;
   Pt reference = { 0.0, 0.0, 0.0 };
   Pt axis = { 1.0, 0.0, 0.0 };
   double angle = PI / 2.0;
   pb2.rotate( reference, axis, angle );

   // Add both polymers to plot
   Stl stl;
   Pt ex = { 1, 0, 0 }, ey = { 0, 1, 0 }, ez = { 0, 0, 1 };
   for ( Polymer P : { pb1, pb2 } )
   {
      Pt start = P.start;
      for ( auto &C : P.cubes )
      {
         Mat R = rotationMatrix( C.first, C.second );
         Pt side1 = P.side * ( R * ex );
         Pt side2 = P.side * ( R * ey );
         Pt side3 = P.side * ( R * ez );
         Pt centre = start + 0.5 * side3;
         Cuboid cub( centre, side1, side2, side3 );
         stl.add( &cub );
         start = start + side3;
      }
   }

   stl.draw( "stl.stl" );
}

[alexas]

I am very, very grateful for your patience. I will start to implement this in my code right away and will come back to report how I managed.

[alexas]

Just small question, how exactly I can extract the new rotations in the axis-angle format from pb2? Using this formulas what you wrote previously?

- if you want the axis (n1,n2,n3) and angle (theta) back from the resultant R, then you can get them from the trace and the antisymmetric part of R:
cos(theta) = (R11+R22+R33 - 1) / 2
n1.sin(theta) = (R32-R23) / 4
n2.sin(theta) = (R13-R31) / 4
n3.sin(theta) = (R21-R12) / 4
(This has two equivalent solutions, depending on which way you wish to orient your axis).

[lastchance]

Look at the key routine Polymer::rotate. This is what actually rotates the polymer.

For each segment in the polymer it calls routine axisAngle() to get back from the 3x3 rotation matrix to the axis vector and angle form.

If you look at that routine (it's near the top of the code) then note that the indices in the code are 0-based rather than the 1-based set in my description earlier. This is simply because c++ arrays index from 0, not 1. So [2][1] would become [1][0] etc.

The shortage of comments in my post is because I struggled to stay within the character limit for a single post.

[alexas]

As my rotations are stored as exponential maps, I wrote such function. Exponential map uses float, so I hope the conversion will be ok. If not, I can maybe rewrite ExponentialMap in double. Also, if I understand correctly, you use "start" only for plotting later.

[lastchance]

Also, if I understand correctly, you use "start" only for plotting later.

No, "start" gives you the absolute first position of the molecule in space. You have to rotate it about the rotation centre like everything else.


I have no idea what connection your problem has with exponential maps.

[alexas]

It works. Finally. After so many frustration hours.
And just because of you.

I want to send you a bottle of a good rum or other good drink you like. Let me know please how I can do that.

Really. I appreciate this so much.

[lastchance]

Well, I'm glad that you got it to work. Made me think, anyway.

Some of those routines could prove useful elsewhere in future, so I'll forgo the alcoholic prize.

[alexas]

Just simple Thanks just does not sound seriously, I was totally lost in this problem for a long time and I was looking in a lot of place for help, but the only solution proposed was to use vertice coordinates with rotational matrixes which would force me to rewrite a lot of code which is already very well optimized for collision check and tested. Your approach was exactly what I needed.

I will need to test the speed of it for very long chains, but anyway, maybe I will be able to parallelize it as well.