Changeset a2f8a9 for src


Ignore:
Timestamp:
Aug 20, 2014, 1:06:16 PM (11 years ago)
Author:
Frederik Heber <heber@…>
Children:
ef3885
Parents:
1cde4e8
git-author:
Frederik Heber <heber@…> (07/09/14 22:08:37)
git-committer:
Frederik Heber <heber@…> (08/20/14 13:06:16)
Message:

Added calculation of center of minimum distance by bisection.

  • this will give us a unique a definite point independent of the (rotational) position of the point set on the unit sphere.
  • added unit test.
Location:
src/Fragmentation/Exporters
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/SphericalPointDistribution.cpp

    r1cde4e8 ra2f8a9  
    4343
    4444#include <algorithm>
    45 #include <boost/math/quaternion.hpp>
     45#include <boost/assign.hpp>
    4646#include <cmath>
    4747#include <functional>
    … …  
    5858#include "LinearAlgebra/Vector.hpp"
    5959
     60using namespace boost::assign;
     61
    6062// static entities
    6163const double SphericalPointDistribution::warn_amplitude = 1e-2;
    … …  
    8385 */
    8486inline
    85 Vector calculateCenter(
     87Vector calculateGeographicMidpoint(
    8688    const SphericalPointDistribution::VectorArray_t &_positions,
    8789    const SphericalPointDistribution::IndexList_t &_indices)
    … …  
    9698
    9799  return Center;
     100}
     101
     102inline
     103double calculateMinimumDistance(
     104    const Vector &_center,
     105    const SphericalPointDistribution::VectorArray_t &_points,
     106    const SphericalPointDistribution::IndexList_t & _indices)
     107{
     108  double MinimumDistance = 0.;
     109  for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
     110      iter != _indices.end(); ++iter) {
     111    const double angle = _center.Angle(_points[*iter]);
     112    MinimumDistance += angle*angle;
     113  }
     114  return sqrt(MinimumDistance);
     115}
     116
     117/** Calculates the center of minimum distance for a given set of points \a _points.
     118 *
     119 * According to http://www.geomidpoint.com/calculation.html this goes a follows:
     120 * -# Let CurrentPoint be the geographic midpoint found in Method A. this is used as the starting point for the search.
     121 * -# Let MinimumDistance be the sum total of all distances from the current point to all locations in 'Your Places'.
     122 * -# Find the total distance between each location in 'Your Places' and all other locations in 'Your Places'. If any one of these locations has a new smallest distance then that location becomes the new CurrentPoint and MinimumDistance.
     123 * -# Let TestDistance be PI/2 radians (6225 miles or 10018 km).
     124 * -# Find the total distance between each of 8 test points and all locations in 'Your Places'. The test points are arranged in a circular pattern around the CurrentPoint at a distance of TestDistance to the north, northeast, east, southeast, south, southwest, west and northwest.
     125 * -# If any of these 8 points has a new smallest distance then that point becomes the new CurrentPoint and MinimumDistance and go back to step 5 using that point.
     126 * -# If none of the 8 test points has a new smallest distance then divide TestDistance by 2 and go back to step 5 using the same point.
     127 * -# Repeat steps 5 to 7 until no new smallest distance can be found or until TestDistance is less than 0.00000002 radians.
     128 *
     129 * \param _points set of points
     130 * \return Center of minimum distance for all these points, is always of length 1
     131 */
     132Vector SphericalPointDistribution::calculateCenterOfMinimumDistance(
     133    const SphericalPointDistribution::VectorArray_t &_positions,
     134    const SphericalPointDistribution::IndexList_t &_indices)
     135{
     136  ASSERT( _positions.size() >= _indices.size(),
     137      "calculateCenterOfMinimumDistance() - less positions than indices given.");
     138  Vector center(1.,0.,0.);
     139
     140  /// first treat some special cases
     141  // no positions given: return x axis vector (NOT zero!)
     142  {
     143    if (_indices.empty())
     144      return center;
     145    // one position given: return it directly
     146    if (_positions.size() == (size_t)1)
     147      return _positions[0];
     148    // two positions on a line given: return closest point to (1.,0.,0.)
     149    if (fabs(_positions[0].ScalarProduct(_positions[1]) + 1.)
     150        < std::numeric_limits<double>::epsilon()*1e4) {
     151      Vector candidate;
     152      if (_positions[0].ScalarProduct(center) > _positions[1].ScalarProduct(center))
     153        candidate = _positions[0];
     154      else
     155        candidate = _positions[1];
     156      // non-uniqueness: all positions on great circle, normal to given line are valid
     157      // so, we just pick one because returning a unique point is topmost priority
     158      Vector normal;
     159      normal.GetOneNormalVector(candidate);
     160      Vector othernormal = candidate;
     161      othernormal.VectorProduct(normal);
     162      // now both normal and othernormal describe the plane containing the great circle
     163      Plane greatcircle(normal, zeroVec, othernormal);
     164      // check which axis is contained and pick the one not
     165      if (greatcircle.isContained(center)) {
     166        center = Vector(0.,1.,0.);
     167        if (greatcircle.isContained(center))
     168          center = Vector(0.,0.,1.);
     169        // now we are done cause a plane cannot contain all three axis vectors
     170      }
     171      center = greatcircle.getClosestPoint(candidate);
     172      // assure length of 1
     173      center.Normalize();
     174    }
     175  }
     176
     177  // start with geographic midpoint
     178  center = calculateGeographicMidpoint(_positions, _indices);
     179  if (!center.IsZero()) {
     180    center.Normalize();
     181    LOG(4, "DEBUG: Starting with geographical midpoint of " << _positions << " under indices "
     182        << _indices << " is " << center);
     183  } else {
     184    // any point is good actually
     185    center = _positions[0];
     186    LOG(4, "DEBUG: Starting with first position " << center);
     187  }
     188
     189  // calculate initial MinimumDistance
     190  double MinimumDistance = calculateMinimumDistance(center, _positions, _indices);
     191  LOG(5, "DEBUG: MinimumDistance to this center is " << MinimumDistance);
     192
     193  // check all present points whether they may serve as a better center
     194  for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
     195      iter != _indices.end(); ++iter) {
     196    const Vector &centerCandidate = _positions[*iter];
     197    const double candidateDistance = calculateMinimumDistance(centerCandidate, _positions, _indices);
     198    if (candidateDistance < MinimumDistance) {
     199      MinimumDistance = candidateDistance;
     200      center = centerCandidate;
     201      LOG(5, "DEBUG: new MinimumDistance to current test point " << MinimumDistance
     202          << " is " << center);
     203    }
     204  }
     205  LOG(5, "DEBUG: new MinimumDistance to center " << center << " is " << MinimumDistance);
     206
     207  // now iterate over TestDistance
     208  double TestDistance = Vector(1.,0.,0.).Angle(Vector(0.,1.,0.));
     209//  LOG(6, "DEBUG: initial TestDistance is " << TestDistance);
     210
     211  const double threshold = sqrt(std::numeric_limits<double>::epsilon());
     212  // check each of eight test points at N, NE, E, SE, S, SW, W, NW
     213  typedef std::vector<double> angles_t;
     214  angles_t testangles;
     215  testangles += 0./180.*M_PI, 45./180.*M_PI, 90./180.*M_PI, 135./180.*M_PI,
     216      180./180.*M_PI, 225./180.*M_PI, 270./180.*M_PI, 315./180.*M_PI;
     217  while (TestDistance > threshold) {
     218    Vector OneNormal;
     219    OneNormal.GetOneNormalVector(center);
     220    Line RotationAxis(zeroVec, OneNormal);
     221    Vector North = RotationAxis.rotateVector(center,TestDistance);
     222    Line CompassRose(zeroVec, center);
     223    bool updatedflag = false;
     224    for (angles_t::const_iterator angleiter = testangles.begin();
     225        angleiter != testangles.end(); ++angleiter) {
     226      Vector centerCandidate = CompassRose.rotateVector(North, *angleiter);
     227//      centerCandidate.Normalize();
     228      const double candidateDistance = calculateMinimumDistance(centerCandidate, _positions, _indices);
     229      if (candidateDistance < MinimumDistance) {
     230        MinimumDistance = candidateDistance;
     231        center = centerCandidate;
     232        updatedflag = true;
     233        LOG(5, "DEBUG: new MinimumDistance to test point at " << *angleiter/M_PI*180.
     234            << "° is " << MinimumDistance);
     235      }
     236    }
     237
     238    // if no new point, decrease TestDistance
     239    if (!updatedflag) {
     240      TestDistance *= 0.5;
     241//      LOG(6, "DEBUG: TestDistance is now " << TestDistance);
     242    }
     243  }
     244  LOG(4, "DEBUG: Final MinimumDistance to center " << center << " is " << MinimumDistance);
     245
     246  return center;
     247}
     248
     249Vector calculateCenterOfMinimumDistance(
     250    const SphericalPointDistribution::PolygonWithIndices &_points)
     251{
     252  return SphericalPointDistribution::calculateCenterOfMinimumDistance(_points.polygon, _points.indices);
     253}
     254
     255/** Calculate the center of a given set of points in \a _positions but only
     256 * for those indicated by \a _indices.
     257 *
     258 */
     259inline
     260Vector calculateCenter(
     261    const SphericalPointDistribution::VectorArray_t &_positions,
     262    const SphericalPointDistribution::IndexList_t &_indices)
     263{
     264//  Vector Center;
     265//  Center.Zero();
     266//  for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
     267//      iter != _indices.end(); ++iter)
     268//    Center += _positions[*iter];
     269//  if (!_indices.empty())
     270//    Center *= 1./(double)_indices.size();
     271//
     272//  return Center;
     273  return SphericalPointDistribution::calculateCenterOfMinimumDistance(_positions, _indices);
    98274}
    99275
    … …  
    744920            oldCenter, rotatednewCenter,
    745921            oldSet, rotatednewSet);
    746         // NOTE: Center must not necessarily lie on the sphere with norm 1, hence, we
    747         // have to normalize it just as before, as oldCenter and newCenter lengths may differ.
    748         if ((!oldCenter.IsZero()) && (!rotatednewCenter.IsZero())) {
    749           oldCenter.Normalize();
    750           rotatednewCenter.Normalize();
    751           LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
    752               << ", oldCenter is " << oldCenter);
    753           const double difference = (rotatednewCenter - oldCenter).NormSquared();
    754           ASSERT( difference < std::numeric_limits<double>::epsilon()*1e4,
    755               "matchSphericalPointDistributions() - rotation does not work as expected by "
    756               +toString(difference)+".");
    757         }
     922        oldCenter.Normalize();
     923        rotatednewCenter.Normalize();
     924        // check whether centers are anti-parallel (factor -1)
     925        // then we have the "non-unique poles" situation: points lie on great circle
     926        // and both poles are valid solution
     927        if (fabs(oldCenter.ScalarProduct(rotatednewCenter) + 1.)
     928            < std::numeric_limits<double>::epsilon()*1e4)
     929          rotatednewCenter *= -1.;
     930        LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
     931            << ", oldCenter is " << oldCenter);
     932        const double difference = (rotatednewCenter - oldCenter).NormSquared();
     933        ASSERT( difference < std::numeric_limits<double>::epsilon()*1e4,
     934            "matchSphericalPointDistributions() - rotation does not work as expected by "
     935            +toString(difference)+".");
    758936      }
    759937#endif

  • src/Fragmentation/Exporters/SphericalPointDistribution.hpp

    r1cde4e8 ra2f8a9  
    121121  };
    122122
     123  static Vector calculateCenterOfMinimumDistance(
     124      const SphericalPointDistribution::VectorArray_t &_positions,
     125      const SphericalPointDistribution::IndexList_t &_indices);
     126
    123127private:
    124128  //!> grant unit tests access to private parts
    … …  
    170174      const PolygonWithIndices &remainingold,
    171175      const PolygonWithIndices &remainingnew);
    172 
    173176};
    174177

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

    r1cde4e8 ra2f8a9  
    6868CPPUNIT_TEST_SUITE_REGISTRATION( SphericalPointDistributionTest );
    6969
     70/** due to root-taking in function we only have limited numerical precision,
     71 * basically half of the double range.
     72 */
     73const double CenterAccuracy = sqrt(std::numeric_limits<double>::epsilon()*1e2);
    7074
    7175void SphericalPointDistributionTest::setUp()
    … …  
    8286}
    8387
    84 
    85 /** UnitTest for matchSphericalPointDistributions() with two points
     88/** UnitTest for calculateCenterOfMinimumDistance()
    8689 */
    87 void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_2()
     90void SphericalPointDistributionTest::calculateCenterOfMinimumDistanceTest()
    8891{
    89   SphericalPointDistribution SPD(1.);
    90   // test with one point, matching trivially
    91   {
    92     SphericalPointDistribution::WeightedPolygon_t polygon;
    93     polygon += std::make_pair(Vector(1.,0.,0.), 1);
    94     SphericalPointDistribution::Polygon_t newpolygon =
    95         SPD.get<2>();
    96     SphericalPointDistribution::Polygon_t expected;
    97     expected += Vector(-1.,0.,0.);
    98     SphericalPointDistribution::Polygon_t remaining =
    99         SphericalPointDistribution::matchSphericalPointDistributions(
    100             polygon,
    101             newpolygon);
    102     CPPUNIT_ASSERT_EQUAL( expected, remaining );
    103   }
    104 
    105   // test with one point, just a flip of axis
    106   {
    107     SphericalPointDistribution::WeightedPolygon_t polygon;
    108     polygon += std::make_pair( Vector(0.,1.,0.), 1);
    109     SphericalPointDistribution::Polygon_t newpolygon =
    110         SPD.get<2>();
    111     SphericalPointDistribution::Polygon_t expected;
    112     expected += Vector(0.,-1.,0.);
    113     SphericalPointDistribution::Polygon_t remaining =
    114         SphericalPointDistribution::matchSphericalPointDistributions(
    115             polygon,
    116             newpolygon);
    117     CPPUNIT_ASSERT_EQUAL( expected, remaining );
    118   }
    119 
    120   // test with one point, just a flip to another axis
    121   {
    122     SphericalPointDistribution::WeightedPolygon_t polygon;
    123     polygon += std::make_pair( Vector(0.,0.,-1.), 1);
    124     SphericalPointDistribution::Polygon_t newpolygon =
    125         SPD.get<2>();
    126     SphericalPointDistribution::Polygon_t expected;
    127     expected += Vector(0.,0.,1.);
    128     SphericalPointDistribution::Polygon_t remaining =
    129         SphericalPointDistribution::matchSphericalPointDistributions(
    130             polygon,
    131             newpolygon);
    132     CPPUNIT_ASSERT_EQUAL( expected, remaining );
    133   }
    134 
    135   // test with one point, full rotation
    136   {
    137     Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    138     SphericalPointDistribution::WeightedPolygon_t polygon;
    139     polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    140     SphericalPointDistribution::Polygon_t newpolygon =
    141         SPD.get<2>();
    142     SphericalPointDistribution::Polygon_t expected;
    143     expected += RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI);
    144     SphericalPointDistribution::Polygon_t remaining =
    145         SphericalPointDistribution::matchSphericalPointDistributions(
    146             polygon,
    147             newpolygon);
    148     CPPUNIT_ASSERT_EQUAL( expected, remaining );
    149   }
    150 }
    151 
    152 void perturbPolygon(
    153     SphericalPointDistribution::WeightedPolygon_t &_polygon,
    154     double _amplitude
    155     )
    156 {
    157   for (SphericalPointDistribution::WeightedPolygon_t::iterator iter = _polygon.begin();
    158       iter != _polygon.end(); ++iter) {
    159     Vector perturber;
    160     perturber.GetOneNormalVector(iter->first);
    161     perturber.Scale(_amplitude);
    162     iter->first = iter->first + perturber;
    163     (iter->first).Normalize();
     92  // single point
     93  {
     94    SphericalPointDistribution::VectorArray_t points;
     95    points +=
     96        Vector(1.,0.,0.);
     97    SphericalPointDistribution::IndexList_t indices;
     98    indices += 0;
     99    const Vector expected = points[0];
     100    const Vector center =
     101        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
     102//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
     103//    CPPUNIT_ASSERT_EQUAL ( expected, center );
     104    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
     105  }
     106
     107  // single point, rotated
     108  {
     109    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
     110    SphericalPointDistribution::VectorArray_t points;
     111    points +=
     112        RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI);
     113    SphericalPointDistribution::IndexList_t indices;
     114    indices += 0;
     115    const Vector expected = points[0];
     116    const Vector center =
     117        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
     118//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
     119//    CPPUNIT_ASSERT_EQUAL ( expected, center );
     120    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
     121  }
     122
     123  // two points
     124  {
     125    SphericalPointDistribution::VectorArray_t points;
     126    points +=
     127        Vector(1.,0.,0.),
     128        Vector(0.,1.,0.);
     129    SphericalPointDistribution::IndexList_t indices;
     130    indices += 0,1;
     131    const Vector expected = Vector(M_SQRT1_2,M_SQRT1_2,0.);
     132    const Vector center =
     133        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
     134//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
     135//    CPPUNIT_ASSERT_EQUAL ( expected, center );
     136    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
     137  }
     138
     139  // two points, rotated
     140  {
     141    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
     142    SphericalPointDistribution::VectorArray_t points;
     143    points +=
     144        RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI),
     145        RotationAxis.rotateVector(Vector(0.,1.,0.), 47.6/180*M_PI);
     146    SphericalPointDistribution::IndexList_t indices;
     147    indices += 0,1;
     148    const Vector expected = RotationAxis.rotateVector(Vector(M_SQRT1_2,M_SQRT1_2,0.), 47.6/180*M_PI);
     149    const Vector center =
     150        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
     151//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
     152//    CPPUNIT_ASSERT_EQUAL ( expected, center );
     153    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
     154  }
     155
     156  // three points in line
     157  {
     158    SphericalPointDistribution::VectorArray_t points;
     159    points +=
     160        Vector(1.,0.,0.),
     161        Vector(0.,1.,0.),
     162        Vector(-1.,0.,0.);
     163    SphericalPointDistribution::IndexList_t indices;
     164    indices += 0,1,2;
     165    const Vector expected = points[1];
     166    const Vector center =
     167        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
     168//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
     169//    CPPUNIT_ASSERT_EQUAL ( expected, center );
     170    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
     171  }
     172
     173  // three points in line, rotated
     174  {
     175    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
     176    SphericalPointDistribution::VectorArray_t points;
     177    points +=
     178        RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI),
     179        RotationAxis.rotateVector(Vector(0.,1.,0.), 47.6/180*M_PI),
     180        RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI);
     181    SphericalPointDistribution::IndexList_t indices;
     182    indices += 0,1,2;
     183    const Vector expected = points[1];
     184    const Vector center =
     185        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
     186//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
     187//    CPPUNIT_ASSERT_EQUAL ( expected, center );
     188    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
    164189  }
    165190}
    … …  
    180205  SphericalPointDistribution::Polygon_t::const_iterator otheriter = _otherpolygon.begin();
    181206  for (; iter != _polygon.end(); ++iter, ++otheriter) {
    182     status &= (*iter - *otheriter).Norm() < _amplitude;
     207    status &= (*iter).IsEqualTo(*otheriter, _amplitude);
    183208  }
    184209  return status;
    … …  
    242267    expected += Vector(0.,1.,0.);
    243268    CPPUNIT_ASSERT( !areEqualToWithinBounds(polygon, expected, 0.05) );
     269  }
     270}
     271
     272/** UnitTest for matchSphericalPointDistributions() with two points
     273 */
     274void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_2()
     275{
     276  SphericalPointDistribution SPD(1.);
     277  // test with one point, matching trivially
     278  {
     279    SphericalPointDistribution::WeightedPolygon_t polygon;
     280    polygon += std::make_pair(Vector(1.,0.,0.), 1);
     281    SphericalPointDistribution::Polygon_t newpolygon =
     282        SPD.get<2>();
     283    SphericalPointDistribution::Polygon_t expected;
     284    expected += Vector(-1.,0.,0.);
     285    SphericalPointDistribution::Polygon_t remaining =
     286        SphericalPointDistribution::matchSphericalPointDistributions(
     287            polygon,
     288            newpolygon);
     289//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     290    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
     291  }
     292
     293  // test with one point, just a flip of axis
     294  {
     295    SphericalPointDistribution::WeightedPolygon_t polygon;
     296    polygon += std::make_pair( Vector(0.,1.,0.), 1);
     297    SphericalPointDistribution::Polygon_t newpolygon =
     298        SPD.get<2>();
     299    SphericalPointDistribution::Polygon_t expected;
     300    expected += Vector(0.,-1.,0.);
     301    SphericalPointDistribution::Polygon_t remaining =
     302        SphericalPointDistribution::matchSphericalPointDistributions(
     303            polygon,
     304            newpolygon);
     305//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     306    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
     307  }
     308
     309  // test with one point, just a flip to another axis
     310  {
     311    SphericalPointDistribution::WeightedPolygon_t polygon;
     312    polygon += std::make_pair( Vector(0.,0.,-1.), 1);
     313    SphericalPointDistribution::Polygon_t newpolygon =
     314        SPD.get<2>();
     315    SphericalPointDistribution::Polygon_t expected;
     316    expected += Vector(0.,0.,1.);
     317    SphericalPointDistribution::Polygon_t remaining =
     318        SphericalPointDistribution::matchSphericalPointDistributions(
     319            polygon,
     320            newpolygon);
     321//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     322    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
     323  }
     324
     325  // test with one point, full rotation
     326  {
     327    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
     328    SphericalPointDistribution::WeightedPolygon_t polygon;
     329    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
     330    SphericalPointDistribution::Polygon_t newpolygon =
     331        SPD.get<2>();
     332    SphericalPointDistribution::Polygon_t expected;
     333    expected += RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI);
     334    SphericalPointDistribution::Polygon_t remaining =
     335        SphericalPointDistribution::matchSphericalPointDistributions(
     336            polygon,
     337            newpolygon);
     338//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     339    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
     340  }
     341}
     342
     343void perturbPolygon(
     344    SphericalPointDistribution::WeightedPolygon_t &_polygon,
     345    double _amplitude
     346    )
     347{
     348  for (SphericalPointDistribution::WeightedPolygon_t::iterator iter = _polygon.begin();
     349      iter != _polygon.end(); ++iter) {
     350    Vector perturber;
     351    perturber.GetOneNormalVector(iter->first);
     352    perturber.Scale(_amplitude);
     353    iter->first = iter->first + perturber;
     354    (iter->first).Normalize();
    244355  }
    245356}
    … …  
    367478            polygon,
    368479            newpolygon);
    369     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     480//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     481    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    370482  }
    371483
    … …  
    387499            polygon,
    388500            newpolygon);
    389     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     501//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     502    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    390503  }
    391504
    … …  
    404517            polygon,
    405518            newpolygon);
    406     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     519//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     520    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    407521    // also slightly perturbed
    408522    const double amplitude = 0.05;
    … …  
    429543            polygon,
    430544            newpolygon);
    431     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     545//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     546    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    432547    // also slightly perturbed
    433548    const double amplitude = 0.05;
    … …  
    449564            polygon,
    450565            newpolygon);
    451     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     566//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     567    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    452568    // also slightly perturbed
    453569    const double amplitude = 0.05;
    … …  
    471587            polygon,
    472588            newpolygon);
    473     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     589//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     590    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    474591    // also slightly perturbed
    475592    const double amplitude = 0.05;
    … …  
    497614            polygon,
    498615            newpolygon);
    499     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     616    //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     617        CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    500618  }
    501619
    … …  
    517635            polygon,
    518636            newpolygon);
    519     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     637//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     638    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    520639  }
    521640
    … …  
    534653            polygon,
    535654            newpolygon);
    536     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     655//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     656    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    537657    // also slightly perturbed
    538658    const double amplitude = 0.05;
    … …  
    555675            polygon,
    556676            newpolygon);
    557     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     677//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     678    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    558679    // also slightly perturbed
    559680    const double amplitude = 0.05;
    … …  
    580701            polygon,
    581702            newpolygon);
    582     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     703//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     704    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    583705    // also slightly perturbed
    584706    const double amplitude = 0.05;
    … …  
    603725            polygon,
    604726            newpolygon);
    605     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     727//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     728    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    606729    // also slightly perturbed
    607730    const double amplitude = 0.05;
    … …  
    630753            polygon,
    631754            newpolygon);
    632     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     755//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     756    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    633757    // also slightly perturbed
    634758    const double amplitude = 0.05;
    … …  
    659783            newpolygon);
    660784//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
    661     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     785//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     786    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    662787  }
    663788
    … …  
    677802            newpolygon);
    678803//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
    679     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     804//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     805    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    680806  }
    681807
    … …  
    696822            newpolygon);
    697823//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
    698     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     824//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     825    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    699826  }
    700827
    … …  
    714841            newpolygon);
    715842//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
    716     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     843//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     844    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    717845  }
    718846}
    … …  
    736864            polygon,
    737865            newpolygon);
    738     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     866//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     867    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    739868  }
    740869
    … …  
    756885            polygon,
    757886            newpolygon);
    758     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     887//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     888    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    759889  }
    760890
    … …  
    773903            polygon,
    774904            newpolygon);
    775     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     905//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     906    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    776907    // also slightly perturbed
    777908    const double amplitude = 0.05;
    … …  
    827958            polygon,
    828959            newpolygon);
    829     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     960//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     961    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    830962    // also slightly perturbed
    831963    const double amplitude = 0.05;
    … …  
    854986            polygon,
    855987            newpolygon);
    856     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     988//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     989    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    857990    // also slightly perturbed
    858991    const double amplitude = 0.05;
    … …  
    8801013            polygon,
    8811014            newpolygon);
    882     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1015//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1016    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    8831017  }
    8841018
    … …  
    9001034            polygon,
    9011035            newpolygon);
    902     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1036//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1037    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    9031038  }
    9041039
    … …  
    9171052            polygon,
    9181053            newpolygon);
    919     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1054//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1055    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    9201056    // also slightly perturbed
    9211057    const double amplitude = 0.05;
    … …  
    9711107            polygon,
    9721108            newpolygon);
    973     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1109//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1110    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    9741111    // also slightly perturbed
    9751112    const double amplitude = 0.05;
    … …  
    9981135            polygon,
    9991136            newpolygon);
    1000     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1137//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1138    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    10011139    // also slightly perturbed
    10021140    const double amplitude = 0.05;
    … …  
    10241162            polygon,
    10251163            newpolygon);
    1026     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1164//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1165    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    10271166  }
    10281167
    … …  
    10441183            polygon,
    10451184            newpolygon);
    1046     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1185//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1186    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    10471187  }
    10481188
    … …  
    10611201            polygon,
    10621202            newpolygon);
    1063     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1203//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1204    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    10641205    // also slightly perturbed
    10651206    const double amplitude = 0.05;
    … …  
    11151256            polygon,
    11161257            newpolygon);
    1117     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1258//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1259    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    11181260    // also slightly perturbed
    11191261    const double amplitude = 0.05;
    … …  
    11421284            polygon,
    11431285            newpolygon);
    1144     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1286//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1287    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    11451288    // also slightly perturbed
    11461289    const double amplitude = 0.05;
    … …  
    11681311            polygon,
    11691312            newpolygon);
    1170     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1313//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1314    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    11711315  }
    11721316
    … …  
    11881332            polygon,
    11891333            newpolygon);
    1190     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1334//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1335    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    11911336  }
    11921337
    … …  
    12051350            polygon,
    12061351            newpolygon);
    1207     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1352//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1353    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    12081354    // also slightly perturbed
    12091355    const double amplitude = 0.05;
    … …  
    12631409            polygon,
    12641410            newpolygon);
    1265     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1411//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1412    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    12661413    // also slightly perturbed
    12671414    const double amplitude = 0.05;
    … …  
    12901437            polygon,
    12911438            newpolygon);
    1292     CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1439//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     1440    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    12931441    // also slightly perturbed
    12941442    const double amplitude = 0.05;

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

    r1cde4e8 ra2f8a9  
    2222{
    2323    CPPUNIT_TEST_SUITE( SphericalPointDistributionTest) ;
     24    CPPUNIT_TEST( calculateCenterOfMinimumDistanceTest );
    2425    CPPUNIT_TEST ( areEqualToWithinBoundsTest );
    2526    CPPUNIT_TEST ( joinPointsTest );
    … …  
    3738      void setUp();
    3839      void tearDown();
     40      void calculateCenterOfMinimumDistanceTest();
    3941      void areEqualToWithinBoundsTest();
    4042      void joinPointsTest();
Note: See TracChangeset for help on using the changeset viewer.