Changeset 3da643


Ignore:
Timestamp:
Aug 20, 2014, 1:04:08 PM (11 years ago)
Author:
Frederik Heber <heber@…>
Children:
0710bf
Parents:
6aa6b7
git-author:
Frederik Heber <heber@…> (06/12/14 07:23:12)
git-committer:
Frederik Heber <heber@…> (08/20/14 13:04:08)
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.
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/SphericalPointDistribution.cpp

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

  • src/Fragmentation/Exporters/SphericalPointDistribution.hpp

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

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

    r6aa6b7 r3da643  
    1717        SphericalPointDistributionUnitTest
    1818
    19 XFAIL_TESTS += SphericalPointDistributionUnitTest
    2019TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
    2120check_PROGRAMS += $(FRAGMENTATIONEXPORTERSTESTS)

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

    r6aa6b7 r3da643  
    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

    r6aa6b7 r3da643  
    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();

  • tests/regression/Fragmentation/FragmentMolecule/testsuite-fragmentation-fragment-molecule.at

    r6aa6b7 r3da643  
    3737AT_SETUP([Fragmentation - Fragmentation with cycles])
    3838AT_KEYWORDS([fragmentation fragment-molecule cycle])
    39 AT_XFAIL_IF([/bin/true])
    4039
    4140file=benzene.pdb

  • tests/regression/Fragmentation/StoreSaturatedFragment/testsuite-fragmentation-store-saturated-fragment.at

    r6aa6b7 r3da643  
    2121AT_SETUP([Fragmentation - Store saturated fragment])
    2222AT_KEYWORDS([fragmentation store-saturated-fragment])
    23 AT_XFAIL_IF([/bin/true])
    2423
    2524file=BondFragment0.xyz
Note: See TracChangeset for help on using the changeset viewer.