Changeset 2cd0a0 for src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Timestamp:

Sep 10, 2016, 4:14:02 PM (9 years ago)

Author:

Frederik Heber <heber@…>

Branches:

SaturateAtoms_singleDegree

Children:

226860

Parents:

408415

git-author:

Frederik Heber <heber@…> (06/05/14 17:42:39)

git-committer:

Frederik Heber <heber@…> (09/10/16 16:14:02)

Message:

Dropped quaternion rotation for simple arbitrary rotation axis and angle.

• with newCenter, oldCenter and the cross product we have all we need: a rotation axis and an angle. We don't need to burden ourselves with those stupid, absolutely not working quaternions.
• removeMatchingPoints() now works on an array.
• Orientation rotation was wrong way round, added check.
• TESTFIX: removed QuaternionTest from SphericalPointDistributionTest, marked FragmentMolecule-cycle and StoreSaturatedFragment regression tests as XFAIL.

File:

• src/Fragmentation/Exporters/SphericalPointDistribution.cpp (modified) (4 diffs)

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -139 +139 @@
 
 SphericalPointDistribution::Polygon_t removeMatchingPoints(
-    const SphericalPointDistribution::Polygon_t &_returnpolygon,
+    const VectorArray_t &_points,
     const IndexList_t &_matchingindices
     )
@@ -147 +147 @@
   std::sort(indices.begin(), indices.end());
   LOG(4, "DEBUG: sorted matching is " << indices);
-  IndexArray_t::const_iterator valueiter = indices.begin();
-  SphericalPointDistribution::Polygon_t::const_iterator pointiter =
-      _returnpolygon.begin();
-  for (unsigned int i=0; i< _returnpolygon.size(); ++i, ++pointiter) {
-    // skip all those in values
-    if (*valueiter == i)
-      ++valueiter;
-    else
-      remainingreturnpolygon.push_back(*pointiter);
+  IndexArray_t remainingindices(_points.size(), -1);
+  std::generate(remainingindices.begin(), remainingindices.end(), UniqueNumber);
+  IndexArray_t::iterator remainiter = std::set_difference(
+      remainingindices.begin(), remainingindices.end(),
+      indices.begin(), indices.end(),
+      remainingindices.begin());
+  remainingindices.erase(remainiter, remainingindices.end());
+  LOG(4, "DEBUG: remaining indices are " << remainingindices);
+  for (IndexArray_t::const_iterator iter = remainingindices.begin();
+      iter != remainingindices.end(); ++iter) {
+    remainingreturnpolygon.push_back(_points[*iter]);
   }
-  LOG(4, "DEBUG: remaining indices are " << remainingreturnpolygon);
 
   return remainingreturnpolygon;
@@ -292 +293 @@
     // determine rotation angles to align the two point distributions with
     // respect to bestmatching
-    SphericalPointDistribution::Polygon_t rotated_newpolygon;
+    VectorArray_t rotated_newpolygon = remainingnew;
     Vector oldCenter;
     {
@@ -309 +310 @@
       if ((oldCenter - newCenter).NormSquared() > std::numeric_limits<double>::epsilon()*1e4) {
         // setup quaternion
-        Vector RotationAxis = newCenter;
-        RotationAxis.VectorProduct(oldCenter);
+        Vector RotationAxis = oldCenter;
+        RotationAxis.VectorProduct(newCenter);
+        Line Axis(zeroVec, RotationAxis);
         RotationAxis.Normalize();
-        const double RotationAngle = oldCenter.Angle(newCenter)/(M_PI/2.);
-//            RotationAxis.Angle(oldCenter) - RotationAxis.Angle(newCenter);
-        boost::math::quaternion<double> q
-            (RotationAngle, RotationAxis[0], RotationAxis[1], RotationAxis[2]);
-        LOG(5, "DEBUG: RotationAxis is " << RotationAxis
-            << ", RotationAngle is " << RotationAngle);
-        LOG(5, "DEBUG: Quaternion describing rotation is " << q);
+        const double RotationAngle = oldCenter.Angle(newCenter); // /(M_PI/2.);
+        LOG(5, "DEBUG: Rotate coordinate system by " << RotationAngle
+            << " around axis " << RotationAxis);
+
+        // apply rotation angles
+        for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
+            iter != rotated_newpolygon.end(); ++iter) {
+          Vector &current = *iter;
+          LOG(5, "DEBUG: Original point is " << current);
+          current =  Axis.rotateVector(current, RotationAngle);
+          LOG(5, "DEBUG: Rotated point is " << current);
+        }
+      }
+    }
+    // rotate triangle/line/point around itself to match orientation
+    if (MCS.bestmatching.size() > 1) {
+      if (oldCenter.NormSquared() > std::numeric_limits<double>::epsilon()*1e4) {
+        // construct RotationAxis and two points on its plane, defining the angle
+        const Line RotationAxis(zeroVec, oldCenter);
+        Vector oldPosition(rotated_newpolygon[*MCS.bestmatching.begin()]);
+        oldPosition.ProjectOntoPlane(RotationAxis.getDirection());
+        Vector newPosition(remainingold[*MCS.bestmatching.begin()]);
+        newPosition.ProjectOntoPlane(RotationAxis.getDirection());
+
+        // construct reference vector to determine direction of rotation
+        Vector dreiBein(oldPosition);
+        dreiBein.VectorProduct(oldCenter);
+        dreiBein.Normalize();
+        const double sign =
+            (dreiBein.ScalarProduct(newPosition) < 0.) ? -1. : +1.;
+        LOG(6, "DEBUG: oldCenter on plane is " << oldPosition
+            << ", newCenter in plane is " << newPosition
+            << ", and dreiBein is " << dreiBein);
+        const double RotationAngle = sign * oldPosition.Angle(newPosition);
+        LOG(5, "DEBUG: Rotating around self is " << RotationAngle
+            << " around axis " << RotationAxis);
+
 #ifndef NDEBUG
+        // check: first bestmatching in rotated_newpolygon and remainingnew
+        // should now equal
         {
-        // check that rotation works in DEBUG mode
-        boost::math::quaternion<double> p(0., newCenter[0], newCenter[1], newCenter[2]);
-        boost::math::quaternion<double> q_inverse =
-            boost::math::conj(q)/(boost::math::norm(q));
-        p = q * p * q_inverse;
-        boost::math::quaternion<double> identity(1,0,0,0);
-        ASSERT( boost::math::norm(q*q_inverse - identity) < std::numeric_limits<double>::epsilon()*1e4,
-            "matchSphericalPointDistributions() - quaternion does not rotate newCenter "
-            +toString(q*q_inverse)+" into oldCenter "+toString(identity)+".");
-        boost::math::quaternion<double> comparison(0., -oldCenter[0], oldCenter[1], oldCenter[2]);
-        ASSERT( boost::math::norm(p - comparison) < std::numeric_limits<double>::epsilon()*1e4,
-            "matchSphericalPointDistributions() - quaternion does not rotate newCenter "
-            +toString(p)+" into oldCenter "+toString(comparison)+".");
+          const IndexList_t::const_iterator iter = MCS.bestmatching.begin();
+          Vector rotatednew = RotationAxis.rotateVector(
+              rotated_newpolygon[*iter],
+              RotationAngle);
+          LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
+              << " while old was " << remainingold[*iter]);
+          ASSERT( (rotatednew - remainingold[*iter]).Norm()
+              < std::numeric_limits<double>::epsilon()*1e4,
+              "matchSphericalPointDistributions() - orientation rotation does not work as expected.");
         }
 #endif
 
-        // rotate spherical distribution
-        rotated_newpolygon = rotatePolygon(_newpolygon, q);
-        LOG(5, "DEBUG: Rotated new polygon is " << rotated_newpolygon);
-      } else {
-        rotated_newpolygon = _newpolygon;
-      }
-    }
-    // rotate triangle/line/point around itself to match orientation
-    const Line RotationAxis(zeroVec, oldCenter);
-    const double RotationAngle =
-        oldCenter.Angle(remainingold[0])
-        - oldCenter.Angle(remainingnew[*MCS.bestmatching.begin()]);
-    LOG(5, "DEBUG: Rotate around self is " << RotationAngle
-        << " around axis " << RotationAxis);
-
-    for (SphericalPointDistribution::Polygon_t::iterator iter = rotated_newpolygon.begin();
-        iter != rotated_newpolygon.end(); ++iter) {
-      RotationAxis.rotateVector(*iter, RotationAngle);
+        for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
+            iter != rotated_newpolygon.end(); ++iter) {
+          Vector &current = *iter;
+          LOG(6, "DEBUG: Original point is " << current);
+          current = RotationAxis.rotateVector(current, RotationAngle);
+          LOG(6, "DEBUG: Rotated point is " << current);
+        }
+      }
     }