Changeset 1ae9aa for src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Timestamp:

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

Author:

Frederik Heber <heber@…>

Children:

23c605

Parents:

260540

git-author:

Frederik Heber <heber@…> (06/29/14 21:20:49)

git-committer:

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

Message:

Extracted joinPoints() function to make it accessible to tests.

• added unit test.
• moved some function definitions around and added documentation.

File:

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

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -62 +62 @@
 typedef std::vector<double> DistanceArray_t;
 
+// class generator: taken from www.cplusplus.com example std::generate
+struct c_unique {
+  int current;
+  c_unique() {current=0;}
+  int operator()() {return current++;}
+} UniqueNumber;
+
 inline
 DistanceArray_t calculatePairwiseDistances(
@@ -82 +89 @@
 }
 
-// class generator: taken from www.cplusplus.com example std::generate
-struct c_unique {
-  int current;
-  c_unique() {current=0;}
-  int operator()() {return current++;}
-} UniqueNumber;
-
+/** 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;
+}
+
+/** Decides by an orthonormal third vector whether the sign of the rotation
+ * angle should be negative or positive.
+ *
+ * \return -1 or 1
+ */
+inline
+double determineSignOfRotation(
+    const Vector &_oldPosition,
+    const Vector &_newPosition,
+    const Vector &_RotationAxis
+    )
+{
+  Vector dreiBein(_oldPosition);
+  dreiBein.VectorProduct(_RotationAxis);
+  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);
+  return sign;
+}
+
+/** Convenience function to recalculate old and new center for determining plane
+ * rotation.
+ */
+inline
+void calculateOldAndNewCenters(
+    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);
+  // 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);
+}
 /** Returns squared L2 error of the given \a _Matching.
  *
@@ -218 +283 @@
 }
 
-/** Decides by an orthonormal third vector whether the sign of the rotation
- * angle should be negative or positive.
- *
- * \return -1 or 1
- */
-inline
-double determineSignOfRotation(
-    const Vector &_oldPosition,
-    const Vector &_newPosition,
-    const Vector &_RotationAxis
-    )
-{
-  Vector dreiBein(_oldPosition);
-  dreiBein.VectorProduct(_RotationAxis);
-  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);
-  return sign;
-}
-
 /** Finds combinatorially the best matching between points in \a _polygon
  * and \a _newpolygon.
@@ -246 +288 @@
  * We find the matching with the smallest L2 error, where we break when we stumble
  * upon a matching with zero error.
+ *
+ * As points in \a _polygon may be have a weight greater 1 we have to match it to
+ * multiple points in \a _newpolygon. Eventually, these multiple points are combined
+ * for their center of weight, which is the only thing follow-up function see of
+ * these multiple points. Hence, we actually modify \a _newpolygon in the process
+ * such that the returned IndexList_t indicates a bijective mapping in the end.
  *
  * \sa recurseMatchings() for going through all matchings
@@ -254 +302 @@
  */
 SphericalPointDistribution::IndexList_t SphericalPointDistribution::findBestMatching(
-    const SphericalPointDistribution::WeightedPolygon_t &_polygon,
-    const SphericalPointDistribution::Polygon_t &_newpolygon
+    const WeightedPolygon_t &_polygon,
+    Polygon_t &_newpolygon
     )
 {
@@ -279 +327 @@
     recurseMatchings(MCS, matching, indices, matchingsize);
   }
-  return MCS.bestmatching;
-}
-
-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;
-}
-
-inline
-void calculateOldAndNewCenters(
-    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);
-  // 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);
+
+  // combine multiple points and create simple IndexList from IndexTupleList
+  IndexTupleList_t bestmatching;
+  for (IndexList_t::const_iterator iter = MCS.bestmatching.begin();
+      iter != MCS.bestmatching.end(); ++iter)
+    bestmatching.push_back(IndexList_t(1, *iter));
+  const SphericalPointDistribution::IndexList_t IndexList =
+      joinPoints(_newpolygon, MCS.newpoints, bestmatching);
+
+  return IndexList;
+}
+
+SphericalPointDistribution::IndexList_t SphericalPointDistribution::joinPoints(
+    Polygon_t &_newpolygon,
+    const VectorArray_t &_newpoints,
+    const IndexTupleList_t &_bestmatching
+    )
+{
+  // combine all multiple points
+  IndexList_t IndexList;
+  IndexArray_t removalpoints;
+  unsigned int UniqueIndex = _newpolygon.size(); // all indices up to size are used right now
+  VectorArray_t newCenters;
+  newCenters.reserve(_bestmatching.size());
+  for (IndexTupleList_t::const_iterator tupleiter = _bestmatching.begin();
+      tupleiter != _bestmatching.end(); ++tupleiter) {
+    ASSERT (tupleiter->size() > 0,
+        "findBestMatching() - encountered tuple in bestmatching with size 0.");
+    if (tupleiter->size() == 1) {
+      // add point and index
+      IndexList.push_back(*tupleiter->begin());
+    } else {
+      // combine into weighted and normalized center
+      Vector Center = calculateCenter(_newpoints, *tupleiter);
+      Center.Normalize();
+      _newpolygon.push_back(Center);
+      LOG(5, "DEBUG: Combining " << tupleiter->size() << "points to weighted center "
+          << Center << " with new index " << UniqueIndex);
+      // mark for removal
+      removalpoints.insert(removalpoints.end(), tupleiter->begin(), tupleiter->end());
+      // add new index
+      IndexList.push_back(UniqueIndex++);
+    }
+  }
+  // IndexList is now our new bestmatching (that is bijective)
+  LOG(4, "DEBUG: Our new bijective IndexList reads as " << IndexList);
+
+  // modifying _newpolygon: remove all points in removalpoints, add those in newCenters
+  Polygon_t allnewpoints = _newpolygon;
+  {
+    _newpolygon.clear();
+    std::sort(removalpoints.begin(), removalpoints.end());
+    size_t i = 0;
+    IndexArray_t::const_iterator removeiter = removalpoints.begin();
+    for (Polygon_t::iterator iter = allnewpoints.begin();
+        iter != allnewpoints.end(); ++iter, ++i) {
+      if ((removeiter != removalpoints.end()) && (i == *removeiter)) {
+        // don't add, go to next remove index
+        ++removeiter;
+      } else {
+        // otherwise add points
+        _newpolygon.push_back(*iter);
+      }
+    }
+  }
+  LOG(4, "DEBUG: The polygon with recentered points removed is " << _newpolygon);
+
+  // map IndexList to new shrinked _newpolygon
+  typedef std::set<unsigned int> IndexSet_t;
+  IndexSet_t SortedIndexList(IndexList.begin(), IndexList.end());
+  IndexList.clear();
+  {
+    size_t offset = 0;
+    IndexSet_t::const_iterator listiter = SortedIndexList.begin();
+    IndexArray_t::const_iterator removeiter = removalpoints.begin();
+    for (size_t i = 0; i < allnewpoints.size(); ++i) {
+      if ((removeiter != removalpoints.end()) && (i == *removeiter)) {
+        ++offset;
+        ++removeiter;
+      } else if ((listiter != SortedIndexList.end()) && (i == *listiter)) {
+        IndexList.push_back(*listiter - offset);
+        ++listiter;
+      }
+    }
+  }
+  LOG(4, "DEBUG: Our new bijective IndexList corrected for removed points reads as "
+      << IndexList);
+
+  return IndexList;
 }
 
@@ -492 +592 @@
 SphericalPointDistribution::matchSphericalPointDistributions(
     const SphericalPointDistribution::WeightedPolygon_t &_polygon,
-    const SphericalPointDistribution::Polygon_t &_newpolygon
+    SphericalPointDistribution::Polygon_t &_newpolygon
     )
 {