Changeset 39616b for src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Timestamp:

Aug 20, 2014, 1:07:11 PM (11 years ago)

Author:

Frederik Heber <heber@…>

Children:

7e45c4

Parents:

e94448

git-author:

Frederik Heber <heber@…> (07/21/14 08:17:42)

git-committer:

Frederik Heber <heber@…> (08/20/14 13:07:11)

Message:

findPointAlignment() calculates mean rotation angle.

• added calculateCenterOfGravity(), as we should use the center of gravity and not the center of sphere as rotational center (i.e. point on the plane to be rotated).

File:

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

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -321 +321 @@
  */
 inline
+Vector calculateCenterOfGravity(
+    const SphericalPointDistribution::VectorArray_t &_positions,
+    const SphericalPointDistribution::IndexList_t &_indices)
+{
+  Vector Center;
+  Center.Zero();
+  for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
+      iter != _indices.end(); ++iter)
+    Center += _positions[*iter];
+  if (!_indices.empty())
+    Center *= 1./(double)_indices.size();
+
+  return Center;
+}
+
+/** 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::VectorArray_t &_positions,
@@ -838 +858 @@
       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)) {
-    // we rotate at oldCenter and around the radial direction, which is again given
-    // by oldCenter.
-    Rotation.first = oldCenter;
-    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);
-
-    // construct reference vector to determine direction of rotation
-    const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
-    Rotation.second = sign * oldPosition.Angle(newPosition);
-  } else {
-    LOG(6, "DEBUG: oldPosition and newPosition are equivalent, hence no orientating rotation.");
-  }
+  // we rotate at oldCenter and around the radial direction, which is again given
+  // by oldCenter.
+  Rotation.first = oldCenter;
+  Rotation.first.Normalize();  // note weighted sum of normalized weight is not normalized
+
+  // calculate center of the rotational plane
+  newCenter = calculateCenterOfGravity(remainingnew.polygon, remainingnew.indices);
+  oldCenter = calculateCenterOfGravity(remainingold.polygon, remainingold.indices);
+  LOG(6, "DEBUG: Using oldCenter " << oldCenter << " as rotation center and "
+      << Rotation.first << " as axis.");
+
+  LOG(6, "DEBUG: old indices are " << remainingold.indices
+      << ", new indices are " << remainingnew.indices);
+  IndexList_t::const_iterator newiter = remainingnew.indices.begin();
+  IndexList_t::const_iterator olditer = remainingold.indices.begin();
+  for (;
+      (newiter != remainingnew.indices.end()) && (olditer != remainingold.indices.end());
+      ++newiter,++olditer) {
+    Vector newPosition = remainingnew.polygon[*newiter];
+    Vector oldPosition = remainingold.polygon[*olditer];
+    LOG(6, "DEBUG: oldPosition is " << oldPosition << " (length: "
+        << oldPosition.Norm() << ") and newPosition is " << newPosition << " length(: "
+        << newPosition.Norm() << ")");
+
+    if (!oldPosition.IsEqualTo(newPosition)) {
+      oldPosition -= oldCenter;
+      newPosition -= newCenter;
+      oldPosition = (oldPosition - oldPosition.Projection(Rotation.first));
+      newPosition = (newPosition - newPosition.Projection(Rotation.first));
+      LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
+
+      // construct reference vector to determine direction of rotation
+      const double sign = determineSignOfRotation(newPosition, oldPosition, Rotation.first);
+      LOG(6, "DEBUG: Determining angle between " << oldPosition << " and " << newPosition);
+      const double angle = sign * newPosition.Angle(oldPosition);
+      LOG(6, "DEBUG: Adding " << angle << " to weighted rotation sum.");
+      Rotation.second += angle;
+    } else {
+      LOG(6, "DEBUG: oldPosition and newPosition are equivalent, hence no orientating rotation.");
+    }
+  }
+  Rotation.second *= 1./(double)remainingnew.indices.size();
 
   return Rotation;
@@ -1062 +1099 @@
             << " while old was " << oldSet.polygon[0]);
         const double difference = (rotatednew - oldSet.polygon[0]).NormSquared();
-        ASSERT( difference < distance+1e-8,
+        ASSERT( difference < distance+warn_amplitude,
             "matchSphericalPointDistributions() - orientation rotation ends up off by "
             +toString(difference)+", more than "
-            +toString(distance+1e-8)+".");
+            +toString(distance+warn_amplitude)+".");
       }
 #endif
@@ -1232 +1269 @@
             << " while old was " << oldSet.polygon[0]);
         const double difference = (rotatednew - oldSet.polygon[0]).NormSquared();
-        ASSERT( difference < distance+1e-8,
+        ASSERT( difference < distance+warn_amplitude,
             "matchSphericalPointDistributions() - orientation rotation ends up off by "
             +toString(difference)+", more than "
-            +toString(distance+1e-8)+".");
+            +toString(distance+warn_amplitude)+".");
       }
 #endif