Ignore:
Timestamp:
Sep 10, 2016, 4:14:02 PM (9 years ago)
Author:
Frederik Heber <heber@…>
Branches:
SaturateAtoms_singleDegree
Children:
899cc9
Parents:
226860
git-author:
Frederik Heber <heber@…> (06/12/14 07:23:12)
git-committer:
Frederik Heber <heber@…> (09/10/16 16:14:02)
Message:

Using the idea of three points giving a triangle to find rotation axis.

  • we calculate the center of either triangle and rotate the center of the ideal point distribution to match the one from the given points.
  • next we have the triangles normals as axis, take the first matching point and rotate align it.
  • we have to deal with a lot of special cases: What if only zero, one, or two points are given ...
  • in general we assume that the triangle lies relatively flat on the sphere's surface but what if the origin is in the triangle plane or even the calculated center is at the origin ...
  • TESTS: SphericalPointDistributionUnitTest working again, regression tests FragmentMolecule-cylces and StoreSaturatedFragment working.
Location:
src/Fragmentation/Exporters
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/SphericalPointDistribution.cpp

    r226860 rbcf13b  
    5353
    5454#include "LinearAlgebra/Line.hpp"
     55#include "LinearAlgebra/Plane.hpp"
    5556#include "LinearAlgebra/RealSpaceMatrix.hpp"
    5657#include "LinearAlgebra/Vector.hpp"
    5758
    58 typedef std::list<unsigned int> IndexList_t;
    59 typedef std::vector<unsigned int> IndexArray_t;
    60 typedef std::vector<Vector> VectorArray_t;
     59// static entities
     60const double SphericalPointDistribution::SQRT_3(sqrt(3.0));
     61const double SphericalPointDistribution::warn_amplitude = 1e-2;
     62
    6163typedef std::vector<double> DistanceArray_t;
    6264
    63 // static instances
    64 const double SphericalPointDistribution::SQRT_3(sqrt(3.0));
    65 
     65inline
    6666DistanceArray_t calculatePairwiseDistances(
    67     const std::vector<Vector> &_returnpolygon,
    68     const IndexList_t &_indices
     67    const std::vector<Vector> &_points,
     68    const SphericalPointDistribution::IndexList_t &_indices
    6969    )
    7070{
    7171  DistanceArray_t result;
    72   for (IndexList_t::const_iterator firstiter = _indices.begin();
     72  for (SphericalPointDistribution::IndexList_t::const_iterator firstiter = _indices.begin();
    7373      firstiter != _indices.end(); ++firstiter) {
    74     for (IndexList_t::const_iterator seconditer = firstiter;
     74    for (SphericalPointDistribution::IndexList_t::const_iterator seconditer = firstiter;
    7575        seconditer != _indices.end(); ++seconditer) {
    7676      if (firstiter == seconditer)
    7777        continue;
    78       const double distance = (_returnpolygon[*firstiter] - _returnpolygon[*seconditer]).NormSquared();
     78      const double distance = (_points[*firstiter] - _points[*seconditer]).NormSquared();
    7979      result.push_back(distance);
    8080    }
    … …  
    101101 * \return pair with L1 and squared L2 error
    102102 */
    103 std::pair<double, double> calculateErrorOfMatching(
     103std::pair<double, double> SphericalPointDistribution::calculateErrorOfMatching(
    104104    const std::vector<Vector> &_old,
    105105    const std::vector<Vector> &_new,
    … …  
    138138}
    139139
    140 SphericalPointDistribution::Polygon_t removeMatchingPoints(
     140SphericalPointDistribution::Polygon_t SphericalPointDistribution::removeMatchingPoints(
    141141    const VectorArray_t &_points,
    142142    const IndexList_t &_matchingindices
    … …  
    163163}
    164164
    165 struct MatchingControlStructure {
    166   bool foundflag;
    167   double bestL2;
    168   IndexList_t bestmatching;
    169   VectorArray_t oldreturnpolygon;
    170   VectorArray_t newreturnpolygon;
    171 };
    172 
    173165/** Recursive function to go through all possible matchings.
    174166 *
    … …  
    178170 * \param _matchingsize
    179171 */
    180 void recurseMatchings(
     172void SphericalPointDistribution::recurseMatchings(
    181173    MatchingControlStructure &_MCS,
    182174    IndexList_t &_matching,
    … …  
    215207      // calculate errors
    216208      std::pair<double, double> errors = calculateErrorOfMatching(
    217           _MCS.oldreturnpolygon, _MCS.newreturnpolygon, _matching);
     209          _MCS.oldpoints, _MCS.newpoints, _matching);
    218210      if (errors.first < L1THRESHOLD) {
    219211        _MCS.bestmatching = _matching;
    … …  
    227219}
    228220
    229 /** Rotates a given polygon around x, y, and z axis by the given angles.
    230  *
    231  * \param _polygon polygon whose points to rotate
    232  * \param _q quaternion specifying the rotation of the coordinate system
     221/** Decides by an orthonormal third vector whether the sign of the rotation
     222 * angle should be negative or positive.
     223 *
     224 * \return -1 or 1
    233225 */
    234 SphericalPointDistribution::Polygon_t rotatePolygon(
     226inline
     227double determineSignOfRotation(
     228    const Vector &_oldPosition,
     229    const Vector &_newPosition,
     230    const Vector &_RotationAxis
     231    )
     232{
     233  Vector dreiBein(_oldPosition);
     234  dreiBein.VectorProduct(_RotationAxis);
     235  dreiBein.Normalize();
     236  const double sign =
     237      (dreiBein.ScalarProduct(_newPosition) < 0.) ? -1. : +1.;
     238  LOG(6, "DEBUG: oldCenter on plane is " << _oldPosition
     239      << ", newCenter in plane is " << _newPosition
     240      << ", and dreiBein is " << dreiBein);
     241  return sign;
     242}
     243
     244/** Finds combinatorially the best matching between points in \a _polygon
     245 * and \a _newpolygon.
     246 *
     247 * We find the matching with the smallest L2 error, where we break when we stumble
     248 * upon a matching with zero error.
     249 *
     250 * \sa recurseMatchings() for going through all matchings
     251 *
     252 * \param _polygon here, we have indices 0,1,2,...
     253 * \param _newpolygon and here we need to find the correct indices
     254 * \return list of indices: first in \a _polygon goes to first index for \a _newpolygon
     255 */
     256SphericalPointDistribution::IndexList_t SphericalPointDistribution::findBestMatching(
    235257    const SphericalPointDistribution::Polygon_t &_polygon,
    236     const boost::math::quaternion<double> &_q)
    237 {
    238   SphericalPointDistribution::Polygon_t rotated_polygon = _polygon;
    239   boost::math::quaternion<double> q_inverse =
    240       boost::math::conj(_q)/(boost::math::norm(_q));
    241 
    242   // apply rotation angles
    243   for (SphericalPointDistribution::Polygon_t::iterator iter = rotated_polygon.begin();
    244       iter != rotated_polygon.end(); ++iter) {
    245     Vector &current = *iter;
    246     boost::math::quaternion<double> p(0, current[0], current[1], current[2]);
    247     p = _q * p * q_inverse;
    248     LOG(5, "DEBUG: Rotated point is " << p);
    249     // i have no idea why but first component comes up with wrong sign
    250     current[0] = -p.R_component_2();
    251     current[1] = p.R_component_3();
    252     current[2] = p.R_component_4();
    253   }
    254 
    255   return rotated_polygon;
     258    const SphericalPointDistribution::Polygon_t &_newpolygon
     259    )
     260{
     261  MatchingControlStructure MCS;
     262  MCS.foundflag = false;
     263  MCS.bestL2 = std::numeric_limits<double>::max();
     264  MCS.oldpoints.insert(MCS.oldpoints.begin(), _polygon.begin(),_polygon.end() );
     265  MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );
     266
     267  // search for bestmatching combinatorially
     268  {
     269    // translate polygon into vector to enable index addressing
     270    IndexList_t indices(_newpolygon.size());
     271    std::generate(indices.begin(), indices.end(), UniqueNumber);
     272    IndexList_t matching;
     273
     274    // walk through all matchings
     275    const unsigned int matchingsize = _polygon.size();
     276    ASSERT( matchingsize <= indices.size(),
     277        "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new points to choose for matching to old ones.");
     278    recurseMatchings(MCS, matching, indices, matchingsize);
     279  }
     280  return MCS.bestmatching;
     281}
     282
     283inline
     284Vector calculateCenter(
     285    const SphericalPointDistribution::VectorArray_t &_positions,
     286    const SphericalPointDistribution::IndexList_t &_indices)
     287{
     288  Vector Center;
     289  Center.Zero();
     290  for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
     291      iter != _indices.end(); ++iter)
     292    Center += _positions[*iter];
     293  if (!_indices.empty())
     294    Center *= 1./(double)_indices.size();
     295
     296  return Center;
     297}
     298
     299inline
     300void calculateOldAndNewCenters(
     301    Vector &_oldCenter,
     302    Vector &_newCenter,
     303    const SphericalPointDistribution::VectorArray_t &_referencepositions,
     304    const SphericalPointDistribution::VectorArray_t &_currentpositions,
     305    const SphericalPointDistribution::IndexList_t &_bestmatching)
     306{
     307  const size_t NumberIds = std::min(_bestmatching.size(), (size_t)3);
     308  SphericalPointDistribution::IndexList_t continuousIds(NumberIds, -1);
     309  std::generate(continuousIds.begin(), continuousIds.end(), UniqueNumber);
     310  _oldCenter = calculateCenter(_referencepositions, continuousIds);
     311  // C++11 defines a copy_n function ...
     312  SphericalPointDistribution::IndexList_t::const_iterator enditer = _bestmatching.begin();
     313  std::advance(enditer, NumberIds);
     314  SphericalPointDistribution::IndexList_t firstbestmatchingIds(NumberIds, -1);
     315  std::copy(_bestmatching.begin(), enditer, firstbestmatchingIds.begin());
     316  _newCenter = calculateCenter( _currentpositions, firstbestmatchingIds);
     317}
     318
     319SphericalPointDistribution::Rotation_t SphericalPointDistribution::findPlaneAligningRotation(
     320    const VectorArray_t &_referencepositions,
     321    const VectorArray_t &_currentpositions,
     322    const IndexList_t &_bestmatching
     323    )
     324{
     325  bool dontcheck = false;
     326  // initialize to no rotation
     327  Rotation_t Rotation;
     328  Rotation.first.Zero();
     329  Rotation.first[0] = 1.;
     330  Rotation.second = 0.;
     331
     332  // calculate center of triangle/line/point consisting of first points of matching
     333  Vector oldCenter;
     334  Vector newCenter;
     335  calculateOldAndNewCenters(
     336      oldCenter, newCenter,
     337      _referencepositions, _currentpositions, _bestmatching);
     338
     339  if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
     340    LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
     341    oldCenter.Normalize();
     342    newCenter.Normalize();
     343    if (!oldCenter.IsEqualTo(newCenter)) {
     344      // calculate rotation axis and angle
     345      Rotation.first = oldCenter;
     346      Rotation.first.VectorProduct(newCenter);
     347      Rotation.second = oldCenter.Angle(newCenter); // /(M_PI/2.);
     348    } else {
     349      // no rotation required anymore
     350    }
     351  } else {
     352    LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
     353    if ((oldCenter.IsZero()) && (newCenter.IsZero())) {
     354      // either oldCenter or newCenter (or both) is directly at origin
     355      if (_bestmatching.size() == 2) {
     356        // line case
     357        Vector oldPosition = _currentpositions[*_bestmatching.begin()];
     358        Vector newPosition = _referencepositions[0];
     359        // check whether we need to rotate at all
     360        if (!oldPosition.IsEqualTo(newPosition)) {
     361          Rotation.first = oldPosition;
     362          Rotation.first.VectorProduct(newPosition);
     363          // orientation will fix the sign here eventually
     364          Rotation.second = oldPosition.Angle(newPosition);
     365        } else {
     366          // no rotation required anymore
     367        }
     368      } else {
     369        // triangle case
     370        // both triangles/planes have same center, hence get axis by
     371        // VectorProduct of Normals
     372        Plane newplane(_referencepositions[0], _referencepositions[1], _referencepositions[2]);
     373        VectorArray_t vectors;
     374        for (IndexList_t::const_iterator iter = _bestmatching.begin();
     375            iter != _bestmatching.end(); ++iter)
     376          vectors.push_back(_currentpositions[*iter]);
     377        Plane oldplane(vectors[0], vectors[1], vectors[2]);
     378        Vector oldPosition = oldplane.getNormal();
     379        Vector newPosition = newplane.getNormal();
     380        // check whether we need to rotate at all
     381        if (!oldPosition.IsEqualTo(newPosition)) {
     382          Rotation.first = oldPosition;
     383          Rotation.first.VectorProduct(newPosition);
     384          Rotation.first.Normalize();
     385
     386          // construct reference vector to determine direction of rotation
     387          const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
     388          Rotation.second = sign * oldPosition.Angle(newPosition);
     389          LOG(5, "DEBUG: Rotating plane normals by " << Rotation.second
     390              << " around axis " << Rotation.first);
     391        } else {
     392          // else do nothing
     393        }
     394      }
     395    } else {
     396      // TODO: we can't do anything here, but this case needs to be dealt with when
     397      // we have no ideal geometries anymore
     398      if ((oldCenter-newCenter).Norm() > warn_amplitude)
     399        ELOG(2, "oldCenter is " << oldCenter << ", yet newCenter is " << newCenter);
     400      // else they are considered close enough
     401      dontcheck = true;
     402    }
     403  }
     404
     405#ifndef NDEBUG
     406  // check: rotation brings newCenter onto oldCenter position
     407  if (!dontcheck) {
     408    Line Axis(zeroVec, Rotation.first);
     409    Vector test = Axis.rotateVector(newCenter, Rotation.second);
     410    LOG(4, "CHECK: rotated newCenter is " << test
     411        << ", oldCenter is " << oldCenter);
     412    ASSERT( (test - oldCenter).NormSquared() < std::numeric_limits<double>::epsilon()*1e4,
     413        "matchSphericalPointDistributions() - rotation does not work as expected by "
     414        +toString((test - oldCenter).NormSquared())+".");
     415  }
     416#endif
     417
     418  return Rotation;
     419}
     420
     421SphericalPointDistribution::Rotation_t SphericalPointDistribution::findPointAligningRotation(
     422    const VectorArray_t &remainingold,
     423    const VectorArray_t &remainingnew,
     424    const IndexList_t &_bestmatching)
     425{
     426  // initialize rotation to zero
     427  Rotation_t Rotation;
     428  Rotation.first.Zero();
     429  Rotation.first[0] = 1.;
     430  Rotation.second = 0.;
     431
     432  // recalculate center
     433  Vector oldCenter;
     434  Vector newCenter;
     435  calculateOldAndNewCenters(
     436      oldCenter, newCenter,
     437      remainingold, remainingnew, _bestmatching);
     438
     439  Vector oldPosition = remainingnew[*_bestmatching.begin()];
     440  Vector newPosition = remainingold[0];
     441  LOG(6, "DEBUG: oldPosition is " << oldPosition << " and newPosition is " << newPosition);
     442  if (!oldPosition.IsEqualTo(newPosition)) {
     443    if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
     444      oldCenter.Normalize();  // note weighted sum of normalized weight is not normalized
     445      Rotation.first = oldCenter;
     446      LOG(6, "DEBUG: Picking normalized oldCenter as Rotation.first " << oldCenter);
     447      oldPosition.ProjectOntoPlane(Rotation.first);
     448      newPosition.ProjectOntoPlane(Rotation.first);
     449      LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
     450    } else {
     451      if (_bestmatching.size() == 2) {
     452        // line situation
     453        try {
     454          Plane oldplane(oldPosition, oldCenter, newPosition);
     455          Rotation.first = oldplane.getNormal();
     456          LOG(6, "DEBUG: Plane is " << oldplane << " and normal is " << Rotation.first);
     457        } catch (LinearDependenceException &e) {
     458          LOG(6, "DEBUG: Vectors defining plane are linearly dependent.");
     459          // oldPosition and newPosition are on a line, just flip when not equal
     460          if (!oldPosition.IsEqualTo(newPosition)) {
     461            Rotation.first.Zero();
     462            Rotation.first.GetOneNormalVector(oldPosition);
     463            LOG(6, "DEBUG: For flipping we use Rotation.first " << Rotation.first);
     464            assert( Rotation.first.ScalarProduct(oldPosition) < std::numeric_limits<double>::epsilon()*1e4);
     465  //              Rotation.second = M_PI;
     466          } else {
     467            LOG(6, "DEBUG: oldPosition and newPosition are equivalent.");
     468          }
     469        }
     470      } else {
     471        // triangle situation
     472        Plane oldplane(remainingold[0], remainingold[1], remainingold[2]);
     473        Rotation.first = oldplane.getNormal();
     474        LOG(6, "DEBUG: oldPlane is " << oldplane << " and normal is " << Rotation.first);
     475        oldPosition.ProjectOntoPlane(Rotation.first);
     476        LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
     477      }
     478    }
     479    // construct reference vector to determine direction of rotation
     480    const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
     481    Rotation.second = sign * oldPosition.Angle(newPosition);
     482  } else {
     483    LOG(6, "DEBUG: oldPosition and newPosition are equivalent, hence no orientating rotation.");
     484  }
     485
     486  return Rotation;
    256487}
    257488
    … …  
    269500      << " with new polygon " << _newpolygon);
    270501
     502  if (_polygon.size() == _newpolygon.size()) {
     503    // same number of points desired as are present? Do nothing
     504    LOG(2, "INFO: There are no vacant points to return.");
     505    return remainingreturnpolygon;
     506  }
     507
    271508  if (_polygon.size() > 0) {
    272     MatchingControlStructure MCS;
    273     MCS.foundflag = false;
    274     MCS.bestL2 = std::numeric_limits<double>::max();
    275     MCS.oldreturnpolygon.insert(MCS.oldreturnpolygon.begin(), _polygon.begin(),_polygon.end() );
    276     MCS.newreturnpolygon.insert(MCS.newreturnpolygon.begin(), _newpolygon.begin(),_newpolygon.end() );
    277 
    278     // search for bestmatching combinatorially
     509    IndexList_t bestmatching = findBestMatching(_polygon, _newpolygon);
     510    LOG(2, "INFO: Best matching is " << bestmatching);
     511
     512    // determine rotation angles to align the two point distributions with
     513    // respect to bestmatching:
     514    // we use the center between the three first matching points
     515    /// the first rotation brings these two centers to coincide
     516    VectorArray_t rotated_newpolygon = remainingnew;
    279517    {
    280       // translate polygon into vector to enable index addressing
    281       IndexList_t indices(_newpolygon.size());
    282       std::generate(indices.begin(), indices.end(), UniqueNumber);
    283       IndexList_t matching;
    284 
    285       // walk through all matchings
    286       const unsigned int matchingsize = _polygon.size();
    287       ASSERT( matchingsize <= indices.size(),
    288           "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new returnpolygon to choose for matching to old ones.");
    289       recurseMatchings(MCS, matching, indices, matchingsize);
    290     }
    291     LOG(2, "INFO: Best matching is " << MCS.bestmatching);
    292 
    293     // determine rotation angles to align the two point distributions with
    294     // respect to bestmatching
    295     VectorArray_t rotated_newpolygon = remainingnew;
    296     Vector oldCenter;
    297     {
    298       // calculate center of triangle/line/point consisting of first points of matching
    299       Vector newCenter;
    300       IndexList_t::const_iterator iter = MCS.bestmatching.begin();
    301       unsigned int i = 0;
    302       for (; (i<3) && (i<MCS.bestmatching.size()); ++i, ++iter) {
    303         oldCenter += remainingold[i];
    304         newCenter += remainingnew[*iter];
    305       }
    306       oldCenter *= 1./(double)i;
    307       newCenter *= 1./(double)i;
    308       LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
    309 
    310       if ((oldCenter - newCenter).NormSquared() > std::numeric_limits<double>::epsilon()*1e4) {
    311         // setup quaternion
    312         Vector RotationAxis = oldCenter;
    313         RotationAxis.VectorProduct(newCenter);
    314         Line Axis(zeroVec, RotationAxis);
    315         RotationAxis.Normalize();
    316         const double RotationAngle = oldCenter.Angle(newCenter); // /(M_PI/2.);
    317         LOG(5, "DEBUG: Rotate coordinate system by " << RotationAngle
    318             << " around axis " << RotationAxis);
    319 
    320         // apply rotation angles
    321         for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
    322             iter != rotated_newpolygon.end(); ++iter) {
    323           Vector &current = *iter;
    324           LOG(5, "DEBUG: Original point is " << current);
    325           current =  Axis.rotateVector(current, RotationAngle);
    326           LOG(5, "DEBUG: Rotated point is " << current);
     518      Rotation_t Rotation = findPlaneAligningRotation(
     519          remainingold,
     520          remainingnew,
     521          bestmatching);
     522      LOG(5, "DEBUG: Rotating coordinate system by " << Rotation.second
     523          << " around axis " << Rotation.first);
     524      Line Axis(zeroVec, Rotation.first);
     525
     526      // apply rotation angle to bring newCenter to oldCenter
     527      for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
     528          iter != rotated_newpolygon.end(); ++iter) {
     529        Vector &current = *iter;
     530        LOG(6, "DEBUG: Original point is " << current);
     531        current =  Axis.rotateVector(current, Rotation.second);
     532        LOG(6, "DEBUG: Rotated point is " << current);
     533      }
     534
     535#ifndef NDEBUG
     536      // check: rotated "newCenter" should now equal oldCenter
     537      {
     538        Vector oldCenter;
     539        Vector rotatednewCenter;
     540        calculateOldAndNewCenters(
     541            oldCenter, rotatednewCenter,
     542            remainingold, rotated_newpolygon, bestmatching);
     543        // NOTE: Center must not necessarily lie on the sphere with norm 1, hence, we
     544        // have to normalize it just as before, as oldCenter and newCenter lengths may differ.
     545        if ((!oldCenter.IsZero()) && (!rotatednewCenter.IsZero())) {
     546          oldCenter.Normalize();
     547          rotatednewCenter.Normalize();
     548          LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
     549              << ", oldCenter is " << oldCenter);
     550          ASSERT( (rotatednewCenter - oldCenter).NormSquared() < std::numeric_limits<double>::epsilon()*1e4,
     551              "matchSphericalPointDistributions() - rotation does not work as expected by "
     552              +toString((rotatednewCenter - oldCenter).NormSquared())+".");
    327553        }
    328554      }
    329     }
    330     // rotate triangle/line/point around itself to match orientation
    331     if (MCS.bestmatching.size() > 1) {
    332       if (oldCenter.NormSquared() > std::numeric_limits<double>::epsilon()*1e4) {
    333         // construct RotationAxis and two points on its plane, defining the angle
    334         const Line RotationAxis(zeroVec, oldCenter);
    335         Vector oldPosition(rotated_newpolygon[*MCS.bestmatching.begin()]);
    336         oldPosition.ProjectOntoPlane(RotationAxis.getDirection());
    337         Vector newPosition(remainingold[*MCS.bestmatching.begin()]);
    338         newPosition.ProjectOntoPlane(RotationAxis.getDirection());
    339 
    340         // construct reference vector to determine direction of rotation
    341         Vector dreiBein(oldPosition);
    342         dreiBein.VectorProduct(oldCenter);
    343         dreiBein.Normalize();
    344         const double sign =
    345             (dreiBein.ScalarProduct(newPosition) < 0.) ? -1. : +1.;
    346         LOG(6, "DEBUG: oldCenter on plane is " << oldPosition
    347             << ", newCenter in plane is " << newPosition
    348             << ", and dreiBein is " << dreiBein);
    349         const double RotationAngle = sign * oldPosition.Angle(newPosition);
    350         LOG(5, "DEBUG: Rotating around self is " << RotationAngle
    351             << " around axis " << RotationAxis);
     555#endif
     556    }
     557    /// the second (orientation) rotation aligns the planes such that the
     558    /// points themselves coincide
     559    if (bestmatching.size() > 1) {
     560      Rotation_t Rotation = findPointAligningRotation(
     561          remainingold,
     562          rotated_newpolygon,
     563          bestmatching);
     564
     565      // construct RotationAxis and two points on its plane, defining the angle
     566      Rotation.first.Normalize();
     567      const Line RotationAxis(zeroVec, Rotation.first);
     568
     569      LOG(5, "DEBUG: Rotating around self is " << Rotation.second
     570          << " around axis " << RotationAxis);
    352571
    353572#ifndef NDEBUG
    354         // check: first bestmatching in rotated_newpolygon and remainingnew
    355         // should now equal
    356         {
    357           const IndexList_t::const_iterator iter = MCS.bestmatching.begin();
    358           Vector rotatednew = RotationAxis.rotateVector(
    359               rotated_newpolygon[*iter],
    360               RotationAngle);
    361           LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
    362               << " while old was " << remainingold[*iter]);
    363           ASSERT( (rotatednew - remainingold[*iter]).Norm()
    364               < std::numeric_limits<double>::epsilon()*1e4,
    365               "matchSphericalPointDistributions() - orientation rotation does not work as expected.");
    366         }
     573      // check: first bestmatching in rotated_newpolygon and remainingnew
     574      // should now equal
     575      {
     576        const IndexList_t::const_iterator iter = bestmatching.begin();
     577        Vector rotatednew = RotationAxis.rotateVector(
     578            rotated_newpolygon[*iter],
     579            Rotation.second);
     580        LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
     581            << " while old was " << remainingold[0]);
     582        ASSERT( (rotatednew - remainingold[0]).Norm() < warn_amplitude,
     583            "matchSphericalPointDistributions() - orientation rotation ends up off by more than "
     584            +toString(warn_amplitude)+".");
     585      }
    367586#endif
    368587
    369         for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
    370             iter != rotated_newpolygon.end(); ++iter) {
    371           Vector &current = *iter;
    372           LOG(6, "DEBUG: Original point is " << current);
    373           current = RotationAxis.rotateVector(current, RotationAngle);
    374           LOG(6, "DEBUG: Rotated point is " << current);
    375         }
     588      for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
     589          iter != rotated_newpolygon.end(); ++iter) {
     590        Vector &current = *iter;
     591        LOG(6, "DEBUG: Original point is " << current);
     592        current = RotationAxis.rotateVector(current, Rotation.second);
     593        LOG(6, "DEBUG: Rotated point is " << current);
    376594      }
    377595    }
    … …  
    379597    // remove all points in matching and return remaining ones
    380598    SphericalPointDistribution::Polygon_t remainingpoints =
    381         removeMatchingPoints(rotated_newpolygon, MCS.bestmatching);
     599        removeMatchingPoints(rotated_newpolygon, bestmatching);
    382600    LOG(2, "INFO: Remaining points are " << remainingpoints);
    383601    return remainingpoints;

  • src/Fragmentation/Exporters/SphericalPointDistribution.hpp

    r226860 rbcf13b  
    7676  //!> precalculated value for root of 3
    7777  static const double SQRT_3;
     78
     79  typedef std::pair<Vector, double> Rotation_t;
     80
     81  typedef std::list<unsigned int> IndexList_t;
     82  typedef std::vector<unsigned int> IndexArray_t;
     83  typedef std::vector<Vector> VectorArray_t;
     84
     85  //!> amplitude up to which deviations in checks of rotations are tolerated
     86  static const double warn_amplitude;
     87
     88private:
     89  static std::pair<double, double> calculateErrorOfMatching(
     90      const std::vector<Vector> &_old,
     91      const std::vector<Vector> &_new,
     92      const IndexList_t &_Matching);
     93
     94  static Polygon_t removeMatchingPoints(
     95      const VectorArray_t &_points,
     96      const IndexList_t &_matchingindices
     97      );
     98
     99  struct MatchingControlStructure {
     100    bool foundflag;
     101    double bestL2;
     102    IndexList_t bestmatching;
     103    VectorArray_t oldpoints;
     104    VectorArray_t newpoints;
     105  };
     106
     107  static void recurseMatchings(
     108      MatchingControlStructure &_MCS,
     109      IndexList_t &_matching,
     110      IndexList_t _indices,
     111      unsigned int _matchingsize);
     112
     113  static IndexList_t findBestMatching(
     114      const Polygon_t &_polygon,
     115      const Polygon_t &_newpolygon
     116      );
     117
     118  static Rotation_t findPlaneAligningRotation(
     119      const VectorArray_t &_referencepositions,
     120      const VectorArray_t &_currentpositions,
     121      const IndexList_t &_bestmatching
     122      );
     123
     124  static Rotation_t findPointAligningRotation(
     125      const VectorArray_t &remainingold,
     126      const VectorArray_t &remainingnew,
     127      const IndexList_t &_bestmatching);
     128
    78129};
    79130

  • src/Fragmentation/Exporters/unittests/Makefile.am

    r226860 rbcf13b  
    2020        SphericalPointDistributionUnitTest
    2121
    22 XFAIL_TESTS += SphericalPointDistributionUnitTest
    2322TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
    2423check_PROGRAMS += $(FRAGMENTATIONEXPORTERSTESTS)

  • src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp

    r226860 rbcf13b  
    150150}
    151151
     152void perturbPolygon(
     153    SphericalPointDistribution::Polygon_t &_polygon,
     154    double _amplitude
     155    )
     156{
     157  for (SphericalPointDistribution::Polygon_t::iterator iter = _polygon.begin();
     158      iter != _polygon.end(); ++iter) {
     159    Vector perturber;
     160    perturber.GetOneNormalVector((*iter));
     161    perturber.Scale(_amplitude);
     162    *iter = *iter + perturber;
     163    (*iter).Normalize();
     164  }
     165}
     166
     167static
     168bool areEqualToWithinBounds(
     169    const SphericalPointDistribution::Polygon_t &_polygon,
     170    const SphericalPointDistribution::Polygon_t &_otherpolygon,
     171    double _amplitude
     172    )
     173{
     174  // same size?
     175  if (_polygon.size() != _otherpolygon.size())
     176    return false;
     177  // same points ? We just check witrh trivial mapping, nothing fancy ...
     178  bool status = true;
     179  SphericalPointDistribution::Polygon_t::const_iterator iter = _polygon.begin();
     180  SphericalPointDistribution::Polygon_t::const_iterator otheriter = _otherpolygon.begin();
     181  for (; iter != _polygon.end(); ++iter, ++otheriter) {
     182    status &= (*iter - *otheriter).Norm() < _amplitude;
     183  }
     184  return status;
     185}
     186
     187/** UnitTest for areEqualToWithinBounds()
     188 */
     189void SphericalPointDistributionTest::areEqualToWithinBoundsTest()
     190{
     191  // test with no points
     192  {
     193    SphericalPointDistribution::Polygon_t polygon;
     194    SphericalPointDistribution::Polygon_t expected = polygon;
     195    CPPUNIT_ASSERT( areEqualToWithinBounds(polygon, expected, std::numeric_limits<double>::epsilon()*1e2) );
     196  }
     197  // test with one point
     198  {
     199    SphericalPointDistribution::Polygon_t polygon;
     200    polygon += Vector(1.,0.,0.);
     201    SphericalPointDistribution::Polygon_t expected = polygon;
     202    CPPUNIT_ASSERT( areEqualToWithinBounds(polygon, expected, std::numeric_limits<double>::epsilon()*1e2) );
     203  }
     204  // test with two points
     205  {
     206    SphericalPointDistribution::Polygon_t polygon;
     207    polygon += Vector(1.,0.,0.);
     208    polygon += Vector(0.,1.,0.);
     209    SphericalPointDistribution::Polygon_t expected = polygon;
     210    CPPUNIT_ASSERT( areEqualToWithinBounds(polygon, expected, std::numeric_limits<double>::epsilon()*1e2) );
     211  }
     212
     213  // test with two points in different order: THIS GOES WRONG: We only check trivially
     214  {
     215    SphericalPointDistribution::Polygon_t polygon;
     216    polygon += Vector(1.,0.,0.);
     217    polygon += Vector(0.,1.,0.);
     218    SphericalPointDistribution::Polygon_t expected;
     219    expected += Vector(0.,1.,0.);
     220    expected += Vector(1.,0.,0.);
     221    CPPUNIT_ASSERT( !areEqualToWithinBounds(polygon, expected, std::numeric_limits<double>::epsilon()*1e2) );
     222  }
     223
     224  // test with two different points
     225  {
     226    SphericalPointDistribution::Polygon_t polygon;
     227    polygon += Vector(1.,0.,0.);
     228    polygon += Vector(0.,1.,0.);
     229    SphericalPointDistribution::Polygon_t expected;
     230    expected += Vector(1.01,0.,0.);
     231    expected += Vector(0.,1.,0.);
     232    CPPUNIT_ASSERT( areEqualToWithinBounds(polygon, expected, 0.05) );
     233    CPPUNIT_ASSERT( !areEqualToWithinBounds(polygon, expected, 0.005) );
     234  }
     235
     236  // test with different number of points
     237  {
     238    SphericalPointDistribution::Polygon_t polygon;
     239    polygon += Vector(1.,0.,0.);
     240    polygon += Vector(0.,1.,0.);
     241    SphericalPointDistribution::Polygon_t expected;
     242    expected += Vector(0.,1.,0.);
     243    CPPUNIT_ASSERT( !areEqualToWithinBounds(polygon, expected, 0.05) );
     244  }
     245}
     246
     247
    152248/** UnitTest for matchSphericalPointDistributions() with three points
    153249 */
    … …  
    205301            newpolygon);
    206302    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     303    // also slightly perturbed
     304    const double amplitude = 0.05;
     305    perturbPolygon(polygon, amplitude);
     306    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    207307  }
    208308
    … …  
    227327            newpolygon);
    228328    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     329    // also slightly perturbed
     330    const double amplitude = 0.05;
     331    perturbPolygon(polygon, amplitude);
     332    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    229333  }
    230334
    … …  
    241345            newpolygon);
    242346    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     347    // also slightly perturbed
     348    const double amplitude = 0.05;
     349    perturbPolygon(polygon, amplitude);
     350    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    243351  }
    244352
    … …  
    260368            newpolygon);
    261369    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     370    // also slightly perturbed
     371    const double amplitude = 0.05;
     372    perturbPolygon(polygon, amplitude);
     373    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    262374  }
    263375}
    … …  
    318430            newpolygon);
    319431    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     432    // also slightly perturbed
     433    const double amplitude = 0.05;
     434    perturbPolygon(polygon, amplitude);
     435    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
     436  }
     437
     438  // test with two points, matching trivially, also with slightly perturbed
     439  {
     440    SphericalPointDistribution::Polygon_t polygon;
     441    polygon += Vector(1.,0.,0.), Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.);
     442    SphericalPointDistribution::Polygon_t newpolygon =
     443        SPD.get<4>();
     444    SphericalPointDistribution::Polygon_t expected = newpolygon;
     445    expected.pop_front(); // remove first point
     446    expected.pop_front(); // remove second point
     447    SphericalPointDistribution::Polygon_t remaining =
     448        SphericalPointDistribution::matchSphericalPointDistributions(
     449            polygon,
     450            newpolygon);
     451    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     452    // also slightly perturbed
     453    const double amplitude = 0.05;
     454    perturbPolygon(polygon, amplitude);
     455    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    320456  }
    321457
    … …  
    340476            newpolygon);
    341477    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     478    // also slightly perturbed
     479    const double amplitude = 0.05;
     480    perturbPolygon(polygon, amplitude);
     481    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    342482  }
    343483
    … …  
    360500            newpolygon);
    361501    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     502    // also slightly perturbed
     503    const double amplitude = 0.05;
     504    perturbPolygon(polygon, amplitude);
     505    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    362506  }
    363507
    … …  
    384528            newpolygon);
    385529    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     530    // also slightly perturbed
     531    const double amplitude = 0.05;
     532    perturbPolygon(polygon, amplitude);
     533    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    386534  }
    387535}
    … …  
    442590            newpolygon);
    443591    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     592    // also slightly perturbed
     593    const double amplitude = 0.05;
     594    perturbPolygon(polygon, amplitude);
     595    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    444596  }
    445597
    … …  
    494646            newpolygon);
    495647    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     648    // also slightly perturbed
     649    const double amplitude = 0.05;
     650    perturbPolygon(polygon, amplitude);
     651    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    496652  }
    497653
    … …  
    518674            newpolygon);
    519675    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     676    // also slightly perturbed
     677    const double amplitude = 0.05;
     678    perturbPolygon(polygon, amplitude);
     679    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    520680  }
    521681}
    … …  
    576736            newpolygon);
    577737    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     738    // also slightly perturbed
     739    const double amplitude = 0.05;
     740    perturbPolygon(polygon, amplitude);
     741    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    578742  }
    579743
    … …  
    628792            newpolygon);
    629793    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     794    // also slightly perturbed
     795    const double amplitude = 0.05;
     796    perturbPolygon(polygon, amplitude);
     797    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    630798  }
    631799
    … …  
    652820            newpolygon);
    653821    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     822    // also slightly perturbed
     823    const double amplitude = 0.05;
     824    perturbPolygon(polygon, amplitude);
     825    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    654826  }
    655827}
    … …  
    710882            newpolygon);
    711883    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     884    // also slightly perturbed
     885    const double amplitude = 0.05;
     886    perturbPolygon(polygon, amplitude);
     887    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    712888  }
    713889
    … …  
    762938            newpolygon);
    763939    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     940    // also slightly perturbed
     941    const double amplitude = 0.05;
     942    perturbPolygon(polygon, amplitude);
     943    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    764944  }
    765945
    … …  
    786966            newpolygon);
    787967    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     968    // also slightly perturbed
     969    const double amplitude = 0.05;
     970    perturbPolygon(polygon, amplitude);
     971    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    788972  }
    789973}
    … …  
    8441028            newpolygon);
    8451029    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1030    // also slightly perturbed
     1031    const double amplitude = 0.05;
     1032    perturbPolygon(polygon, amplitude);
     1033    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    8461034  }
    8471035
    … …  
    9001088            newpolygon);
    9011089    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1090    // also slightly perturbed
     1091    const double amplitude = 0.05;
     1092    perturbPolygon(polygon, amplitude);
     1093    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    9021094  }
    9031095
    … …  
    9241116            newpolygon);
    9251117    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1118    // also slightly perturbed
     1119    const double amplitude = 0.05;
     1120    perturbPolygon(polygon, amplitude);
     1121    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
    9261122  }
    9271123}

  • src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp

    r226860 rbcf13b  
    2222{
    2323    CPPUNIT_TEST_SUITE( SphericalPointDistributionTest) ;
     24    CPPUNIT_TEST ( areEqualToWithinBoundsTest );
    2425    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_2 );
    2526    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_3 );
    … …  
    3435      void setUp();
    3536      void tearDown();
     37      void areEqualToWithinBoundsTest();
    3638      void matchSphericalPointDistributionsTest_2();
    3739      void matchSphericalPointDistributionsTest_3();
Note: See TracChangeset for help on using the changeset viewer.