/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010 University of Bonn. All rights reserved. * Please see the LICENSE file or "Copyright notice" in builder.cpp for details. */ /* * Line.cpp * * Created on: Apr 30, 2010 * Author: crueger */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "CodePatterns/MemDebug.hpp" #include "LinearAlgebra/Line.hpp" #include #include #include "LinearAlgebra/Vector.hpp" #include "CodePatterns/Log.hpp" #include "CodePatterns/Verbose.hpp" #include "LinearAlgebra/MatrixContent.hpp" #include "CodePatterns/Info.hpp" #include "Exceptions/LinearDependenceException.hpp" #include "Exceptions/SkewException.hpp" #include "LinearAlgebra/Plane.hpp" using namespace std; Line::Line(const Vector &_origin, const Vector &_direction) : direction(new Vector(_direction)) { direction->Normalize(); origin.reset(new Vector(_origin.partition(*direction).second)); } Line::Line(const Line &src) : origin(new Vector(*src.origin)), direction(new Vector(*src.direction)) {} Line::~Line() {} Line &Line::operator=(const Line& rhs){ if(this!=&rhs){ origin.reset(new Vector(*rhs.origin)); direction.reset(new Vector(*rhs.direction)); } return *this; } double Line::distance(const Vector &point) const{ // get any vector from line to point Vector helper = point - *origin; // partition this vector along direction // the residue points from the line to the point return helper.partition(*direction).second.Norm(); } Vector Line::getClosestPoint(const Vector &point) const{ // get any vector from line to point Vector helper = point - *origin; // partition this vector along direction // add only the part along the direction return *origin + helper.partition(*direction).first; } Vector Line::getDirection() const{ return *direction; } Vector Line::getOrigin() const{ return *origin; } vector Line::getPointsOnLine() const{ vector res; res.reserve(2); res.push_back(*origin); res.push_back(*origin+*direction); return res; } /** Calculates the intersection of the two lines that are both on the same plane. * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html * \param *out output stream for debugging * \param *Line1a first vector of first line * \param *Line1b second vector of first line * \param *Line2a first vector of second line * \param *Line2b second vector of second line * \return true - \a this will contain the intersection on return, false - lines are parallel */ Vector Line::getIntersection(const Line& otherLine) const{ Info FunctionInfo(__func__); pointset line1Points = getPointsOnLine(); Vector Line1a = line1Points[0]; Vector Line1b = line1Points[1]; pointset line2Points = otherLine.getPointsOnLine(); Vector Line2a = line2Points[0]; Vector Line2b = line2Points[1]; Vector res; auto_ptr M = auto_ptr(new MatrixContent(4,4)); M->setValue(1.); for (int i=0;i<3;i++) { M->set(0, i, Line1a[i]); M->set(1, i, Line1b[i]); M->set(2, i, Line2a[i]); M->set(3, i, Line2b[i]); } //Log() << Verbose(1) << "Coefficent matrix is:" << endl; //for (int i=0;i<4;i++) { // for (int j=0;j<4;j++) // cout << "\t" << M->Get(i,j); // cout << endl; //} if (fabs(M->Determinant()) > MYEPSILON) { Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl; throw SkewException(__FILE__,__LINE__); } Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl; // constuct a,b,c Vector a = Line1b - Line1a; Vector b = Line2b - Line2a; Vector c = Line2a - Line1a; Vector d = Line2b - Line1b; Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl; if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) { res.Zero(); Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl; throw LinearDependenceException(__FILE__,__LINE__); } // check for parallelity Vector parallel; double factor = 0.; if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) { parallel = Line1a - Line2a; factor = parallel.ScalarProduct(a)/a.Norm(); if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) { res = Line2a; Log() << Verbose(1) << "Lines conincide." << endl; return res; } else { parallel = Line1a - Line2b; factor = parallel.ScalarProduct(a)/a.Norm(); if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) { res = Line2b; Log() << Verbose(1) << "Lines conincide." << endl; return res; } } Log() << Verbose(1) << "Lines are parallel." << endl; res.Zero(); throw LinearDependenceException(__FILE__,__LINE__); } // obtain s double s; Vector temp1, temp2; temp1 = c; temp1.VectorProduct(b); temp2 = a; temp2.VectorProduct(b); Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl; if (fabs(temp2.NormSquared()) > MYEPSILON) s = temp1.ScalarProduct(temp2)/temp2.NormSquared(); else s = 0.; Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl; // construct intersection res = a; res.Scale(s); res += Line1a; Log() << Verbose(1) << "Intersection is at " << res << "." << endl; return res; } /** Rotates the vector by an angle of \a alpha around this line. * \param rhs Vector to rotate * \param alpha rotation angle in radian */ Vector Line::rotateVector(const Vector &rhs, double alpha) const{ Vector helper = rhs; // translate the coordinate system so that the line goes through (0,0,0) helper -= *origin; // partition the vector into a part that gets rotated and a part that lies along the line pair parts = helper.partition(*direction); // we just keep anything that is along the axis Vector res = parts.first; // the rest has to be rotated Vector a = parts.second; // we only have to do the rest, if we actually could partition the vector if(!a.IsZero()){ // construct a vector that is orthogonal to a and direction and has length |a| Vector y = a; // direction is normalized, so the result has length |a| y.VectorProduct(*direction); res += cos(alpha) * a + sin(alpha) * y; } // translate the coordinate system back res += *origin; return res; } Line Line::rotateLine(const Line &rhs, double alpha) const{ Vector lineOrigin = rotateVector(rhs.getOrigin(),alpha); Vector helper = rhs.getDirection(); // rotate the direction without considering the ofset pair parts = helper.partition(*direction); Vector lineDirection = parts.first; Vector a = parts.second; if(!a.IsZero()){ // construct a vector that is orthogonal to a and direction and has length |a| Vector y = a; // direction is normalized, so the result has length |a| y.VectorProduct(*direction); lineDirection += cos(alpha) * a + sin(alpha) * y; } return Line(lineOrigin,lineDirection); } Plane Line::rotatePlane(const Plane &rhs, double alpha) const{ vector points = rhs.getPointsOnPlane(); transform(points.begin(), points.end(), points.begin(), boost::bind(&Line::rotateVector,this,_1,alpha)); return Plane(points[0],points[1],points[2]); } Plane Line::getOrthogonalPlane(const Vector &origin) const{ return Plane(getDirection(),origin); } std::vector Line::getSphereIntersections() const{ std::vector res; // line is kept in normalized form, so we can skip a lot of calculations double discriminant = 1-origin->NormSquared(); // we might have 2, 1 or 0 solutions, depending on discriminant if(discriminant>=0){ if(discriminant==0){ res.push_back(*origin); } else{ Vector helper = sqrt(discriminant)*(*direction); res.push_back(*origin+helper); res.push_back(*origin-helper); } } return res; } LinePoint Line::getLinePoint(const Vector &point) const{ ASSERT(isContained(point),"Line point queried for point not on line"); Vector helper = point - (*origin); double param = helper.ScalarProduct(*direction); return LinePoint(*this,param); } LinePoint Line::posEndpoint() const{ return LinePoint(*this, numeric_limits::infinity()); } LinePoint Line::negEndpoint() const{ return LinePoint(*this,-numeric_limits::infinity()); } bool operator==(const Line &x,const Line &y){ return *x.origin == *y.origin && *x.direction == *y.direction; } Line makeLineThrough(const Vector &x1, const Vector &x2){ if(x1==x2){ throw LinearDependenceException(__FILE__,__LINE__); } return Line(x1,x1-x2); } /******************************** Points on the line ********************/ LinePoint::LinePoint(const LinePoint &src) : line(src.line),param(src.param) {} LinePoint::LinePoint(const Line &_line, double _param) : line(_line),param(_param) {} LinePoint& LinePoint::operator=(const LinePoint &src){ line=src.line; param=src.param; return *this; } Vector LinePoint::getPoint() const{ ASSERT(!isInfinite(),"getPoint() on infinite LinePoint called"); return (*line.origin)+param*(*line.direction); } Line LinePoint::getLine() const{ return line; } bool LinePoint::isInfinite() const{ return isPosInfinity() || isNegInfinity(); } bool LinePoint::isPosInfinity() const{ return param == numeric_limits::infinity(); } bool LinePoint::isNegInfinity() const{ return param ==-numeric_limits::infinity(); } bool operator==(const LinePoint &x, const LinePoint &y){ ASSERT(x.line==y.line,"Operation on two points of different lines"); return x.param == y.param; } bool operator<(const LinePoint &x, const LinePoint &y){ ASSERT(x.line==y.line,"Operation on two points of different lines"); return x.param ("; for (int i=0;i