/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2014 Frederik Heber. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * SphericalPointDistributionUnitTest.cpp
 *
 *  Created on: May 29, 2014
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

using namespace std;

#include <cppunit/CompilerOutputter.h>
#include <cppunit/extensions/TestFactoryRegistry.h>
#include <cppunit/ui/text/TestRunner.h>

// include headers that implement a archive in simple text format
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

#include "SphericalPointDistributionUnitTest.hpp"

#include <algorithm>
#include <boost/assign.hpp>
#include <boost/bind.hpp>
#include <numeric>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"

#include "LinearAlgebra/Line.hpp"

#include "Atom/TesselPoint.hpp"
#include "Fragmentation/Exporters/SphericalPointDistribution.hpp"
#include "LinkedCell/linkedcell.hpp"
#include "LinkedCell/PointCloudAdaptor.hpp"
#include "Tesselation/BoundaryLineSet.hpp"
#include "Tesselation/tesselation.hpp"

#ifdef HAVE_TESTRUNNER
#include "UnitTestMain.hpp"
#endif /*HAVE_TESTRUNNER*/

using namespace boost::assign;

/********************************************** Test classes **************************************/

// Registers the fixture into the 'registry'
CPPUNIT_TEST_SUITE_REGISTRATION( SphericalPointDistributionTest );

/** due to root-taking in function we only have limited numerical precision,
 * basically half of the double range.
 */
const double CenterAccuracy = sqrt(std::numeric_limits<double>::epsilon()*1e2);

void SphericalPointDistributionTest::setUp()
{
  // failing asserts should be thrown
  ASSERT_DO(Assert::Throw);

  setVerbosity(2);
}


void SphericalPointDistributionTest::tearDown()
{
}

/** UnitTest for calculateCenterOfMinimumDistance()
 */
void SphericalPointDistributionTest::calculateCenterOfMinimumDistanceTest()
{
  // single point
  {
    SphericalPointDistribution::VectorArray_t points;
    points +=
        Vector(1.,0.,0.);
    SphericalPointDistribution::IndexList_t indices;
    indices += 0;
    const Vector expected = points[0];
    const Vector center =
        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
//    CPPUNIT_ASSERT_EQUAL ( expected, center );
    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
  }

  // single point, rotated
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::VectorArray_t points;
    points +=
        RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI);
    SphericalPointDistribution::IndexList_t indices;
    indices += 0;
    const Vector expected = points[0];
    const Vector center =
        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
//    CPPUNIT_ASSERT_EQUAL ( expected, center );
    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
  }

  // two points
  {
    SphericalPointDistribution::VectorArray_t points;
    points +=
        Vector(1.,0.,0.),
        Vector(0.,1.,0.);
    SphericalPointDistribution::IndexList_t indices;
    indices += 0,1;
    const Vector expected = Vector(M_SQRT1_2,M_SQRT1_2,0.);
    const Vector center =
        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
//    CPPUNIT_ASSERT_EQUAL ( expected, center );
    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
  }

  // two points, rotated
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::VectorArray_t points;
    points +=
        RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI),
        RotationAxis.rotateVector(Vector(0.,1.,0.), 47.6/180*M_PI);
    SphericalPointDistribution::IndexList_t indices;
    indices += 0,1;
    const Vector expected = RotationAxis.rotateVector(Vector(M_SQRT1_2,M_SQRT1_2,0.), 47.6/180*M_PI);
    const Vector center =
        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
//    CPPUNIT_ASSERT_EQUAL ( expected, center );
    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
  }

  // three points in line
  {
    SphericalPointDistribution::VectorArray_t points;
    points +=
        Vector(1.,0.,0.),
        Vector(0.,1.,0.),
        Vector(-1.,0.,0.);
    SphericalPointDistribution::IndexList_t indices;
    indices += 0,1,2;
    const Vector expected = points[1];
    const Vector center =
        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
//    CPPUNIT_ASSERT_EQUAL ( expected, center );
    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
  }

  // three points in line, rotated
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::VectorArray_t points;
    points +=
        RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI),
        RotationAxis.rotateVector(Vector(0.,1.,0.), 47.6/180*M_PI),
        RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI);
    SphericalPointDistribution::IndexList_t indices;
    indices += 0,1,2;
    const Vector expected = points[1];
    const Vector center =
        SphericalPointDistribution::calculateCenterOfMinimumDistance(points, indices);
//    std::cout << " Difference is " << (expected - center).Norm() << std::endl;
//    CPPUNIT_ASSERT_EQUAL ( expected, center );
    CPPUNIT_ASSERT( expected.IsEqualTo(center, CenterAccuracy));
  }
}

static
bool areEqualToWithinBounds(
    const SphericalPointDistribution::Polygon_t &_polygon,
    const SphericalPointDistribution::Polygon_t &_otherpolygon,
    double _amplitude
    )
{
  // same size?
  if (_polygon.size() != _otherpolygon.size())
    return false;
  // same points ? We just check witrh trivial mapping, nothing fancy ...
  bool status = true;
  SphericalPointDistribution::Polygon_t::const_iterator iter = _polygon.begin();
  SphericalPointDistribution::Polygon_t::const_iterator otheriter = _otherpolygon.begin();
  for (; iter != _polygon.end(); ++iter, ++otheriter) {
    status &= (*iter).IsEqualTo(*otheriter, _amplitude);
  }
  return status;
}

/** UnitTest for areEqualToWithinBounds()
 */
void SphericalPointDistributionTest::areEqualToWithinBoundsTest()
{
  // test with no points
  {
    SphericalPointDistribution::Polygon_t polygon;
    SphericalPointDistribution::Polygon_t expected = polygon;
    CPPUNIT_ASSERT( areEqualToWithinBounds(polygon, expected, std::numeric_limits<double>::epsilon()*1e2) );
  }
  // test with one point
  {
    SphericalPointDistribution::Polygon_t polygon;
    polygon += Vector(1.,0.,0.);
    SphericalPointDistribution::Polygon_t expected = polygon;
    CPPUNIT_ASSERT( areEqualToWithinBounds(polygon, expected, std::numeric_limits<double>::epsilon()*1e2) );
  }
  // test with two points
  {
    SphericalPointDistribution::Polygon_t polygon;
    polygon += Vector(1.,0.,0.);
    polygon += Vector(0.,1.,0.);
    SphericalPointDistribution::Polygon_t expected = polygon;
    CPPUNIT_ASSERT( areEqualToWithinBounds(polygon, expected, std::numeric_limits<double>::epsilon()*1e2) );
  }

  // test with two points in different order: THIS GOES WRONG: We only check trivially
  {
    SphericalPointDistribution::Polygon_t polygon;
    polygon += Vector(1.,0.,0.);
    polygon += Vector(0.,1.,0.);
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(0.,1.,0.);
    expected += Vector(1.,0.,0.);
    CPPUNIT_ASSERT( !areEqualToWithinBounds(polygon, expected, std::numeric_limits<double>::epsilon()*1e2) );
  }

  // test with two different points
  {
    SphericalPointDistribution::Polygon_t polygon;
    polygon += Vector(1.,0.,0.);
    polygon += Vector(0.,1.,0.);
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(1.01,0.,0.);
    expected += Vector(0.,1.,0.);
    CPPUNIT_ASSERT( areEqualToWithinBounds(polygon, expected, 0.05) );
    CPPUNIT_ASSERT( !areEqualToWithinBounds(polygon, expected, 0.005) );
  }

  // test with different number of points
  {
    SphericalPointDistribution::Polygon_t polygon;
    polygon += Vector(1.,0.,0.);
    polygon += Vector(0.,1.,0.);
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(0.,1.,0.);
    CPPUNIT_ASSERT( !areEqualToWithinBounds(polygon, expected, 0.05) );
  }
}

/** The getConnectionTest test functions are templated. We only implement
 * special cases here where we need to override the behavior found in there.
 */

/** UnitTest for getConnections()
 */
template <>
void SphericalPointDistributionTest_assistant::getConnectionTest<0>()
{
  const int N=0;
  SphericalPointDistribution SPD(1.);

  // create empty adjacency
  SphericalPointDistribution::adjacency_t adjacency;

  // get the implemented connections
  SphericalPointDistribution::adjacency_t expected =
      SPD.getConnections<N>();

  // and compare the two
  CPPUNIT_ASSERT_EQUAL( expected, adjacency );
}

/** UnitTest for getConnections()
 */
template <>
void SphericalPointDistributionTest_assistant::getConnectionTest<1>()
{
  const int N=1;
  SphericalPointDistribution SPD(1.);

  // create empty adjacency
  SphericalPointDistribution::adjacency_t adjacency;

  // get the implemented connections
  SphericalPointDistribution::adjacency_t expected =
      SPD.getConnections<N>();

  // and compare the two
  CPPUNIT_ASSERT_EQUAL( expected, adjacency );
}

/** UnitTest for getConnections()
 */
template <>
void SphericalPointDistributionTest_assistant::getConnectionTest<2>()
{
  const int N=2;
  SphericalPointDistribution SPD(1.);

  // create empty adjacency
  SphericalPointDistribution::adjacency_t adjacency;
  adjacency +=
      make_pair<
      unsigned int,
      SphericalPointDistribution::IndexSet_t >
        (0, list_of<unsigned int>(1));
  adjacency +=
      make_pair<
      unsigned int,
      SphericalPointDistribution::IndexSet_t >
        (1, list_of<unsigned int>(0));

  // get the implemented connections
  SphericalPointDistribution::adjacency_t expected =
      SPD.getConnections<N>();

  // and compare the two
  CPPUNIT_ASSERT_EQUAL( expected, adjacency );
}

void freeTesselPointSTLList(TesselPointSTLList &_Corners)
{
  for (TesselPointSTLList::iterator iter = _Corners.begin();
      !_Corners.empty(); iter = _Corners.begin()) {
    delete *iter;
    _Corners.erase(iter);
  }
}

template <int N> SphericalPointDistribution::adjacency_t getAdjacencyConnections()
{
  SphericalPointDistribution SPD(1.);
  // get the points and convert into TesselPoint list
  SphericalPointDistribution::Polygon_t newpolygon = SPD.get<N>();
  TesselPointSTLList Corners;
  SphericalPointDistribution::IndexList_t indices(N);
  std::generate(indices.begin(), indices.end(), UniqueNumber);
  std::transform(
      newpolygon.begin(), newpolygon.end(),
      indices.begin(),
      std::back_inserter(Corners),
      VectorToTesselPoint());

  // create the tesselation
  const double SPHERERADIUS = 1.5;
  Tesselation TesselStruct;
  PointCloudAdaptor<TesselPointSTLList> cloud(&Corners, "TesselPointSTLList");
  TesselStruct(cloud, SPHERERADIUS);

  // create a adjacency list from a tesselation of the (convex set of) points
  SphericalPointDistribution::adjacency_t adjacency;
  for (LineMap::const_iterator iter = TesselStruct.LinesOnBoundary.begin();
      iter != TesselStruct.LinesOnBoundary.end(); ++iter) {
    const BoundaryLineSet * const line = iter->second;
    {
      std::pair< SphericalPointDistribution::adjacency_t::iterator, bool > inserter =
          adjacency.insert(
                std::make_pair(
                    line->endpoints[0]->Nr,
                    SphericalPointDistribution::IndexSet_t() ));
      inserter.first->second.insert(line->endpoints[1]->Nr);
      LOG(6, "DEBUG: Inserting " << line->endpoints[0]->Nr << "," << line->endpoints[1]->Nr);
    }
    {
      std::pair< SphericalPointDistribution::adjacency_t::iterator, bool > inserter =
          adjacency.insert(
                std::make_pair(
                    line->endpoints[1]->Nr,
                    SphericalPointDistribution::IndexSet_t() ));
      inserter.first->second.insert(line->endpoints[0]->Nr);
      LOG(6, "DEBUG: Inserting " << line->endpoints[1]->Nr << "," << line->endpoints[0]->Nr);
    }
  }

  // free allocated TesselPoints
  freeTesselPointSTLList(Corners);

  LOG(2, "INFO: adjacency is " << adjacency);

  return adjacency;
}

template <int N> SphericalPointDistribution::adjacency_t getExpectedConnections()
{
  SphericalPointDistribution SPD(1.);

  // get the implemented connections
  SphericalPointDistribution::adjacency_t expected =
      SPD.getConnections<N>();

  LOG(2, "INFO: expected is " << expected);

  return expected;
}

/** UnitTest for getConnections()
 */
template <>
void SphericalPointDistributionTest_assistant::getConnectionTest<8>()
{
  const int N=8;

  SphericalPointDistribution::adjacency_t adjacency = getAdjacencyConnections<N>();
  SphericalPointDistribution::adjacency_t expected = getExpectedConnections<N>();

  // and compare the two
//  CPPUNIT_ASSERT_EQUAL( expected, adjacency );

  // with eight points we obtain a cube. The problem is that each side
  // is a polygon with four corners that is ambiguous in the tesselation
  // it receives. Hence, we cannot directly test the linking but only
  // the properties.
  size_t NumberEdges_expected = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = expected.begin();
      iter != expected.begin(); ++iter) {
    NumberEdges_expected += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 3 );
    CPPUNIT_ASSERT( iter->second.size() <= 6 );
  }
  size_t NumberEdges_adjacency = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = adjacency.begin();
      iter != adjacency.begin(); ++iter) {
    NumberEdges_adjacency += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 3 );
    CPPUNIT_ASSERT( iter->second.size() <= 6 );
  }
  CPPUNIT_ASSERT_EQUAL( NumberEdges_expected, NumberEdges_adjacency);
}

/** UnitTest for getConnections()
 */
template <>
void SphericalPointDistributionTest_assistant::getConnectionTest<9>()
{
  const int N=9;

  SphericalPointDistribution::adjacency_t adjacency = getAdjacencyConnections<N>();
  SphericalPointDistribution::adjacency_t expected = getExpectedConnections<N>();

  // and compare the two
//  CPPUNIT_ASSERT_EQUAL( expected, adjacency );

  // with nine points we have a square on one end and an pentagon on the hand with
  // some ambiguity. The problem is that each side
  // is a polygon with four/five corners that is ambiguous in the tesselation
  // it receives. Hence, we cannot directly test the linking but only
  // the properties.
  size_t NumberEdges_expected = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = expected.begin();
      iter != expected.begin(); ++iter) {
    NumberEdges_expected += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 4 );
    CPPUNIT_ASSERT( iter->second.size() <= 5 );
  }
  size_t NumberEdges_adjacency = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = adjacency.begin();
      iter != adjacency.begin(); ++iter) {
    NumberEdges_adjacency += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 4 );
    CPPUNIT_ASSERT( iter->second.size() <= 5 );
  }
  CPPUNIT_ASSERT_EQUAL( NumberEdges_expected, NumberEdges_adjacency);
}

/** UnitTest for getConnections()
 */
template <>
void SphericalPointDistributionTest_assistant::getConnectionTest<10>()
{
  const int N=10;

  SphericalPointDistribution::adjacency_t adjacency = getAdjacencyConnections<N>();
  SphericalPointDistribution::adjacency_t expected = getExpectedConnections<N>();

  // and compare the two
//  CPPUNIT_ASSERT_EQUAL( expected, adjacency );

  // with ten points we have two pentagons with some ambiguity. The problem is
  // that each side is a polygon with five corners that is ambiguous in the
  // tesselation it receives. Hence, we cannot directly test the linking but only
  // the properties.
  size_t NumberEdges_expected = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = expected.begin();
      iter != expected.begin(); ++iter) {
    NumberEdges_expected += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 4 );
    CPPUNIT_ASSERT( iter->second.size() <= 5 );
  }
  size_t NumberEdges_adjacency = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = adjacency.begin();
      iter != adjacency.begin(); ++iter) {
    NumberEdges_adjacency += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 4 );
    CPPUNIT_ASSERT( iter->second.size() <= 5 );
  }
  CPPUNIT_ASSERT_EQUAL( NumberEdges_expected, NumberEdges_adjacency);
}

/** UnitTest for getConnections()
 */
template <>
void SphericalPointDistributionTest_assistant::getConnectionTest<11>()
{
  const int N=11;

  SphericalPointDistribution::adjacency_t adjacency = getAdjacencyConnections<N>();
  SphericalPointDistribution::adjacency_t expected = getExpectedConnections<N>();

  // and compare the two
//  CPPUNIT_ASSERT_EQUAL( expected, adjacency );

  // again, we only check properties as tesselation has ambiguities.
  size_t NumberEdges_expected = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = expected.begin();
      iter != expected.begin(); ++iter) {
    NumberEdges_expected += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 4 );
    CPPUNIT_ASSERT( iter->second.size() <= 6 );
  }
  size_t NumberEdges_adjacency = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = adjacency.begin();
      iter != adjacency.begin(); ++iter) {
    NumberEdges_adjacency += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 4 );
    CPPUNIT_ASSERT( iter->second.size() <= 6 );
  }
  CPPUNIT_ASSERT_EQUAL( NumberEdges_expected, NumberEdges_adjacency);
}

/** UnitTest for getConnections()
 */
template <>
void SphericalPointDistributionTest_assistant::getConnectionTest<12>()
{
  const int N=12;

  SphericalPointDistribution::adjacency_t adjacency = getAdjacencyConnections<N>();
  SphericalPointDistribution::adjacency_t expected = getExpectedConnections<N>();

  // and compare the two
//  CPPUNIT_ASSERT_EQUAL( expected, adjacency );

  // again, we only check properties as tesselation has ambiguities.
  size_t NumberEdges_expected = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = expected.begin();
      iter != expected.begin(); ++iter) {
    NumberEdges_expected += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 5 );
    CPPUNIT_ASSERT( iter->second.size() <= 5 );
  }
  size_t NumberEdges_adjacency = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = adjacency.begin();
      iter != adjacency.begin(); ++iter) {
    NumberEdges_adjacency += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 5 );
    CPPUNIT_ASSERT( iter->second.size() <= 5 );
  }
  CPPUNIT_ASSERT_EQUAL( NumberEdges_expected, NumberEdges_adjacency);
}

/** UnitTest for getConnections()
 */
template <>
void SphericalPointDistributionTest_assistant::getConnectionTest<14>()
{
  const int N=14;

  SphericalPointDistribution::adjacency_t adjacency = getAdjacencyConnections<N>();
  SphericalPointDistribution::adjacency_t expected = getExpectedConnections<N>();

  // and compare the two
//  CPPUNIT_ASSERT_EQUAL( expected, adjacency );

  // again, we only check properties as tesselation has ambiguities.
  size_t NumberEdges_expected = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = expected.begin();
      iter != expected.begin(); ++iter) {
    NumberEdges_expected += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 5 );
    CPPUNIT_ASSERT( iter->second.size() <= 6 );
  }
  size_t NumberEdges_adjacency = 0;
  for (SphericalPointDistribution::adjacency_t::const_iterator iter = adjacency.begin();
      iter != adjacency.begin(); ++iter) {
    NumberEdges_adjacency += iter->second.size();
    CPPUNIT_ASSERT( iter->second.size() >= 5 );
    CPPUNIT_ASSERT( iter->second.size() <= 6 );
  }
  CPPUNIT_ASSERT_EQUAL( NumberEdges_expected, NumberEdges_adjacency);
}

void perturbPolygon(
    SphericalPointDistribution::WeightedPolygon_t &_polygon,
    double _amplitude
    )
{
  for (SphericalPointDistribution::WeightedPolygon_t::iterator iter = _polygon.begin();
      iter != _polygon.end(); ++iter) {
    Vector perturber;
    perturber.GetOneNormalVector(iter->first);
    perturber.Scale(_amplitude);
    iter->first = iter->first + perturber;
    (iter->first).Normalize();
  }
}

/** UnitTest for joinPoints()
 */
void SphericalPointDistributionTest::joinPointsTest()
{
  // test with simple configuration of three points
  {
    SphericalPointDistribution::Polygon_t newpolygon;
    newpolygon += Vector(1.,0.,0.);
    newpolygon += Vector(0.,1.,0.);
    newpolygon += Vector(0.,0.,1.);
    SphericalPointDistribution::Polygon_t expectedpolygon = newpolygon;
    SphericalPointDistribution::IndexTupleList_t matching;
    matching += SphericalPointDistribution::IndexList_t(1,0);
    matching += SphericalPointDistribution::IndexList_t(1,1);
    matching += SphericalPointDistribution::IndexList_t(1,2);
    SphericalPointDistribution::IndexList_t IndexList =
        SphericalPointDistribution::joinPoints(
            newpolygon,
            SphericalPointDistribution::VectorArray_t(newpolygon.begin(), newpolygon.end()),
            matching);
    SphericalPointDistribution::IndexList_t expected;
    expected += 0,1,2;
    CPPUNIT_ASSERT_EQUAL( expected, IndexList );
    CPPUNIT_ASSERT_EQUAL( expectedpolygon, newpolygon );
  }

  // test with simple configuration of three points, only two are picked
  {
    SphericalPointDistribution::Polygon_t newpolygon;
    newpolygon += Vector(1.,0.,0.);
    newpolygon += Vector(0.,1.,0.);
    newpolygon += Vector(0.,0.,1.);
    SphericalPointDistribution::Polygon_t expectedpolygon = newpolygon;
    SphericalPointDistribution::IndexTupleList_t matching;
    matching += SphericalPointDistribution::IndexList_t(1,1);
    matching += SphericalPointDistribution::IndexList_t(1,2);
    SphericalPointDistribution::IndexList_t IndexList =
        SphericalPointDistribution::joinPoints(
            newpolygon,
            SphericalPointDistribution::VectorArray_t(newpolygon.begin(), newpolygon.end()),
            matching);
    SphericalPointDistribution::IndexList_t expected;
    expected += 1,2;
    CPPUNIT_ASSERT_EQUAL( expected, IndexList );
    CPPUNIT_ASSERT_EQUAL( expectedpolygon, newpolygon );
  }

  // test with simple configuration of three points, two are joined
  {
    SphericalPointDistribution::Polygon_t newpolygon;
    newpolygon += Vector(1.,0.,0.);
    newpolygon += Vector(0.,1.,0.);
    newpolygon += Vector(0.,0.,1.);
    SphericalPointDistribution::Polygon_t expectedpolygon;
    expectedpolygon += Vector(1.,0.,0.);
    expectedpolygon += Vector(0.,M_SQRT1_2,M_SQRT1_2);
    SphericalPointDistribution::IndexTupleList_t matching;
    SphericalPointDistribution::IndexList_t joined;
    joined += 1,2;
    matching += SphericalPointDistribution::IndexList_t(1,0);
    matching += joined;
    SphericalPointDistribution::IndexList_t IndexList =
        SphericalPointDistribution::joinPoints(
            newpolygon,
            SphericalPointDistribution::VectorArray_t(newpolygon.begin(), newpolygon.end()),
            matching);
    SphericalPointDistribution::IndexList_t expected;
    expected += 0,1;
    CPPUNIT_ASSERT_EQUAL( expected, IndexList );
    CPPUNIT_ASSERT_EQUAL( expectedpolygon, newpolygon );
  }

  // test with simple configuration of six points, two are joined, jumbled indices
  {
    SphericalPointDistribution::Polygon_t newpolygon;
    newpolygon += Vector(1.,0.,1.);
    newpolygon += Vector(1.,0.,0.);
    newpolygon += Vector(1.,1.,0.);
    newpolygon += Vector(0.,1.,0.);
    newpolygon += Vector(0.,0.,1.);
    newpolygon += Vector(1.,0.,1.);
    SphericalPointDistribution::Polygon_t expectedpolygon;
    expectedpolygon += Vector(1.,0.,1.);
    expectedpolygon += Vector(1.,0.,0.);
    expectedpolygon += Vector(1.,1.,0.);
    expectedpolygon += Vector(1.,0.,1.);
    expectedpolygon += Vector(0.,M_SQRT1_2,M_SQRT1_2); // new centers go last
    SphericalPointDistribution::IndexTupleList_t matching;
    SphericalPointDistribution::IndexList_t joined;
    joined += 3,4;
    matching += SphericalPointDistribution::IndexList_t(1,1);
    matching += joined;
    SphericalPointDistribution::IndexList_t IndexList =
        SphericalPointDistribution::joinPoints(
            newpolygon,
            SphericalPointDistribution::VectorArray_t(newpolygon.begin(), newpolygon.end()),
            matching);
    SphericalPointDistribution::IndexList_t expected;
    expected += 1,4;
    CPPUNIT_ASSERT_EQUAL( expected, IndexList );
    CPPUNIT_ASSERT_EQUAL( expectedpolygon, newpolygon );
  }
}

/** UnitTest for getRemainingPoints() with two points
 */
void SphericalPointDistributionTest::getRemainingPointsTest_2()
{
  SphericalPointDistribution SPD(1.);
  // test with one point, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(Vector(1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(-1.,0.,0.);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 2);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with one point, just a flip of axis
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(0.,1.,0.), 1);
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(0.,-1.,0.);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 2);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with one point, just a flip to another axis
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(0.,0.,-1.), 1);
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(0.,0.,1.);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 2);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with one point, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected;
    expected += RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 2);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }
}

/** UnitTest for getRemainingPoints() with three points
 */
void SphericalPointDistributionTest::getRemainingPointsTest_3()
{
  SphericalPointDistribution SPD(1.);

  // test with one point, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<3>();
    expected.pop_front(); // remove first point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 3);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with one point, just a flip of x and y axis
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(0.,1.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<3>();
    expected.pop_front(); // remove first point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter) {
      std::swap((*iter)[0], (*iter)[1]);
      (*iter)[0] *= -1.;
    }
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 3);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with two points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-0.5, sqrt(3)*0.5,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<3>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 3);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with two points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-0.5, sqrt(3)*0.5,0.), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<3>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 3);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with three points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-0.5, sqrt(3)*0.5,0.), 1);
    polygon += std::make_pair( Vector(-0.5, -sqrt(3)*0.5,0.), 1);
    SphericalPointDistribution::Polygon_t newpolygon =
        SPD.get<3>();
    SphericalPointDistribution::Polygon_t expected; // empty cause none are vacant
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 3);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }


  // test with three points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-0.5, sqrt(3)*0.5,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-0.5, -sqrt(3)*0.5,0.), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t newpolygon =
        SPD.get<3>();
    SphericalPointDistribution::Polygon_t expected; // empty cause none are vacant
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 3);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }
}

/** UnitTest for getRemainingPoints() with four points
 */
void SphericalPointDistributionTest::getRemainingPointsTest_4()
{
  SphericalPointDistribution SPD(1.);

  // test with one point, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<4>();
    expected.pop_front(); // remove first point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 4);
    //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
        CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with one point, just a flip of axis
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(0.,1.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<4>();
    expected.pop_front(); // remove first point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter) {
      std::swap((*iter)[0], (*iter)[1]);
      (*iter)[0] *= -1.;
    }
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 4);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with two points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<4>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 4);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with two points, matching trivially, also with slightly perturbed
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<4>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 4);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with two points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<4>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 4);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with three points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 1);
    polygon += std::make_pair( Vector(-1./3.0, -M_SQRT2/3.0, M_SQRT2/sqrt(3)), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<4>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 4);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with three points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1./3.0, -M_SQRT2/3.0, M_SQRT2/sqrt(3)), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<4>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 4);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }
}

/** UnitTest for getRemainingPoints() with four points and weights
 * not all equal to one.
 */
void SphericalPointDistributionTest::getRemainingPointsTest_multiple()
{
  SphericalPointDistribution SPD(1.);

  // test with four points: one point having weight of two
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 2);
    SphericalPointDistribution::Polygon_t newpolygon =
        SPD.get<4>();
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(-0.5773502691896,-5.551115123126e-17,0.8164965809277);
    expected += Vector(-0.5773502691896,-5.551115123126e-17,-0.8164965809277);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 4);
//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with five points: one point having weight of two
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 2);
    SphericalPointDistribution::Polygon_t newpolygon =
        SPD.get<5>();
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(-0.7071067811865,0.7071067811865,0);
    expected += Vector(-0.3535533905933,-0.3535533905933,0.8660254037844);
    expected += Vector(-0.3535533905933,-0.3535533905933,-0.8660254037844);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 5);
    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }


  // test with five points: one point having weight of two, one weight of one
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(M_SQRT1_2,M_SQRT1_2,0.), 2);
    polygon += std::make_pair( Vector(-1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t newpolygon =
        SPD.get<5>();
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(0.3535533786708,-0.3535533955317,-0.8660254066357);
    expected += Vector(0.3535534025157,-0.3535533856548,0.8660254009332);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 5);
//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with six points: two points each having weight of two
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(M_SQRT1_2,-M_SQRT1_2,0.), 2);
    polygon += std::make_pair( Vector(-M_SQRT1_2,M_SQRT1_2,0.), 2);
    SphericalPointDistribution::Polygon_t newpolygon =
        SPD.get<6>();
    SphericalPointDistribution::Polygon_t expected;
    expected += Vector(0.,0.,1.);
    expected += Vector(0.,0.,-1.);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 6);
//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }
}

/** UnitTest for getRemainingPoints() with five points
 */
void SphericalPointDistributionTest::getRemainingPointsTest_5()
{
  SphericalPointDistribution SPD(1.);

  // test with one point, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<5>();
    expected.pop_front(); // remove first point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 5);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with one point, just a flip of axis
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(0.,1.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<5>();
    expected.pop_front(); // remove first point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter) {
      std::swap((*iter)[0], (*iter)[1]);
      (*iter)[0] *= -1.;
    }
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 5);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with two points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<5>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 5);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with two points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180.*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180.*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<5>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180.*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 5);
    // the three remaining points sit on a plane that may be rotated arbitrarily
    // so we cannot simply check for equality between expected and remaining
    // hence, we just check that they are orthogonal to the first two points
    CPPUNIT_ASSERT_EQUAL( expected.size(), remaining.size() );
    for (SphericalPointDistribution::WeightedPolygon_t::const_iterator fixiter = polygon.begin();
        fixiter != polygon.end(); ++fixiter) {
      for (SphericalPointDistribution::Polygon_t::const_iterator iter = remaining.begin();
          iter != remaining.end(); ++iter) {
        CPPUNIT_ASSERT( (fixiter->first).IsNormalTo(*iter) );
      }
    }
  }

  // test with three points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1., 0.0, 0.0), 1);
    polygon += std::make_pair( Vector(0.0, 1., 0.0), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<5>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 5);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with three points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1., 0.0, 0.0), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(0.0, 1., 0.0), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<5>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 5);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }
}

/** UnitTest for getRemainingPoints() with six points
 */
void SphericalPointDistributionTest::getRemainingPointsTest_6()
{
  SphericalPointDistribution SPD(1.);

  // test with one point, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<6>();
    expected.pop_front(); // remove first point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 6);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with one point, just a flip of axis
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(0.,1.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<6>();
    expected.pop_front(); // remove first point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter) {
      std::swap((*iter)[0], (*iter)[1]);
      (*iter)[0] *= -1.;
    }
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 6);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with two points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<6>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second spoint
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 6);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with two points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<6>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second spoint
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 6);
    // the four remaining points sit on a plane that may have been rotated arbitrarily
    // so we cannot simply check for equality between expected and remaining
    // hence, we just check that they are orthogonal to the first two points
    CPPUNIT_ASSERT_EQUAL( expected.size(), remaining.size() );
    for (SphericalPointDistribution::WeightedPolygon_t::const_iterator fixiter = polygon.begin();
        fixiter != polygon.end(); ++fixiter) {
      for (SphericalPointDistribution::Polygon_t::const_iterator iter = remaining.begin();
          iter != remaining.end(); ++iter) {
        CPPUNIT_ASSERT( (fixiter->first).IsNormalTo(*iter) );
      }
    }
  }

  // test with three points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1., 0.0, 0.0), 1);
    polygon += std::make_pair( Vector(0.0, 1., 0.0), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<6>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 6);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with three points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1., 0.0, 0.0), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(0.0, 1., 0.0), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<6>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 6);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }
}

/** UnitTest for getRemainingPoints() with seven points
 */
void SphericalPointDistributionTest::getRemainingPointsTest_7()
{
  SphericalPointDistribution SPD(1.);

  // test with one point, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<7>();
    expected.pop_front(); // remove first point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 7);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with one point, just a flip of axis
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(0.,1.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<7>();
    expected.pop_front(); // remove first point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter) {
      std::swap((*iter)[0], (*iter)[1]);
      (*iter)[0] *= -1.;
    }
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 7);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with two points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<7>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 7);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with two points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<7>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 7);
    // the five remaining points sit on a plane that may have been rotated arbitrarily
    // so we cannot simply check for equality between expected and remaining
    // hence, we just check that they are orthogonal to the first two points
    CPPUNIT_ASSERT_EQUAL( expected.size(), remaining.size() );
    for (SphericalPointDistribution::WeightedPolygon_t::const_iterator fixiter = polygon.begin();
        fixiter != polygon.end(); ++fixiter) {
      for (SphericalPointDistribution::Polygon_t::const_iterator iter = remaining.begin();
          iter != remaining.end(); ++iter) {
        CPPUNIT_ASSERT( (fixiter->first).IsNormalTo(*iter) );
      }
    }
  }

  // test with three points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1., 0.0, 0.0), 1);
    polygon += std::make_pair( Vector(0.0, 1., 0.0), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<7>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 7);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with three points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1., 0.0, 0.0), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(0.0, 1., 0.0), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<7>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 7);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }
}

/** UnitTest for getRemainingPoints() with eight points
 */
void SphericalPointDistributionTest::getRemainingPointsTest_8()
{
  SphericalPointDistribution SPD(1.);

  // test with one point, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<8>();
    expected.pop_front(); // remove first point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 8);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with one point, just a flip of axis
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(0.,1.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<8>();
    expected.pop_front(); // remove first point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter) {
      std::swap((*iter)[0], (*iter)[1]);
      (*iter)[0] *= -1.;
    }
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 8);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
  }

  // test with two points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1.,0.,0.), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<8>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 8);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with two points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<8>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 8);
    // the six remaining points sit on two planes that may have been rotated arbitrarily
    // so we cannot simply check for equality between expected and remaining
    // hence, we just check that they are orthogonal to the first two points
    CPPUNIT_ASSERT_EQUAL( expected.size(), remaining.size() );
    for (SphericalPointDistribution::WeightedPolygon_t::const_iterator fixiter = polygon.begin();
        fixiter != polygon.end(); ++fixiter) {
      SphericalPointDistribution::Polygon_t::const_iterator expectiter = expected.begin();
      SphericalPointDistribution::Polygon_t::const_iterator remainiter = remaining.begin();
      for (;remainiter != remaining.end(); ++expectiter, ++remainiter) {
        // check that points in expected/remaining have same angle to the given ones
//        CPPUNIT_ASSERT_EQUAL( (*expectiter).Angle(*fixiter), (*remainiter).Angle(*fixiter) );
        CPPUNIT_ASSERT( fabs( (*expectiter).Angle(fixiter->first) - (*remainiter).Angle(fixiter->first) )
            < std::numeric_limits<double>::epsilon()*1e4 );
      }
    }
  }

  // test with three points, matching trivially
  {
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair( Vector(1.,0.,0.), 1);
    polygon += std::make_pair( Vector(-1., 0.0, 0.0), 1);
    polygon += std::make_pair( Vector(-1./3.0, 2.0*M_SQRT2/3.0, 0.0), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<8>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 8);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }

  // test with three points, full rotation
  {
    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
    SphericalPointDistribution::WeightedPolygon_t polygon;
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1., 0.0, 0.0), 47.6/180*M_PI), 1);
    polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1./3.0, 2.0*M_SQRT2/3.0, 0.0), 47.6/180*M_PI), 1);
    SphericalPointDistribution::Polygon_t expected =
        SPD.get<8>();
    expected.pop_front(); // remove first point
    expected.pop_front(); // remove second point
    expected.pop_front(); // remove third point
    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
        iter != expected.end(); ++iter)
      *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
    SphericalPointDistribution::Polygon_t remaining =
        SPD.getRemainingPoints(polygon, 8);
//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
    // also slightly perturbed
    const double amplitude = 0.05;
    perturbPolygon(polygon, amplitude);
    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
  }
}