source: ThirdParty/LinearAlgebra/src/LinearAlgebra/Line.cpp@ aa91de0

/*
 * 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/>.
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include "Line.hpp"

#include <cmath>
#include <iostream>
#include <limits>

#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "defs.hpp"
#include "Exceptions.hpp"
#include "LinePoint.hpp"
#include "MatrixContent.hpp"
#include "Plane.hpp"
#include "RealSpaceMatrix.hpp"
#include "Vector.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<Vector> Line::getPointsOnLine() const{
  vector<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{
  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<MatrixContent> M = auto_ptr<MatrixContent>(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()) > LINALG_MYEPSILON()) {
    eLog() << Verbose(2) << "Determinant of coefficient matrix is NOT zero." << endl;
    throw SkewException() << LinearAlgebraDeterminant(M->Determinant()) << LinearAlgebraMatrixContent(&(*M));
  }

  Log() << Verbose(6) << "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(7) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
  if ((a.NormSquared() <= LINALG_MYEPSILON()) || (b.NormSquared() <= LINALG_MYEPSILON())) {
   res.Zero();
   eLog() << Verbose(2) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
   throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&a, &b) );
  }

  // check for parallelity
  Vector parallel;
  double factor = 0.;
  if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) <= LINALG_MYEPSILON()) {
    parallel = Line1a - Line2a;
    factor = parallel.ScalarProduct(a)/a.Norm();
    if ((factor > -LINALG_MYEPSILON()) && (factor - 1. <= LINALG_MYEPSILON())) {
      res = Line2a;
      Log() << Verbose(5) << "Lines conincide." << endl;
      return res;
    } else {
      parallel = Line1a - Line2b;
      factor = parallel.ScalarProduct(a)/a.Norm();
      if ((factor > -LINALG_MYEPSILON()) && (factor - 1. <= LINALG_MYEPSILON())) {
        res = Line2b;
        Log() << Verbose(5) << "Lines conincide." << endl;
        return res;
      }
    }
    eLog() << Verbose(2) << "Lines are parallel." << endl;
    res.Zero();
    throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&a, &b) );
  }

  // obtain s
  double s;
  Vector temp1, temp2;
  temp1 = c;
  temp1.VectorProduct(b);
  temp2 = a;
  temp2.VectorProduct(b);
  Log() << Verbose(7) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
  if (fabs(temp2.NormSquared()) > LINALG_MYEPSILON())
    s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
  else
    s = 0.;
  Log() << Verbose(6) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;

  // construct intersection
  res = a;
  res.Scale(s);
  res += Line1a;
  Log() << Verbose(5) << "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<Vector,Vector> 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<Vector,Vector> 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<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]);
}

RealSpaceMatrix Line::getRotationMatrix(double alpha) const
{
  RealSpaceMatrix M;
  M.setRotation(getDirection(), alpha);

  return M;
}

Plane Line::getOrthogonalPlane(const Vector &origin) const{
  return Plane(getDirection(),origin);
}

std::vector<Vector> Line::getSphereIntersections() const{
  std::vector<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<double>::infinity());
}
LinePoint Line::negEndpoint() const{
  return LinePoint(*this,-numeric_limits<double>::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() << LinearAlgebraVectorPair( make_pair(&x1, &x2) );
  }
  return Line(x1,x1-x2);
}

std::ostream& operator<<(std::ostream& ost, const Line& m)
{
  const Vector origin = m.getOrigin();
  const Vector direction = m.getDirection();
  ost << "(";
  for (int i=0;i<NDIM;i++) {
    ost << origin[i];
    if (i != 2)
      ost << ",";
  }
  ost << ") -> (";
  for (int i=0;i<NDIM;i++) {
    ost << direction[i];
    if (i != 2)
      ost << ",";
  }
  ost << ")";
  return ost;
};


/******************************** Points on the line ********************/