source: ThirdParty/LinearAlgebra/src/unittests/PlaneUnitTest.cpp@ 6a5921

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. 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/>.
 */

/*
 * PlaneUnittest.cpp
 *
 *  Created on: Apr 30, 2010
 *      Author: crueger
 */

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

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

#include <cmath>
#include <limits>

#include "defs.hpp"
#include "Exceptions.hpp"
#include "Line.hpp"
#include "Vector.hpp"

#include "PlaneUnitTest.hpp"

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

CPPUNIT_TEST_SUITE_REGISTRATION( PlaneUnittest );

void PlaneUnittest::setUp(){
  p1 = new Plane(unitVec[0],unitVec[1],unitVec[2]);
  p2 = new Plane(unitVec[0],unitVec[1],zeroVec);
  p3 = new Plane(unitVec[0],zeroVec,unitVec[2]);
  p4 = new Plane(zeroVec,unitVec[1],unitVec[2]);
}

void PlaneUnittest::tearDown(){
  delete p1;
  delete p2;
  delete p3;
  delete p4;
}

void PlaneUnittest::constructionErrorTest(){
  // try several method of construction..
  // see if error checking works

  // three points
  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],unitVec[1],unitVec[2]));
  // when only two points are differnt this gives an error
  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],unitVec[1],unitVec[1]),LinearDependenceException);
  // same with only one point
  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],unitVec[0],unitVec[0]),LinearDependenceException);

  // use two vector giving two directions
  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],unitVec[1],0));
  // and again this is actually only one vector
  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],unitVec[0],0),LinearDependenceException);
  // Zero vector does not give a good direction
  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],zeroVec,0),ZeroVectorException);

  // use a normalvector and an scalar offset
  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],0));
  // The zero vector is no good as a normalvector
  CPPUNIT_ASSERT_THROW(Plane(zeroVec,0),ZeroVectorException);

  // use a normalvector and an offset vector
  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],zeroVec));
  // and the bad zeroVector again
  CPPUNIT_ASSERT_THROW(Plane(zeroVec,zeroVec),ZeroVectorException);
}


// we need to test normals independent of the direction
bool testNormal(const Vector &normal1, const Vector &normal2){
  return (normal1==normal2) || (normal1==-1*normal2);
}

void PlaneUnittest::constructionResultTest(){
  {
    // construct with three points on plane
    Plane p1(unitVec[0],unitVec[1],zeroVec);
    CPPUNIT_ASSERT(testNormal(unitVec[2],p1.getNormal()));
    CPPUNIT_ASSERT_EQUAL(0.,p1.getOffset());

    Plane p2(unitVec[0],unitVec[2],zeroVec);
    CPPUNIT_ASSERT(testNormal(unitVec[1],p2.getNormal()));
    CPPUNIT_ASSERT_EQUAL(0.,p2.getOffset());

    Plane p3(unitVec[1],unitVec[2],zeroVec);
    CPPUNIT_ASSERT(testNormal(unitVec[0],p3.getNormal()));
    CPPUNIT_ASSERT_EQUAL(0.,p3.getOffset());
  }
  {
    // construct with two directions + offset
    Plane p1(unitVec[0],unitVec[1],0);
    CPPUNIT_ASSERT(testNormal(unitVec[2],p1.getNormal()));
    CPPUNIT_ASSERT_EQUAL(0.,p1.getOffset());

    Plane p2(unitVec[0],unitVec[2],0);
    CPPUNIT_ASSERT(testNormal(unitVec[1],p2.getNormal()));
    CPPUNIT_ASSERT_EQUAL(0.,p2.getOffset());

    Plane p3(unitVec[1],unitVec[2],0);
    CPPUNIT_ASSERT(testNormal(unitVec[0],p3.getNormal()));
    CPPUNIT_ASSERT_EQUAL(0.,p3.getOffset());
  }
}

void PlaneUnittest::pointsTest(){
  std::vector<Vector> points1 = p1->getPointsOnPlane();
  CPPUNIT_ASSERT(p1->isContained(points1[0]));
  CPPUNIT_ASSERT(p1->isContained(points1[1]));
  CPPUNIT_ASSERT(p1->isContained(points1[2]));
  // check that the three points differ
  CPPUNIT_ASSERT(points1[0]!=points1[1]);
  CPPUNIT_ASSERT(points1[0]!=points1[2]);
  CPPUNIT_ASSERT(points1[1]!=points1[2]);


  std::vector<Vector> points2 = p2->getPointsOnPlane();
  CPPUNIT_ASSERT(p2->isContained(points2[0]));
  CPPUNIT_ASSERT(p2->isContained(points2[1]));
  CPPUNIT_ASSERT(p2->isContained(points2[2]));
  // check that the three points differ
  CPPUNIT_ASSERT(points2[0]!=points2[1]);
  CPPUNIT_ASSERT(points2[0]!=points2[2]);
  CPPUNIT_ASSERT(points2[1]!=points2[2]);

  std::vector<Vector> points3 = p3->getPointsOnPlane();
  CPPUNIT_ASSERT(p3->isContained(points3[0]));
  CPPUNIT_ASSERT(p3->isContained(points3[1]));
  CPPUNIT_ASSERT(p3->isContained(points3[2]));
  // check that the three points differ
  CPPUNIT_ASSERT(points3[0]!=points3[1]);
  CPPUNIT_ASSERT(points3[0]!=points3[2]);
  CPPUNIT_ASSERT(points3[1]!=points3[2]);

  std::vector<Vector> points4 = p4->getPointsOnPlane();
  CPPUNIT_ASSERT(p4->isContained(points4[0]));
  CPPUNIT_ASSERT(p4->isContained(points4[1]));
  CPPUNIT_ASSERT(p4->isContained(points4[2]));
  // check that the three points differ
  CPPUNIT_ASSERT(points4[0]!=points4[1]);
  CPPUNIT_ASSERT(points4[0]!=points4[2]);
  CPPUNIT_ASSERT(points4[1]!=points4[2]);
}


void PlaneUnittest::operationsTest(){
  {
    Vector t = (1./3.)*(unitVec[0]+unitVec[1]+unitVec[2]);
    CPPUNIT_ASSERT(fabs(p1->distance(zeroVec)-t.Norm()) <= LINALG_MYEPSILON());
    CPPUNIT_ASSERT_EQUAL(t,p1->getClosestPoint(zeroVec));
  }

  CPPUNIT_ASSERT(fabs(p2->distance(unitVec[2])-1) <= LINALG_MYEPSILON());
  CPPUNIT_ASSERT_EQUAL(zeroVec,p2->getClosestPoint(unitVec[2]));
  CPPUNIT_ASSERT(fabs(p3->distance(unitVec[1])-1) <= LINALG_MYEPSILON());
  CPPUNIT_ASSERT_EQUAL(zeroVec,p3->getClosestPoint(unitVec[1]));
  CPPUNIT_ASSERT(fabs(p4->distance(unitVec[0])-1) <= LINALG_MYEPSILON());
  CPPUNIT_ASSERT_EQUAL(zeroVec,p4->getClosestPoint(unitVec[0]));
}

void PlaneUnittest::mirrorTest(){
  Vector fixture;

  // some Vectors that lie on the planes
  fixture = p1->mirrorVector(unitVec[0]);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = p1->mirrorVector(unitVec[1]);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = p1->mirrorVector(unitVec[2]);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  fixture = p2->mirrorVector(zeroVec);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = p2->mirrorVector(unitVec[0]);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = p2->mirrorVector(unitVec[1]);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);

  fixture = p3->mirrorVector(zeroVec);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = p3->mirrorVector(unitVec[0]);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = p3->mirrorVector(unitVec[2]);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  fixture = p4->mirrorVector(zeroVec);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = p4->mirrorVector(unitVec[1]);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = p4->mirrorVector(unitVec[2]);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  // some Vectors outside of the planes
  {
    Vector t = (2./3.)*(unitVec[0]+unitVec[1]+unitVec[2]);
    fixture = p1->mirrorVector(zeroVec);
    CPPUNIT_ASSERT_EQUAL(fixture,t);
  }

  fixture = p2->mirrorVector(unitVec[2]);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
  fixture = p3->mirrorVector(unitVec[1]);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
  fixture = p4->mirrorVector(unitVec[0]);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
}

void PlaneUnittest::LineIntersectionTest(){
  Vector fixture;
  // plane at (0,0,0) normal to (1,0,0) cuts line from (0,0,0) to (2,1,0) at ???
  Line l1 = makeLineThrough(zeroVec,Vector(2,1,0));
  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(unitVec[0], zeroVec).GetIntersection(l1) );
  CPPUNIT_ASSERT_EQUAL( zeroVec, fixture );

  // plane at (2,1,0) normal to (0,1,0) cuts line from (1,0,0) to (0,1,1) at ???
  Line l2 = makeLineThrough(unitVec[0],Vector(0,1,1));
  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(unitVec[1], Vector(2,1,0)).GetIntersection(l2) );
  CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
}