ideas for rotation of multiple cuboids

[alexas]

Hello,
I am asking such maybe more mathematical question here in this forum as it somehow gives me the best ideas in my project. I apologize if this is not acceptable.

I have the cuboids which are rotated at each other top point, see the picture:
https://ibb.co/XsDTTm7

The rotations are described as AxisAngle structures {x, y, z, theta}.

I can also using conversions to the quaternions if needed.
And now I need to rotate the whole structure according its origin points, let's say 45 degrees. I cannot simply multiply the existing rotations by new one as the original structure will be lost. Only the initial cuboid should be rotated by 45 degrees, and the rest of them must must remain in the same structure but rotated as a whole. So I guess the directions of rotation vectors must be recalculated. However, I have no clue currently how to start it or where to look for the information. Could someone please advise?

[lastchance]

Write down the matrix for the given rotation (about the origin). See below.
Multiply all coordinates by that matrix.

If the overall rotation isn't about the origin then translate the rotation point to the origin, do the rotation, translate it back again.


Rotation Matrix

If you have a rotation of angle θ about a unit vector (n_i) passing through the origin, then the rotated coordinates are given by (with summation convention)
M_ijx_j
and the rotation matrix (M_ij) is (with summation convention, δ_ij is Kronecker delta, ε_ijk is alternating tensor):
M_ij=δ_ijcos θ+n_in_j(1-cos θ)-ε_ijkn_ksin θ

(That formula can be worked out by splitting a position vector into components parallel and perpendicular to n, then rotating the latter part.)

[alexas]

if I understood you correctly, then you just rotate each cube individually.
I need to rotate them as a whole which means that every cube should get each own individual rotation matrix

[lastchance]

alexas wrote:

if I understood you correctly,

You didn't.

[alexas]

Oh I see, sorry, you are saying to multiply the coordinates by the rotation matrix, not the existing rotations. Sorry and thank you for your input.

Unfortunately my mathematic background is too low to understand that formula. Do you mean the standard rotation matrix like that? https://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToMatrix/

Thank you,

[lastchance]

You can find it written out element by element at
https://en.m.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
The elements correspond exactly to the one-line formula, although I derived it from vector geometry and I think the derivation in Wikipedia is grossly over-complicated.

[alexas]

no, it does not work. It is my fault that probably I cannot properly explain what I want.

I dont need to rotate the individual cuboid. I need to rotate the whole structure about its point of origin. The cube can be represented simply as a point.

I need somehow rotate the structure, not the individual cubes - they must remain in the same position in regards with other cubes as they were.

So I am rotating points by utilising axis/angle notation:

ExponentialMap::ExponentialMap(const float coord_x,
                               const float coord_y,
                               const float coord_z,
                               const float theta)
    : _coord_x(coord_x)

    , _coord_y(coord_y)

    , _coord_z(coord_z)

    , _theta(theta)
{
    _normalize();
}

void
ExponentialMap::rotate (float& x,
                       float& y,
                       float& z,
                       const float rx,
                       const float ry,
                       const float rz) const
{
    const float tsin = std::sin(_theta);

    const float tcos = std::cos(_theta);

    x = rx * tcos;

    y = ry * tcos;

    z = rz * tcos;

    const float vx = _coord_y * rz - _coord_z * ry;

    const float vy = _coord_z * rx - _coord_x * rz;

    const float vz = _coord_x * ry  - _coord_y * rx;

    x += tsin * vx;

    y += tsin * vy;

    z += tsin * vz;

    const float fp = _coord_x * rx + _coord_y * ry + _coord_z * rz;

    const float tfact = (1.0f - tcos) * fp;

    x += tfact * _coord_x;

    y += tfact * _coord_y;

    z += tfact * _coord_z;
}

//creating random rotation

ExponentialMap
StructureGenerator:: _random_rotation() const
{
    const auto m_pi = getPiValue();

    float x = cos(19.0f * m_pi/180);

    float z = sin(19.0f * m_pi/180);

    float y = 0.0f;

    auto prng = RandomGenerator(std::chrono::high_resolution_clock::now().time_since_epoch().count());

    auto r = ExponentialMap(0.0f, 0.0f, 1.0f, prng.get_float(0.0f, 2 * m_pi));

    float refx, refy, refz;

    r.rotate(refx, refy, refz, x, y, z);

    return ExponentialMap(refx, refy, refz, prng.get_float(0.0f, 2 * m_pi));
}

[jonnin]

rotating the whole thing is the same as rotating each one about a point. not in place, but about a point, possibly the origin or moving it as if the origin.

[alexas]

How?? I really cannot get it. It simply does not work when I multiply the points by the new rotation matrix. Every point then it rotates per new angle, but they are losing their original positions in respect with each other.

[keskiverto]

If you do multiply every point with the same rotation matrix, then they all should rotate around same axis the same angle and maintain the point-to-point distances.

[jonnin]

then your new RM is wrong. I forget how to do them, always have to look it up...
page is kinda spammy but this is a simple breakdown of it all:
https://www.gatevidyalay.com/3d-rotation-in-computer-graphics-definition-examples/

[keskiverto]

Given point P (x y z) and matrix R

a b c
d e f
g h i

The rotated point P' (xn yn zn) = R*P

xn = a*x + b*y + c*z;
yn = d*x + e*y + f*z;
zn = g*x + h*y + i*z;

The trick is to create the R correctly.

Jonnin's link shows matrices for "rotate around X", "... Y", and "... Z". That is usually possible, but can lead to "gimbal lock".

Alternative is quaternion. It is possible to create R from quaternion:
https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Quaternion-derived_rotation_matrix

It is also possible to create quaternion from axis of rotation and angle.
Close to AxisAngle, but one is not the other.

[jonnin]

if its just visual graphics, you can kick the offending values just a little bit, eg 0.01 radians, and it will clear up, you just have to throw in a detection for near zero angles (if I remember it right) and do the fudgery. If you need it *precise*, then quaternions are the way to go. In a scene in visual graphics, no one is going to notice something being rotated a nanodegree or two extra if the distances are smallish. If you are rotating pluto around the sun, you may notice it, and need more precision.

quaternions are amazing, but a little heavy on the trig function calls. What they do, in nontechnical non mathbabble, is simply generate a 4th dimension of redundant information so when the 3-d information fails the 4th dimension still has the 'lost' information so you can 'correct' the problem. (This is all done automatically in the provided equations).

[lastchance]

Well, you can try this code @alexas. It creates an STL image file (stl.stl). You can probably open that directly in Windows 10. If not then just use an online viewer like
https://www.viewstl.com/

The code creates a line of cubes and cylinders, then rotates them. Have a play.

#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 Vec
{
   double x, y, z;
   Vec( double x = 0, double y = 0, double z = 0 ) : x( x ), y( y ), z( z ) {}
};

//----------------------------

Vec operator + ( const Vec &a, const Vec &b ){ return { a.x + b.x, a.y + b.y, a.z + b.z }; }
Vec operator - ( const Vec &a, const Vec &b ){ return { a.x - b.x, a.y - b.y, a.z - b.z }; }
Vec operator / ( const Vec &a, double d     ){ return { a.x / d, a.y / d, a.z / d }; }
Vec operator * ( double d    , const Vec &a ){ return { d * a.x, d * a.y, d * a.z }; }
Vec cross      ( const Vec &a, const Vec &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 Vec &a, const Vec &b ){ return a.x * b.x + a.y * b.y + a.z * b.z; }
double normsq  ( const Vec &a               ){ return dot( a, a ); }
double len     ( const Vec &a               ){ return sqrt( normsq( a ) ); }
ostream & operator << ( ostream &out, const Vec a ){ return out << a.x << " " << a.y << " " << a.z << " "; }

Vec operator * ( const Mat &M, const Vec &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 rotationMatrix( const Vec &axis, double radians )
{
   Mat result( 3, vector<double>( 3, 0.0 ) );
   Vec 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;
}

//======================================================================

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

//======================================================================

class Stl
{
   vector<Triangle> triangles;
public:
   void add( Shape *S ) { S->addToPlot( *this ); }
   void addTriangle( const Vec &v1, const Vec &v2, const Vec &v3 ) { triangles.push_back( Triangle( v1, v2, v3 ) ); }
   void addRectangle( const Vec &v1, const Vec &v2, const Vec &v3, const Vec &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 Vec &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 Vec &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 Vec &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 )
   {
      Vec n = cross( t.v2 - t.v1, t.v3 - t.v1 );
      n = n / len( n );                // unit normal vector;

      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
{
   Vec centre;                         // centre
   Vec side1, side2, side3;            // side vectors

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

   void addToPlot( Stl &stl );
};

//----------------------------

void Cuboid::addToPlot( Stl &stl )
{
   Vec v1 = centre - 0.5 * ( side1 + side2 + side3 );
   Vec v2 = v1 + side1, v3 = v2 + side2, v4 = v3 - side1;
   Vec 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 Cylinder : public Shape
{
   Vec centre;                                             // centre (NOT base centre)
   Vec side1, side2, side3;                                // vectors along radius, along radius perpendicular, along axis
   int nface;                                              // number of rectangular faces

public:
   Cylinder( const Vec &c, double r, double h, int n = 30 )
      : centre( c ), side1( Vec( r, 0, 0 ) ), side2( Vec( 0, r, 0 ) ), side3( Vec( 0, 0, h ) ), nface( n ) {}
   Cylinder( const Vec &c, const Vec &s1, const Vec &s3, int n = 30 )
      : centre( c ), side1( s1 ), side2( cross( s3, s1 ) / len( s3 ) ), side3( s3 ), nface( n ) {}

   void addToPlot( Stl &stl );
};

//----------------------------

void Cylinder::addToPlot( Stl &stl )
{
   Vec bottom = centre - 0.5 * side3;                      // centre of base
   Vec top    = bottom + side3;                            // centre of top

   double dtheta = 2.0 * PI / nface;
   Vec v1, v2, v3, v4;
   v2 = bottom + side1;
   v3 = v2 + side3;
   for ( int n = 1; n <= nface; n++ )
   {
      double theta = n * dtheta;
      v1 = v2;
      v4 = v3;
      v2 = bottom + cos( theta ) * side1 + sin( theta ) * side2;
      v3 = v2 + side3;
      stl.addRectangle( v1, v2, v3, v4 );                  // add sides as a series of rectangles
      stl.addTriangle( v2, v1, bottom );                   // add triangles for bottom
      stl.addTriangle( v4, v3, top    );                   // add triangles for top
   }
}

//======================================================================

int main()
{
   Stl stl;
   for ( int i = 1; i <= 4; i++ )                // Create a line of cylinders and cubes
   {
      Cylinder cyl( 2 * i, 0.4, 0.6 );
      Cuboid cub( Vec( 2 * i - 1, 0, 0 ), 0.6 );
      stl.add( &cyl );
      stl.add( &cub );
   }

   // Duplicate and rotate
   Stl extra = stl;
   Vec axis( 0.0, 0.0, 1.0 );
   double degrees = 135.0;
   Mat R = rotationMatrix( axis, degrees * PI / 180.0 );
   extra.rotate( R );

   // Add to original
   stl.insert( extra );

   // Output to file
   stl.draw( "stl.stl" );
}