Changeset 0d5ca7 for src/Fragmentation/Exporters

Timestamp:

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

Author:

Frederik Heber <heber@…>

Children:

e139180

Parents:

90426a

git-author:

Frederik Heber <heber@…> (06/02/14 14:51:33)

git-committer:

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

Message:

FIX: SphericalPointDistribution succeeds with unit test.

• removed XFAIL from FragmentationAutomation regression tests.

Location:

src/Fragmentation/Exporters

Files:

• SphericalPointDistribution.cpp (modified) (11 diffs)
• unittests/Makefile.am (modified) (1 diff)

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -44 +44 @@
 #include <algorithm>
 #include <cmath>
+#include <functional>
+#include <iterator>
 #include <limits>
 #include <list>
@@ -81 +83 @@
   int current;
   c_unique() {current=0;}
-  int operator()() {return ++current;}
+  int operator()() {return current++;}
 } UniqueNumber;
 
@@ -125 +127 @@
     }
   } else
-    ELOG(2, "calculateErrorOfMatching() - Given matching is empty.");
+    ELOG(3, "calculateErrorOfMatching() - Given matching's size is less than 2.");
   LOG(3, "INFO: Resulting errors for matching (L1, L2): "
       << errors.first << "," << errors.second << ".");
@@ -156 +158 @@
 }
 
-/** Rotates a given polygon around x, y, and z axis by the given angles.
- *
- * Essentially, we concentrate on the three points of the polygon to rotate
- * to the correct position. First, we rotate its center via \a angles,
- * then we rotate the "triangle" around itself/\a _RotationAxis by
- * \a _RotationAngle.
- *
- * \param _polygon polygon whose points to rotate
- * \param _angles vector with rotation angles for x,y,z axis
- * \param _RotationAxis
- * \param _RotationAngle
- */
-SphericalPointDistribution::Polygon_t rotatePolygon(
-    const SphericalPointDistribution::Polygon_t &_polygon,
-    const std::vector<double> &_angles,
-    const Line &_RotationAxis,
-    const double _RotationAngle)
-{
-  SphericalPointDistribution::Polygon_t rotated_polygon = _polygon;
-  RealSpaceMatrix rotation;
-  ASSERT( _angles.size() == 3,
-      "rotatePolygon() - not exactly "+toString(3)+" angles given.");
-  rotation.setRotation(_angles[0] * M_PI/180., _angles[1] * M_PI/180., _angles[2] * M_PI/180.);
-  LOG(4, "DEBUG: Rotation matrix is " << rotation);
-
-  // apply rotation angles
-  for (SphericalPointDistribution::Polygon_t::iterator iter = rotated_polygon.begin();
-      iter != rotated_polygon.end(); ++iter) {
-    *iter = rotation * (*iter);
-    _RotationAxis.rotateVector(*iter, _RotationAngle);
-  }
-
-  return rotated_polygon;
-}
-
 struct MatchingControlStructure {
   bool foundflag;
@@ -214 +181 @@
   LOG(4, "DEBUG: Recursing with current matching " << _matching
       << ", remaining indices " << _indices
-      << ", and sought size " << _matchingsize);
+      << ", and sought size " << _matching.size()+_matchingsize);
   //!> threshold for L1 error below which matching is immediately acceptable
   const double L1THRESHOLD = 1e-2;
   if (!_MCS.foundflag) {
-    LOG(3, "INFO: Current matching has size " << _matching.size() << " of " << _matchingsize);
-    if (_matching.size() < _matchingsize) {
+    LOG(4, "DEBUG: Current matching has size " << _matching.size() << ", places left " << _matchingsize);
+    if (_matchingsize > 0) {
       // go through all indices
       for (IndexList_t::iterator iter = _indices.begin();
-          iter != _indices.end();) {
+          (iter != _indices.end()) && (!_MCS.foundflag);) {
         // add index to matching
         _matching.push_back(*iter);
-        LOG(4, "DEBUG: Adding " << *iter << " to matching.");
+        LOG(5, "DEBUG: Adding " << *iter << " to matching.");
         // remove index but keep iterator to position (is the next to erase element)
         IndexList_t::iterator backupiter = _indices.erase(iter);
@@ -237 +204 @@
       }
       // gone through all indices then exit recursion
-      _MCS.foundflag = true;
+      if (_matching.empty())
+        _MCS.foundflag = true;
     } else {
       LOG(3, "INFO: Found matching " << _matching);
@@ -246 +214 @@
         _MCS.bestmatching = _matching;
         _MCS.foundflag = true;
-      }
-      if (_MCS.bestL2 > errors.second) {
+      } else if (_MCS.bestL2 > errors.second) {
         _MCS.bestmatching = _matching;
         _MCS.bestL2 = errors.second;
@@ -254 +221 @@
   }
 }
+
+/** Convert cartesian to polar coordinates.
+ *
+ * \param _cartesian vector in cartesian coordinates
+ * \return vector containing \f$ (r,\theta, \varphi \f$ tuple for polar coordinates
+ */
+std::vector<double> getPolarCoordinates(const Vector &_cartesian)
+{
+  std::vector<double> polar(3,0.);
+  const double xsqr = _cartesian[0] * _cartesian[0];
+  const double ysqr = _cartesian[1] * _cartesian[1];
+  polar[0] = sqrt(xsqr + ysqr + _cartesian[2]*_cartesian[2]);
+  if (fabs(_cartesian[2]) < std::numeric_limits<double>::epsilon()*1e4) {
+    if (fabs(xsqr + ysqr) < std::numeric_limits<double>::epsilon()*1e4) {
+      polar[1] = 0.;
+    } else {
+      // xsqr + ysqr is always non-negative
+      polar[1] = M_PI/2.;
+    }
+  } else {
+    polar[1] = atan( sqrt(xsqr + ysqr)/_cartesian[2]);
+    if (_cartesian[2] <= -std::numeric_limits<double>::epsilon()*1e4)
+      polar[1] += M_PI;
+  }
+
+  if (fabs(_cartesian[0]) < std::numeric_limits<double>::epsilon()*1e4) {
+    if (fabs(_cartesian[1]) < std::numeric_limits<double>::epsilon()*1e4) {
+      polar[2] = 0.;
+    } else if (_cartesian[1] > std::numeric_limits<double>::epsilon()*1e4) {
+      polar[2] = M_PI/2.;
+    } else {
+      polar[2] = -M_PI/2.;
+    }
+  } else {
+    polar[2] = atan ( _cartesian[1]/_cartesian[0] );
+    if (_cartesian[0] <= -std::numeric_limits<double>::epsilon()*1e4)
+      polar[2] += M_PI;
+  }
+  return polar;
+}
+
+/** Calculate cartesian coordinates from given polar ones.
+ *
+ * \param _polar vector with polar coordinates
+ * \return cartesian coordinates
+ */
+Vector getCartesianCoordinates(const std::vector<double> &_polar)
+{
+  Vector cartesian;
+  ASSERT( _polar.size() == 3,
+      "convertToCartesianCoordinates() - tuples has insufficient components.");
+  cartesian[0] = _polar[0] * sin(_polar[1]) * cos(_polar[2]);
+  cartesian[1] = _polar[0] * sin(_polar[1]) * sin(_polar[2]);
+  cartesian[2] = _polar[0] * cos(_polar[1]);
+  return cartesian;
+}
+
+/** Rotates a given polygon around x, y, and z axis by the given angles.
+ *
+ * \param _polygon polygon whose points to rotate
+ * \param _angles vector with rotation angles for x,y,z axis
+ */
+SphericalPointDistribution::Polygon_t rotatePolygon(
+    const SphericalPointDistribution::Polygon_t &_polygon,
+    const std::vector<double> &_angles)
+{
+  SphericalPointDistribution::Polygon_t rotated_polygon = _polygon;
+  RealSpaceMatrix rotation;
+  ASSERT( _angles.size() == 3,
+      "rotatePolygon() - not exactly "+toString(3)+" components given.");
+
+  // apply rotation angles
+  for (SphericalPointDistribution::Polygon_t::iterator iter = rotated_polygon.begin();
+      iter != rotated_polygon.end(); ++iter) {
+    // transform to polar
+    std::vector<double> polar = getPolarCoordinates(*iter);
+    LOG(5, "DEBUG: Converting " << *iter << " to " << polar);
+    // sum up difference
+    std::transform(
+        _angles.begin(), _angles.end(),
+        polar.begin(),
+        polar.begin(),
+        std::plus<double>());
+    // convert back
+    *iter = getCartesianCoordinates(polar);
+    LOG(5, "DEBUG: Converting modified " << polar << " to " << *iter);
+  }
+
+  return rotated_polygon;
+}
+
 
 SphericalPointDistribution::Polygon_t
@@ -264 +322 @@
   VectorArray_t remainingold(_polygon.begin(), _polygon.end());
   VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
-  LOG(3, "INFO: Matching old polygon " << _polygon
+  LOG(2, "INFO: Matching old polygon " << _polygon
       << " with new polygon " << _newpolygon);
 
@@ -287 +345 @@
       recurseMatchings(MCS, matching, indices, matchingsize);
     }
-    LOG(3, "INFO: Best matching is " << MCS.bestmatching);
+    LOG(2, "INFO: Best matching is " << MCS.bestmatching);
 
     // determine rotation angles to align the two point distributions with
     // respect to bestmatching
-    std::vector<double> angles(3);
+    std::vector<double> angles(NDIM);
     Vector newCenter;
+    Vector oldCenter;
     {
       // calculate center of triangle/line/point consisting of first points of matching
-      Vector oldCenter;
       IndexList_t::const_iterator iter = MCS.bestmatching.begin();
       unsigned int i = 0;
@@ -304 +362 @@
       oldCenter *= 1./(double)i;
       newCenter *= 1./(double)i;
-      LOG(3, "INFO: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
-
-      Vector direction(0.,0.,0.);
-      for(size_t i=0;i<NDIM;++i) {
-        // create new rotation axis
-        direction[i] = 1.;
-        const Line axis (zeroVec, direction);
-        // calculate rotation angle for this axis
-        const double alpha = direction.Angle(oldCenter) - direction.Angle(newCenter);
-        // perform rotation
-        axis.rotateVector(newCenter, alpha);
-        // store angle
-        angles[i] = alpha;
-        // reset direction component for next iteration
-        direction[i] = 0.;
-      }
-    }
-    LOG(3, "INFO: (x,y,z) angles are" << angles);
-    const Line RotationAxis(zeroVec, newCenter);
-    const double RotationAngle =
-        newCenter.Angle(remainingold[0])
-        - newCenter.Angle(remainingnew[*MCS.bestmatching.begin()]);
-    LOG(3, "INFO: Rotate around self is " << RotationAngle
-        << " around axis " << RotationAxis);
-
+      LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
+
+      // transform to polar coordinates and note difference in angular parts
+      std::vector<double> oldpolar = getPolarCoordinates(oldCenter);
+      std::vector<double> newpolar = getPolarCoordinates(newCenter);
+      std::vector<double> differencepolar;
+      std::transform(
+          oldpolar.begin(), oldpolar.end(),
+          newpolar.begin(),
+          std::back_inserter(differencepolar),
+          std::minus<double>());
+      LOG(3, "INFO: (r,theta,phi) angles are" << differencepolar);
+    }
     // rotate _newpolygon
     SphericalPointDistribution::Polygon_t rotated_newpolygon =
-        rotatePolygon(_newpolygon, angles, RotationAxis, RotationAngle);
-    LOG(3, "INFO: Rotated new polygon is " << rotated_newpolygon);
+        rotatePolygon(_newpolygon, angles);
+    LOG(5, "DEBUG: Rotated new polygon is " << rotated_newpolygon);
+
+    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);
+    }
 
     // remove all points in matching and return remaining ones
-    return removeMatchingPoints(rotated_newpolygon, MCS.bestmatching);
+    SphericalPointDistribution::Polygon_t remainingpoints =
+        removeMatchingPoints(rotated_newpolygon, MCS.bestmatching);
+    LOG(2, "INFO: Remaining points are " << remainingpoints);
+    return remainingpoints;
   } else
     return _newpolygon;

src/Fragmentation/Exporters/unittests/Makefile.am

@@ -18 +18 @@
 
 TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
-XFAIL_TESTS += SphericalPointDistributionUnitTest
 check_PROGRAMS += $(FRAGMENTATIONEXPORTERSTESTS) 
 noinst_PROGRAMS += $(FRAGMENTATIONEXPORTERSTESTS)