Ignore:
Timestamp:
Sep 12, 2016, 11:48:34 PM (9 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable
Children:
2199c2
Parents:
158ecb
git-author:
Frederik Heber <heber@…> (06/12/14 07:23:12)
git-committer:
Frederik Heber <heber@…> (09/12/16 23:48:34)
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

    r158ecb rb67d89  
    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#ifndef NDEBUG
     326  bool dontcheck = false;
     327#endif
     328  // initialize to no rotation
     329  Rotation_t Rotation;
     330  Rotation.first.Zero();
     331  Rotation.first[0] = 1.;
     332  Rotation.second = 0.;
     333
     334  // calculate center of triangle/line/point consisting of first points of matching
     335  Vector oldCenter;
     336  Vector newCenter;
     337  calculateOldAndNewCenters(
     338      oldCenter, newCenter,
     339      _referencepositions, _currentpositions, _bestmatching);
     340
     341  if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
     342    LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
     343    oldCenter.Normalize();
     344    newCenter.Normalize();
     345    if (!oldCenter.IsEqualTo(newCenter)) {
     346      // calculate rotation axis and angle
     347      Rotation.first = oldCenter;
     348      Rotation.first.VectorProduct(newCenter);
     349      Rotation.second = oldCenter.Angle(newCenter); // /(M_PI/2.);
     350    } else {
     351      // no rotation required anymore
     352    }
     353  } else {
     354    LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
     355    if ((oldCenter.IsZero()) && (newCenter.IsZero())) {
     356      // either oldCenter or newCenter (or both) is directly at origin
     357      if (_bestmatching.size() == 2) {
     358        // line case
     359        Vector oldPosition = _currentpositions[*_bestmatching.begin()];
     360        Vector newPosition = _referencepositions[0];
     361        // check whether we need to rotate at all
     362        if (!oldPosition.IsEqualTo(newPosition)) {
     363          Rotation.first = oldPosition;
     364          Rotation.first.VectorProduct(newPosition);
     365          // orientation will fix the sign here eventually
     366          Rotation.second = oldPosition.Angle(newPosition);
     367        } else {
     368          // no rotation required anymore
     369        }
     370      } else {
     371        // triangle case
     372        // both triangles/planes have same center, hence get axis by
     373        // VectorProduct of Normals
     374        Plane newplane(_referencepositions[0], _referencepositions[1], _referencepositions[2]);
     375        VectorArray_t vectors;
     376        for (IndexList_t::const_iterator iter = _bestmatching.begin();
     377            iter != _bestmatching.end(); ++iter)
     378          vectors.push_back(_currentpositions[*iter]);
     379        Plane oldplane(vectors[0], vectors[1], vectors[2]);
     380        Vector oldPosition = oldplane.getNormal();
     381        Vector newPosition = newplane.getNormal();
     382        // check whether we need to rotate at all
     383        if (!oldPosition.IsEqualTo(newPosition)) {
     384          Rotation.first = oldPosition;
     385          Rotation.first.VectorProduct(newPosition);
     386          Rotation.first.Normalize();
     387
     388          // construct reference vector to determine direction of rotation
     389          const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
     390          Rotation.second = sign * oldPosition.Angle(newPosition);
     391          LOG(5, "DEBUG: Rotating plane normals by " << Rotation.second
     392              << " around axis " << Rotation.first);
     393        } else {
     394          // else do nothing
     395        }
     396      }
     397    } else {
     398      // TODO: we can't do anything here, but this case needs to be dealt with when
     399      // we have no ideal geometries anymore
     400      if ((oldCenter-newCenter).Norm() > warn_amplitude)
     401        ELOG(2, "oldCenter is " << oldCenter << ", yet newCenter is " << newCenter);
     402#ifndef NDEBUG
     403      // else they are considered close enough
     404      dontcheck = true;
     405#endif
     406    }
     407  }
     408
     409#ifndef NDEBUG
     410  // check: rotation brings newCenter onto oldCenter position
     411  if (!dontcheck) {
     412    Line Axis(zeroVec, Rotation.first);
     413    Vector test = Axis.rotateVector(newCenter, Rotation.second);
     414    LOG(4, "CHECK: rotated newCenter is " << test
     415        << ", oldCenter is " << oldCenter);
     416    ASSERT( (test - oldCenter).NormSquared() < std::numeric_limits<double>::epsilon()*1e4,
     417        "matchSphericalPointDistributions() - rotation does not work as expected by "
     418        +toString((test - oldCenter).NormSquared())+".");
     419  }
     420#endif
     421
     422  return Rotation;
     423}
     424
     425SphericalPointDistribution::Rotation_t SphericalPointDistribution::findPointAligningRotation(
     426    const VectorArray_t &remainingold,
     427    const VectorArray_t &remainingnew,
     428    const IndexList_t &_bestmatching)
     429{
     430  // initialize rotation to zero
     431  Rotation_t Rotation;
     432  Rotation.first.Zero();
     433  Rotation.first[0] = 1.;
     434  Rotation.second = 0.;
     435
     436  // recalculate center
     437  Vector oldCenter;
     438  Vector newCenter;
     439  calculateOldAndNewCenters(
     440      oldCenter, newCenter,
     441      remainingold, remainingnew, _bestmatching);
     442
     443  Vector oldPosition = remainingnew[*_bestmatching.begin()];
     444  Vector newPosition = remainingold[0];
     445  LOG(6, "DEBUG: oldPosition is " << oldPosition << " and newPosition is " << newPosition);
     446  if (!oldPosition.IsEqualTo(newPosition)) {
     447    if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
     448      oldCenter.Normalize();  // note weighted sum of normalized weight is not normalized
     449      Rotation.first = oldCenter;
     450      LOG(6, "DEBUG: Picking normalized oldCenter as Rotation.first " << oldCenter);
     451      oldPosition.ProjectOntoPlane(Rotation.first);
     452      newPosition.ProjectOntoPlane(Rotation.first);
     453      LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
     454    } else {
     455      if (_bestmatching.size() == 2) {
     456        // line situation
     457        try {
     458          Plane oldplane(oldPosition, oldCenter, newPosition);
     459          Rotation.first = oldplane.getNormal();
     460          LOG(6, "DEBUG: Plane is " << oldplane << " and normal is " << Rotation.first);
     461        } catch (LinearDependenceException &e) {
     462          LOG(6, "DEBUG: Vectors defining plane are linearly dependent.");
     463          // oldPosition and newPosition are on a line, just flip when not equal
     464          if (!oldPosition.IsEqualTo(newPosition)) {
     465            Rotation.first.Zero();
     466            Rotation.first.GetOneNormalVector(oldPosition);
     467            LOG(6, "DEBUG: For flipping we use Rotation.first " << Rotation.first);
     468            assert( Rotation.first.ScalarProduct(oldPosition) < std::numeric_limits<double>::epsilon()*1e4);
     469  //              Rotation.second = M_PI;
     470          } else {
     471            LOG(6, "DEBUG: oldPosition and newPosition are equivalent.");
     472          }
     473        }
     474      } else {
     475        // triangle situation
     476        Plane oldplane(remainingold[0], remainingold[1], remainingold[2]);
     477        Rotation.first = oldplane.getNormal();
     478        LOG(6, "DEBUG: oldPlane is " << oldplane << " and normal is " << Rotation.first);
     479        oldPosition.ProjectOntoPlane(Rotation.first);
     480        LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
     481      }
     482    }
     483    // construct reference vector to determine direction of rotation
     484    const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
     485    Rotation.second = sign * oldPosition.Angle(newPosition);
     486  } else {
     487    LOG(6, "DEBUG: oldPosition and newPosition are equivalent, hence no orientating rotation.");
     488  }
     489
     490  return Rotation;
    256491}
    257492
    … …  
    269504      << " with new polygon " << _newpolygon);
    270505
     506  if (_polygon.size() == _newpolygon.size()) {
     507    // same number of points desired as are present? Do nothing
     508    LOG(2, "INFO: There are no vacant points to return.");
     509    return remainingreturnpolygon;
     510  }
     511
    271512  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
     513    IndexList_t bestmatching = findBestMatching(_polygon, _newpolygon);
     514    LOG(2, "INFO: Best matching is " << bestmatching);
     515
     516    // determine rotation angles to align the two point distributions with
     517    // respect to bestmatching:
     518    // we use the center between the three first matching points
     519    /// the first rotation brings these two centers to coincide
     520    VectorArray_t rotated_newpolygon = remainingnew;
    279521    {
    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);
     522      Rotation_t Rotation = findPlaneAligningRotation(
     523          remainingold,
     524          remainingnew,
     525          bestmatching);
     526      LOG(5, "DEBUG: Rotating coordinate system by " << Rotation.second
     527          << " around axis " << Rotation.first);
     528      Line Axis(zeroVec, Rotation.first);
     529
     530      // apply rotation angle to bring newCenter to oldCenter
     531      for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
     532          iter != rotated_newpolygon.end(); ++iter) {
     533        Vector &current = *iter;
     534        LOG(6, "DEBUG: Original point is " << current);
     535        current =  Axis.rotateVector(current, Rotation.second);
     536        LOG(6, "DEBUG: Rotated point is " << current);
     537      }
     538
     539#ifndef NDEBUG
     540      // check: rotated "newCenter" should now equal oldCenter
     541      {
     542        Vector oldCenter;
     543        Vector rotatednewCenter;
     544        calculateOldAndNewCenters(
     545            oldCenter, rotatednewCenter,
     546            remainingold, rotated_newpolygon, bestmatching);
     547        // NOTE: Center must not necessarily lie on the sphere with norm 1, hence, we
     548        // have to normalize it just as before, as oldCenter and newCenter lengths may differ.
     549        if ((!oldCenter.IsZero()) && (!rotatednewCenter.IsZero())) {
     550          oldCenter.Normalize();
     551          rotatednewCenter.Normalize();
     552          LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
     553              << ", oldCenter is " << oldCenter);
     554          ASSERT( (rotatednewCenter - oldCenter).NormSquared() < std::numeric_limits<double>::epsilon()*1e4,
     555              "matchSphericalPointDistributions() - rotation does not work as expected by "
     556              +toString((rotatednewCenter - oldCenter).NormSquared())+".");
    327557        }
    328558      }
    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);
     559#endif
     560    }
     561    /// the second (orientation) rotation aligns the planes such that the
     562    /// points themselves coincide
     563    if (bestmatching.size() > 1) {
     564      Rotation_t Rotation = findPointAligningRotation(
     565          remainingold,
     566          rotated_newpolygon,
     567          bestmatching);
     568
     569      // construct RotationAxis and two points on its plane, defining the angle
     570      Rotation.first.Normalize();
     571      const Line RotationAxis(zeroVec, Rotation.first);
     572
     573      LOG(5, "DEBUG: Rotating around self is " << Rotation.second
     574          << " around axis " << RotationAxis);
    352575
    353576#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         }
     577      // check: first bestmatching in rotated_newpolygon and remainingnew
     578      // should now equal
     579      {
     580        const IndexList_t::const_iterator iter = bestmatching.begin();
     581        Vector rotatednew = RotationAxis.rotateVector(
     582            rotated_newpolygon[*iter],
     583            Rotation.second);
     584        LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
     585            << " while old was " << remainingold[0]);
     586        ASSERT( (rotatednew - remainingold[0]).Norm() < warn_amplitude,
     587            "matchSphericalPointDistributions() - orientation rotation ends up off by more than "
     588            +toString(warn_amplitude)+".");
     589      }
    367590#endif
    368591
    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         }
     592      for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
     593          iter != rotated_newpolygon.end(); ++iter) {
     594        Vector &current = *iter;
     595        LOG(6, "DEBUG: Original point is " << current);
     596        current = RotationAxis.rotateVector(current, Rotation.second);
     597        LOG(6, "DEBUG: Rotated point is " << current);
    376598      }
    377599    }
    … …  
    379601    // remove all points in matching and return remaining ones
    380602    SphericalPointDistribution::Polygon_t remainingpoints =
    381         removeMatchingPoints(rotated_newpolygon, MCS.bestmatching);
     603        removeMatchingPoints(rotated_newpolygon, bestmatching);
    382604    LOG(2, "INFO: Remaining points are " << remainingpoints);
    383605    return remainingpoints;

  • src/Fragmentation/Exporters/SphericalPointDistribution.hpp

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

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

    r158ecb rb67d89  
    2020        SphericalPointDistributionUnitTest
    2121
    22 XFAIL_TESTS += SphericalPointDistributionUnitTest
    2322TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
    2423check_PROGRAMS += $(FRAGMENTATIONEXPORTERSTESTS)

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

    r158ecb rb67d89  
    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

    r158ecb rb67d89  
    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.