Changeset 4b96da for src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Timestamp:

May 25, 2016, 7:13:59 AM (9 years ago)

Author:

Frederik Heber <heber@…>

Children:

66700f2

Parents:

4d611d

git-author:

Frederik Heber <heber@…> (07/09/14 21:14:53)

git-committer:

Frederik Heber <heber@…> (05/25/16 07:13:59)

Message:

Refactoring SphericalPointDistribution to combine point list and subset of indices.

• added struct PolygonWithIndices.

File:

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

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -99 +99 @@
 }
 
+/** Calculate the center of a given set of points in \a _positions but only
+ * for those indicated by \a _indices.
+ *
+ */
+inline
+Vector calculateCenter(
+    const SphericalPointDistribution::PolygonWithIndices &_polygon)
+{
+  return calculateCenter(_polygon.polygon, _polygon.indices);
+}
+
 inline
 DistanceArray_t calculatePairwiseDistances(
@@ -177 +188 @@
     Vector &_oldCenter,
     Vector &_newCenter,
-    const SphericalPointDistribution::VectorArray_t &_referencepositions,
-    const SphericalPointDistribution::VectorArray_t &_currentpositions,
-    const SphericalPointDistribution::IndexList_t &_bestmatching)
-{
-  const size_t NumberIds = std::min(_bestmatching.size(), (size_t)3);
-  SphericalPointDistribution::IndexList_t continuousIds(NumberIds, -1);
-  std::generate(continuousIds.begin(), continuousIds.end(), UniqueNumber);
-  _oldCenter = calculateCenter(_referencepositions, continuousIds);
+    const SphericalPointDistribution::PolygonWithIndices &_referencepositions,
+    const SphericalPointDistribution::PolygonWithIndices &_currentpositions)
+{
+  _oldCenter = calculateCenter(_referencepositions.polygon, _referencepositions.indices);
   // C++11 defines a copy_n function ...
-  SphericalPointDistribution::IndexList_t::const_iterator enditer = _bestmatching.begin();
-  std::advance(enditer, NumberIds);
-  SphericalPointDistribution::IndexList_t firstbestmatchingIds(NumberIds, -1);
-  std::copy(_bestmatching.begin(), enditer, firstbestmatchingIds.begin());
-  _newCenter = calculateCenter( _currentpositions, firstbestmatchingIds);
+  _newCenter = calculateCenter( _currentpositions.polygon, _currentpositions.indices);
 }
 /** Returns squared L2 error of the given \a _Matching.
@@ -251 +254 @@
 
 SphericalPointDistribution::Polygon_t SphericalPointDistribution::removeMatchingPoints(
-    const VectorArray_t &_points,
-    const IndexList_t &_matchingindices
+    const PolygonWithIndices &_points
     )
 {
-  SphericalPointDistribution::Polygon_t remainingreturnpolygon;
-  IndexArray_t indices(_matchingindices.begin(), _matchingindices.end());
+  SphericalPointDistribution::Polygon_t remainingpoints;
+  IndexArray_t indices(_points.indices.begin(), _points.indices.end());
   std::sort(indices.begin(), indices.end());
   LOG(4, "DEBUG: sorted matching is " << indices);
-  IndexArray_t remainingindices(_points.size(), -1);
+  IndexArray_t remainingindices(_points.polygon.size(), -1);
   std::generate(remainingindices.begin(), remainingindices.end(), UniqueNumber);
   IndexArray_t::iterator remainiter = std::set_difference(
@@ -269 +271 @@
   for (IndexArray_t::const_iterator iter = remainingindices.begin();
       iter != remainingindices.end(); ++iter) {
-          remainingreturnpolygon.push_back(_points[*iter]);
-  }
-
-  return remainingreturnpolygon;
+    remainingpoints.push_back(_points.polygon[*iter]);
+  }
+
+  return remainingpoints;
 }
 
@@ -501 +503 @@
 
 SphericalPointDistribution::Rotation_t SphericalPointDistribution::findPlaneAligningRotation(
-    const VectorArray_t &_referencepositions,
-    const VectorArray_t &_currentpositions,
-    const IndexList_t &_bestmatching
+    const PolygonWithIndices &_referencepositions,
+    const PolygonWithIndices &_currentpositions
     )
 {
@@ -518 +519 @@
   calculateOldAndNewCenters(
       oldCenter, newCenter,
-      _referencepositions, _currentpositions, _bestmatching);
+      _referencepositions, _currentpositions);
 
   if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
@@ -536 +537 @@
     if ((oldCenter.IsZero()) && (newCenter.IsZero())) {
       // either oldCenter or newCenter (or both) is directly at origin
-      if (_bestmatching.size() == 2) {
+      if (_currentpositions.indices.size() == 2) {
         // line case
-        Vector oldPosition = _currentpositions[*_bestmatching.begin()];
-        Vector newPosition = _referencepositions[0];
+        Vector oldPosition = _currentpositions.polygon[*_currentpositions.indices.begin()];
+        Vector newPosition = _referencepositions.polygon[0];
         // check whether we need to rotate at all
         if (!oldPosition.IsEqualTo(newPosition)) {
@@ -553 +554 @@
         // both triangles/planes have same center, hence get axis by
         // VectorProduct of Normals
-        Plane newplane(_referencepositions[0], _referencepositions[1], _referencepositions[2]);
+        Plane newplane(_referencepositions.polygon[0], _referencepositions.polygon[1], _referencepositions.polygon[2]);
         VectorArray_t vectors;
-        for (IndexList_t::const_iterator iter = _bestmatching.begin();
-            iter != _bestmatching.end(); ++iter)
-          vectors.push_back(_currentpositions[*iter]);
+        for (IndexList_t::const_iterator iter = _currentpositions.indices.begin();
+            iter != _currentpositions.indices.end(); ++iter)
+          vectors.push_back(_currentpositions.polygon[*iter]);
         Plane oldplane(vectors[0], vectors[1], vectors[2]);
         Vector oldPosition = oldplane.getNormal();
@@ -603 +604 @@
 
 SphericalPointDistribution::Rotation_t SphericalPointDistribution::findPointAligningRotation(
-    const VectorArray_t &remainingold,
-    const VectorArray_t &remainingnew,
-    const IndexList_t &_bestmatching)
+    const PolygonWithIndices &remainingold,
+    const PolygonWithIndices &remainingnew)
 {
   // initialize rotation to zero
@@ -618 +618 @@
   calculateOldAndNewCenters(
       oldCenter, newCenter,
-      remainingold, remainingnew, _bestmatching);
-
-  Vector oldPosition = remainingnew[*_bestmatching.begin()];
-  Vector newPosition = remainingold[0];
-  LOG(6, "DEBUG: oldPosition is " << oldPosition << " and newPosition is " << newPosition);
+      remainingold, remainingnew);
+
+  Vector oldPosition = remainingnew.polygon[*remainingnew.indices.begin()];
+  Vector newPosition = remainingold.polygon[0];
+  LOG(6, "DEBUG: oldPosition is " << oldPosition << " (length: "
+      << oldPosition.Norm() << ") and newPosition is " << newPosition << " length(: "
+      << newPosition.Norm() << ")");
   if (!oldPosition.IsEqualTo(newPosition)) {
     if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
-      oldCenter.Normalize();  // note weighted sum of normalized weight is not normalized
+      // we rotate at oldCenter and around the radial direction, which is again given
+      // by oldCenter.
       Rotation.first = oldCenter;
-      LOG(6, "DEBUG: Picking normalized oldCenter as Rotation.first " << oldCenter);
-      oldPosition.ProjectOntoPlane(Rotation.first);
-      newPosition.ProjectOntoPlane(Rotation.first);
+      Rotation.first.Normalize();  // note weighted sum of normalized weight is not normalized
+      LOG(6, "DEBUG: Using oldCenter " << oldCenter << " as rotation center and "
+          << Rotation.first << " as axis.");
+      oldPosition -= oldCenter;
+      newPosition -= oldCenter;
+      oldPosition = (oldPosition - oldPosition.Projection(Rotation.first));
+      newPosition = (newPosition - newPosition.Projection(Rotation.first));
       LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
     } else {
-      if (_bestmatching.size() == 2) {
+      if (remainingnew.indices.size() == 2) {
         // line situation
         try {
@@ -653 +660 @@
       } else {
         // triangle situation
-        Plane oldplane(remainingold[0], remainingold[1], remainingold[2]);
+        Plane oldplane(remainingold.polygon[0], remainingold.polygon[1], remainingold.polygon[2]);
         Rotation.first = oldplane.getNormal();
         LOG(6, "DEBUG: oldPlane is " << oldplane << " and normal is " << Rotation.first);
@@ -678 +685 @@
 {
   SphericalPointDistribution::Polygon_t remainingpoints;
-  VectorArray_t remainingold;
-  for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
-      iter != _polygon.end(); ++iter)
-    remainingold.push_back(iter->first);
+
   LOG(2, "INFO: Matching old polygon " << _polygon
       << " with new polygon " << _newpolygon);
@@ -695 +699 @@
     LOG(2, "INFO: Best matching is " << bestmatching);
 
-    // _newpolygon has changed, so now convert to array
-    VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
+    const size_t NumberIds = std::min(bestmatching.size(), (size_t)3);
+    // create old set
+    PolygonWithIndices oldSet;
+    oldSet.indices.resize(NumberIds, -1);
+    std::generate(oldSet.indices.begin(), oldSet.indices.end(), UniqueNumber);
+    for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
+        iter != _polygon.end(); ++iter)
+      oldSet.polygon.push_back(iter->first);
+
+    // _newpolygon has changed, so now convert to array with matched indices
+    PolygonWithIndices newSet;
+    SphericalPointDistribution::IndexList_t::const_iterator beginiter = bestmatching.begin();
+    SphericalPointDistribution::IndexList_t::const_iterator enditer = bestmatching.begin();
+    std::advance(enditer, NumberIds);
+    newSet.indices.resize(NumberIds, -1);
+    std::copy(beginiter, enditer, newSet.indices.begin());
+    std::copy(_newpolygon.begin(),_newpolygon.end(), std::back_inserter(newSet.polygon));
 
     // determine rotation angles to align the two point distributions with
@@ -702 +721 @@
     // we use the center between the three first matching points
     /// the first rotation brings these two centers to coincide
-    VectorArray_t rotated_newpolygon = remainingnew;
+    PolygonWithIndices rotatednewSet = newSet;
     {
-      Rotation_t Rotation = findPlaneAligningRotation(
-          remainingold,
-          remainingnew,
-          bestmatching);
+      Rotation_t Rotation = findPlaneAligningRotation(oldSet, newSet);
       LOG(5, "DEBUG: Rotating coordinate system by " << Rotation.second
           << " around axis " << Rotation.first);
@@ -713 +729 @@
 
       // apply rotation angle to bring newCenter to oldCenter
-      for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
-          iter != rotated_newpolygon.end(); ++iter) {
+      for (VectorArray_t::iterator iter = rotatednewSet.polygon.begin();
+          iter != rotatednewSet.polygon.end(); ++iter) {
         Vector &current = *iter;
         LOG(6, "DEBUG: Original point is " << current);
@@ -728 +744 @@
         calculateOldAndNewCenters(
             oldCenter, rotatednewCenter,
-            remainingold, rotated_newpolygon, bestmatching);
+            oldSet, rotatednewSet);
         // NOTE: Center must not necessarily lie on the sphere with norm 1, hence, we
         // have to normalize it just as before, as oldCenter and newCenter lengths may differ.
@@ -736 +752 @@
           LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
               << ", oldCenter is " << oldCenter);
-          ASSERT( (rotatednewCenter - oldCenter).NormSquared() < std::numeric_limits<double>::epsilon()*1e4,
+          const double difference = (rotatednewCenter - oldCenter).NormSquared();
+          ASSERT( difference < std::numeric_limits<double>::epsilon()*1e4,
               "matchSphericalPointDistributions() - rotation does not work as expected by "
-              +toString((rotatednewCenter - oldCenter).NormSquared())+".");
+              +toString(difference)+".");
         }
       }
@@ -746 +763 @@
     /// points themselves coincide
     if (bestmatching.size() > 1) {
-      Rotation_t Rotation = findPointAligningRotation(
-          remainingold,
-          rotated_newpolygon,
-          bestmatching);
+      Rotation_t Rotation = findPointAligningRotation(oldSet, rotatednewSet);
 
       // construct RotationAxis and two points on its plane, defining the angle
@@ -763 +777 @@
       {
         const IndexList_t::const_iterator iter = bestmatching.begin();
+
+        // check whether both old and newPosition are at same distance to oldCenter
+        Vector oldCenter = calculateCenter(oldSet);
+        const double distance = fabs(
+              (oldSet.polygon[0] - oldCenter).NormSquared()
+              - (rotatednewSet.polygon[*iter] - oldCenter).NormSquared()
+            );
+        LOG(4, "CHECK: Squared distance between oldPosition and newPosition "
+            << " with respect to oldCenter " << oldCenter << " is " << distance);
+//        ASSERT( distance < warn_amplitude,
+//            "matchSphericalPointDistributions() - old and newPosition's squared distance to oldCenter differs by "
+//            +toString(distance));
+
         Vector rotatednew = RotationAxis.rotateVector(
-            rotated_newpolygon[*iter],
+            rotatednewSet.polygon[*iter],
             Rotation.second);
         LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
-            << " while old was " << remainingold[0]);
-        ASSERT( (rotatednew - remainingold[0]).Norm() < warn_amplitude,
-            "matchSphericalPointDistributions() - orientation rotation ends up off by more than "
-            +toString(warn_amplitude)+".");
+            << " while old was " << oldSet.polygon[0]);
+        const double difference = (rotatednew - oldSet.polygon[0]).NormSquared();
+        ASSERT( difference < distance+1e-8,
+            "matchSphericalPointDistributions() - orientation rotation ends up off by "
+            +toString(difference)+", more than "
+            +toString(distance+1e-8)+".");
       }
 #endif
 
-      for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
-          iter != rotated_newpolygon.end(); ++iter) {
+      for (VectorArray_t::iterator iter = rotatednewSet.polygon.begin();
+          iter != rotatednewSet.polygon.end(); ++iter) {
         Vector &current = *iter;
         LOG(6, "DEBUG: Original point is " << current);
@@ -785 +814 @@
     // remove all points in matching and return remaining ones
     SphericalPointDistribution::Polygon_t remainingpoints =
-        removeMatchingPoints(rotated_newpolygon, bestmatching);
+        removeMatchingPoints(rotatednewSet);
     LOG(2, "INFO: Remaining points are " << remainingpoints);
     return remainingpoints;