Is it logical wanting to combine CGAL with Quaternions

Question

I am creating a Discrete Element Method simulation program and I am using CGAL to describe the polyhedrons. From reading literature I was planning to do my differential equations for rotation with Quaternions due to the better numerical stability and lack of gimbal lock. However CGAL does not seem to support rotation based on quaternions. (Please tell me if I am incorrect here) I find it a bit surprising that this seems to be missing, certainly since CGAL likes to be absolute in its accuracy which seems to fit well with the numerical stability of quaternions.

Question: Can I somehow combine Boost Quaternions with CGAL or is there any easy way to implement this. And if so, would this be a logical idea to try?

The other options I think I have are:

• writing my differential equations for the affine rotation used is CGAL and deal with the downsides there.
• store the orientation as an affine rotation matrix and convert it to Quaternions and use this in the diff. equations. Obviously I am worried about the needed conversion step here every timestep.

Any suggestions or other options that I might think of are greatly appreciated.

Answer 1

First Option: Use the Aff_transformation_3 class

CGAL does not provide a quaternion class, it does provide the Aff_transformation_3 class though. Which you could easily use like this:

CGAL::Surface_mesh<Kernel> P;
std::transform( P.points_begin(), P.points_end(), P.points_begin(), yourAffineTransformation);

for defining the transformation matrix see this.

Second Option: Use Quaternions

If you want to use quaternions you would need to construct one with an external library. For example you could use Eigen:

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> //or whichever kernel suits your needs
#include <CGAL/Surface_mesh.h>
#include <Eigen/Geometry>

using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Polyhedron = CGAL::Surface_mesh<Kernel>;
using Point = CGAL::Point_3<Kernel>;

// define the function that rotates your mesh
template <typename Vect, typename Quaternion>
void rotateCGALPolyhedron(Polyhedron P, Vect to_rotation_center,
                          Quaternion quat) {
  for (auto vi : P.vertices()) {
    Point p = P.point(vi);
    // translate your point to the rotation center. In your case this would be
    // the center of mass of the Polyhderon
    Vect V(p[0] - to_rotation_center[0], p[1] - to_rotation_center[1],
           p[2] - to_rotation_center[2]);
    // construct the translation vector that moves your point to the rotated
    // position
    Vect v = quat * V; //the Vect operator*(Quaternion, Vect) must be implemented!! If you use Eigen::Quaternion you could use Eigen::Vector3d
    // retranslate the point back to its initial position and translate it using
    // the previously created translation vector
    P.point(size_t(vi)) =
        Point(to_rotation_center[0] + v[0], to_rotation_center[1] + v[1],
              to_rotation_center[2] + v[2]);
  }
}

int main() {

  // define your rotation using eigen's quaternion class
  Eigen::Quaternion<double> quad(..);
  Eigen::Vector_3d centerOfMass; //find the center of mass of the mesh you want to rotate
  rotateCGALPolyhedron(P.vertices.begin(), P.vertices.end(), centerOfMass,
                       quad);

return 0;
}

As you can see since cgal does not have an implementation for quaternions if you want to use quaternions the code is lengthy compared to the Aff_transformation_3 case.