Changeset a2f8a9 for src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Timestamp:

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

Author:

Frederik Heber <heber@…>

Children:

ef3885

Parents:

1cde4e8

git-author:

Frederik Heber <heber@…> (07/09/14 22:08:37)

git-committer:

Frederik Heber <heber@…> (08/20/14 13:06:16)

Message:

Added calculation of center of minimum distance by bisection.

• this will give us a unique a definite point independent of the (rotational) position of the point set on the unit sphere.
• added unit test.

File:

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

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -43 +43 @@
 
 #include <algorithm>
-#include <boost/math/quaternion.hpp>
+#include <boost/assign.hpp>
 #include <cmath>
 #include <functional>
@@ -58 +58 @@
 #include "LinearAlgebra/Vector.hpp"
 
+using namespace boost::assign;
+
 // static entities
 const double SphericalPointDistribution::warn_amplitude = 1e-2;
@@ -83 +85 @@
  */
 inline
-Vector calculateCenter(
+Vector calculateGeographicMidpoint(
     const SphericalPointDistribution::VectorArray_t &_positions,
     const SphericalPointDistribution::IndexList_t &_indices)
@@ -96 +98 @@
 
   return Center;
+}
+
+inline
+double calculateMinimumDistance(
+    const Vector &_center,
+    const SphericalPointDistribution::VectorArray_t &_points,
+    const SphericalPointDistribution::IndexList_t & _indices)
+{
+  double MinimumDistance = 0.;
+  for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
+      iter != _indices.end(); ++iter) {
+    const double angle = _center.Angle(_points[*iter]);
+    MinimumDistance += angle*angle;
+  }
+  return sqrt(MinimumDistance);
+}
+
+/** Calculates the center of minimum distance for a given set of points \a _points.
+ *
+ * According to http://www.geomidpoint.com/calculation.html this goes a follows:
+ * -# Let CurrentPoint be the geographic midpoint found in Method A. this is used as the starting point for the search.
+ * -# Let MinimumDistance be the sum total of all distances from the current point to all locations in 'Your Places'.
+ * -# Find the total distance between each location in 'Your Places' and all other locations in 'Your Places'. If any one of these locations has a new smallest distance then that location becomes the new CurrentPoint and MinimumDistance.
+ * -# Let TestDistance be PI/2 radians (6225 miles or 10018 km).
+ * -# Find the total distance between each of 8 test points and all locations in 'Your Places'. The test points are arranged in a circular pattern around the CurrentPoint at a distance of TestDistance to the north, northeast, east, southeast, south, southwest, west and northwest.
+ * -# If any of these 8 points has a new smallest distance then that point becomes the new CurrentPoint and MinimumDistance and go back to step 5 using that point.
+ * -# If none of the 8 test points has a new smallest distance then divide TestDistance by 2 and go back to step 5 using the same point.
+ * -# Repeat steps 5 to 7 until no new smallest distance can be found or until TestDistance is less than 0.00000002 radians.
+ *
+ * \param _points set of points
+ * \return Center of minimum distance for all these points, is always of length 1
+ */
+Vector SphericalPointDistribution::calculateCenterOfMinimumDistance(
+    const SphericalPointDistribution::VectorArray_t &_positions,
+    const SphericalPointDistribution::IndexList_t &_indices)
+{
+  ASSERT( _positions.size() >= _indices.size(),
+      "calculateCenterOfMinimumDistance() - less positions than indices given.");
+  Vector center(1.,0.,0.);
+
+  /// first treat some special cases
+  // no positions given: return x axis vector (NOT zero!)
+  {
+    if (_indices.empty())
+      return center;
+    // one position given: return it directly
+    if (_positions.size() == (size_t)1)
+      return _positions[0];
+    // two positions on a line given: return closest point to (1.,0.,0.)
+    if (fabs(_positions[0].ScalarProduct(_positions[1]) + 1.)
+        < std::numeric_limits<double>::epsilon()*1e4) {
+      Vector candidate;
+      if (_positions[0].ScalarProduct(center) > _positions[1].ScalarProduct(center))
+        candidate = _positions[0];
+      else
+        candidate = _positions[1];
+      // non-uniqueness: all positions on great circle, normal to given line are valid
+      // so, we just pick one because returning a unique point is topmost priority
+      Vector normal;
+      normal.GetOneNormalVector(candidate);
+      Vector othernormal = candidate;
+      othernormal.VectorProduct(normal);
+      // now both normal and othernormal describe the plane containing the great circle
+      Plane greatcircle(normal, zeroVec, othernormal);
+      // check which axis is contained and pick the one not
+      if (greatcircle.isContained(center)) {
+        center = Vector(0.,1.,0.);
+        if (greatcircle.isContained(center))
+          center = Vector(0.,0.,1.);
+        // now we are done cause a plane cannot contain all three axis vectors
+      }
+      center = greatcircle.getClosestPoint(candidate);
+      // assure length of 1
+      center.Normalize();
+    }
+  }
+
+  // start with geographic midpoint
+  center = calculateGeographicMidpoint(_positions, _indices);
+  if (!center.IsZero()) {
+    center.Normalize();
+    LOG(4, "DEBUG: Starting with geographical midpoint of " << _positions << " under indices "
+        << _indices << " is " << center);
+  } else {
+    // any point is good actually
+    center = _positions[0];
+    LOG(4, "DEBUG: Starting with first position " << center);
+  }
+
+  // calculate initial MinimumDistance
+  double MinimumDistance = calculateMinimumDistance(center, _positions, _indices);
+  LOG(5, "DEBUG: MinimumDistance to this center is " << MinimumDistance);
+
+  // check all present points whether they may serve as a better center
+  for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
+      iter != _indices.end(); ++iter) {
+    const Vector &centerCandidate = _positions[*iter];
+    const double candidateDistance = calculateMinimumDistance(centerCandidate, _positions, _indices);
+    if (candidateDistance < MinimumDistance) {
+      MinimumDistance = candidateDistance;
+      center = centerCandidate;
+      LOG(5, "DEBUG: new MinimumDistance to current test point " << MinimumDistance
+          << " is " << center);
+    }
+  }
+  LOG(5, "DEBUG: new MinimumDistance to center " << center << " is " << MinimumDistance);
+
+  // now iterate over TestDistance
+  double TestDistance = Vector(1.,0.,0.).Angle(Vector(0.,1.,0.));
+//  LOG(6, "DEBUG: initial TestDistance is " << TestDistance);
+
+  const double threshold = sqrt(std::numeric_limits<double>::epsilon());
+  // check each of eight test points at N, NE, E, SE, S, SW, W, NW
+  typedef std::vector<double> angles_t;
+  angles_t testangles;
+  testangles += 0./180.*M_PI, 45./180.*M_PI, 90./180.*M_PI, 135./180.*M_PI,
+      180./180.*M_PI, 225./180.*M_PI, 270./180.*M_PI, 315./180.*M_PI;
+  while (TestDistance > threshold) {
+    Vector OneNormal;
+    OneNormal.GetOneNormalVector(center);
+    Line RotationAxis(zeroVec, OneNormal);
+    Vector North = RotationAxis.rotateVector(center,TestDistance);
+    Line CompassRose(zeroVec, center);
+    bool updatedflag = false;
+    for (angles_t::const_iterator angleiter = testangles.begin();
+        angleiter != testangles.end(); ++angleiter) {
+      Vector centerCandidate = CompassRose.rotateVector(North, *angleiter);
+//      centerCandidate.Normalize();
+      const double candidateDistance = calculateMinimumDistance(centerCandidate, _positions, _indices);
+      if (candidateDistance < MinimumDistance) {
+        MinimumDistance = candidateDistance;
+        center = centerCandidate;
+        updatedflag = true;
+        LOG(5, "DEBUG: new MinimumDistance to test point at " << *angleiter/M_PI*180.
+            << "° is " << MinimumDistance);
+      }
+    }
+
+    // if no new point, decrease TestDistance
+    if (!updatedflag) {
+      TestDistance *= 0.5;
+//      LOG(6, "DEBUG: TestDistance is now " << TestDistance);
+    }
+  }
+  LOG(4, "DEBUG: Final MinimumDistance to center " << center << " is " << MinimumDistance);
+
+  return center;
+}
+
+Vector calculateCenterOfMinimumDistance(
+    const SphericalPointDistribution::PolygonWithIndices &_points)
+{
+  return SphericalPointDistribution::calculateCenterOfMinimumDistance(_points.polygon, _points.indices);
+}
+
+/** 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,
+    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;
+  return SphericalPointDistribution::calculateCenterOfMinimumDistance(_positions, _indices);
 }
 
@@ -744 +920 @@
             oldCenter, rotatednewCenter,
             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.
-        if ((!oldCenter.IsZero()) && (!rotatednewCenter.IsZero())) {
-          oldCenter.Normalize();
-          rotatednewCenter.Normalize();
-          LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
-              << ", oldCenter is " << oldCenter);
-          const double difference = (rotatednewCenter - oldCenter).NormSquared();
-          ASSERT( difference < std::numeric_limits<double>::epsilon()*1e4,
-              "matchSphericalPointDistributions() - rotation does not work as expected by "
-              +toString(difference)+".");
-        }
+        oldCenter.Normalize();
+        rotatednewCenter.Normalize();
+        // check whether centers are anti-parallel (factor -1)
+        // then we have the "non-unique poles" situation: points lie on great circle
+        // and both poles are valid solution
+        if (fabs(oldCenter.ScalarProduct(rotatednewCenter) + 1.)
+            < std::numeric_limits<double>::epsilon()*1e4)
+          rotatednewCenter *= -1.;
+        LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
+            << ", oldCenter is " << oldCenter);
+        const double difference = (rotatednewCenter - oldCenter).NormSquared();
+        ASSERT( difference < std::numeric_limits<double>::epsilon()*1e4,
+            "matchSphericalPointDistributions() - rotation does not work as expected by "
+            +toString(difference)+".");
       }
 #endif