计算几何算法二：向量/线段

向量的基本运算（加、减），点乘、单位叉乘模、点与线段的关系、线段上离某点最近的点、二线段是否相交、点与线段的面积、点到线段所在直线上的垂足。

// RtVector.h: interface for the CRtVector class.
//
//

#if !defined(AFX_RTVECTOR_H__150FF33B_0C72_4C7D_97FE_2D9ADA5AA897__INCLUDED_)
#define AFX_RTVECTOR_H__150FF33B_0C72_4C7D_97FE_2D9ADA5AA897__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "RtPoint.h"

template<class T,int dim = 2 >
class CRtVector
{
public:
 CRtVector()
 {
 }
 CRtVector( CRtPoint<T, dim>& start,CRtPoint<T, dim>& stop ): mvo_start( start ),mvo_stop( stop )
 {
 }
 virtual ~CRtVector()
 {
 }

 //copy constructor
 CRtVector(const CRtVector& in)
 {
  mvo_start = in.mvo_start;
  mvo_stop = in.mvo_stop;
 }

//operators overload:
 //Accumulative addition of two vectors
 CRtVector& operator += (const CRtVector& in )
 {
  mvo_start += in.mvo_start;
  mvo_stop += in.mvo_stop;

  return *this;
 }
 // Accumulative subtraction of two vectors
 CRtVector& operator -= (const CRtVector& in )
 {
  mvo_start -= in.mvo_start;
  mvo_stop -= in.mvo_stop;

  return *this;
 }

 //Vector Equality, epsilon used due to numerical imprecision
 bool operator == (const CRtVector& in )
 {
  if ((mvo_start == in.mvo_start) && (mvo_stop == in.mvo_stop))
  {
   return true;
  }
  
  return false;
 }
 //evaluate to point
 CRtVector& operator = (const CRtVector& in)
 {
  mvo_start = in.mvo_start;
  mvo_stop = in.mvo_stop;

  return *this;
 }
 //get a  element
 CRtPoint<T, dim> operator [] (const int index)
 {
  return index == 0 ? mvo_start:mvo_stop;
 }
 
 //Reassign a vectors without making a temporary structure
 void SetData( CRtPoint<T, dim>& start,CRtPoint<T, dim>& stop )
 {
  mvo_start = start;
  mvo_stop = stop;
 }
 
 //Dot Product
 double DotProduct( CRtVector& in );
 
 //Cross Product Mode
 double CrossProductMode( CRtPoint<T, dim>& in);

 //the relation a point and line segment
 double RelationPoint( CRtPoint<T, dim>& in);

 //Point C to the straight line segment AB where the pedal
 const CRtPoint<T, dim> PedalLine( CRtPoint<T, dim>& in);

 //Seek the shortest distance from point to line segment and return segment from that point on the nearest point of
 //Segment to the point nearest point, not necessarily Pedal
 double NearestPoint( CRtPoint<T, dim>& in,CRtPoint<T, dim>& out);

 //Whether the two line segments intersect
 bool IsIntersect( CRtVector& in );

 //the area of a point and line segments
 double AreaTriangle( CRtPoint<T>& in);

private:
 CRtPoint<T, dim> mvo_start;//vector start point
 CRtPoint<T, dim> mvo_stop;//vector stop point
};

//Adds two vectors together
template<class T, int dim>
inline const CRtVector<T, dim> operator + (const CRtVector<T, dim>& src0, const CRtVector<T, dim>& src1 )
{
 return CRtVector<T, dim>( src0.mvo_start + src1.mvo_start , src0.mvo_stop + src1.mvo_stop );
}

//Subtracts two vectors
template<class T, int dim>
inline const CRtVector<T, dim> operator - (const CRtVector<T, dim>& src0, const CRtVector<T, dim>& src1 )
{
 return CRtVector<T, dim>( src0.mvo_start - src1.mvo_start , src0.mvo_stop - src1.mvo_stop );
}

/***************************************************************************** 
Dot Product:
If two vectors are non-zero vector 
 return value < 0：The angle between two vectors is acute angle
 return value = 0：The angle between two vectors for the rectangular
 return value > 0：Angle between two vectors as the obtuse angle
*******************************************************************************/
template<class T,int dim>
double CRtVector<T, dim>::DotProduct( CRtVector& in )
{
 CRtPoint<T,dim> point;
 double val = 0;
 point = (mvo_stop - mvo_start) * (in.mvo_stop - in.mvo_start);
 T *data = point.GetData();

 for (int i = 0; i < dim; i++)
  val += *(data + i);

 return val;
}

/****************************************************************************** 
Cross Product Mode:
Unit vector of the mode is 1, the available cross product of the modulus:
Cross product of a very important property can be judged by its symbol between  the two vectors
of cis-anti-clockwise relationship:
 If the P × OQ < 0, then P in the OQ in a clockwise direction.
 If the P × OQ > 0, then P in the OQ of the counter-clockwise.
 If the P × OQ = 0, then P and OQ collinear, but may in the same direction may reverse.
*******************************************************************************/ 
template<class T,int dim>
double CRtVector<T, dim>::CrossProductMode( CRtPoint<T, dim>& in)
{
 CRtPoint<T, dim> p0,p1;
 p0 = mvo_stop - mvo_start;
 p1 = in - mvo_start;

 if (2 == dim) return p0[0] * p1[1] - p0[1] * p1[0];
 if (3 == dim) return p0[1] * p1[2] + p0[2] * p1[0] + p0[0] * p1[1]
                - p0[2] * p1[1] + p0[0] * p1[2] + p0[1] * p1[0];
}

/******************************************************************************   
   AC dot AB    (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay) 
  r = --------- =  -------------------------------
   ||AB||^2                  L^2 
 
  r has the following meaning: 
   r=0      P = A 
   r=1      P = B 
   r<0    P is on the backward extension of AB 
   r>1      P is on the forward extension of AB 
   0<r<1  P is interior to AB 
*******************************************************************************/ 
template<class T,int dim>
double CRtVector<T, dim>::RelationPoint( CRtPoint<T, dim>& in) 
{  
 CRtVector<T, dim> vector;
 vector.SetData(mvo_start, in);
 return DotProduct(vector) / pow( mvo_start.DistEuclidean(mvo_stop) , 2 ); 
}

//Point C to the straight line segment AB where the pedal
template<class T,int dim>
const CRtPoint<T, dim> CRtVector<T, dim>::PedalLine( CRtPoint<T, dim>& in) 
{ 
 double r = RelationPoint( in ); 
 CRtPoint<T, dim> point = mvo_stop - mvo_start;
 point *= r;
 point += mvo_start;
 return point; 
} 

//Seek the shortest distance from point to line segment and return segment from that point on the nearest point of
//Segment to the point nearest point, not necessarily Pedal
template<class T,int dim>
double CRtVector<T, dim>::NearestPoint( CRtPoint<T, dim>& in,CRtPoint<T, dim>& out) 
{ 
 double val = RelationPoint( in );
 
 if( val < 0 )  out = mvo_start; 
 if( val > 1 )  out = mvo_stop;  
 if( val >= 0 && val <= 1) out = PedalLine( in);
 
 return in.DistEuclidean(out);
}
/******************************************************************************
1) Fast exclusion test
 Is located in the diagonal line segment P1P2 rectangle of R,
 Is located in the diagonal line segment Q1Q2 rectangle of T,
 If the R and T do not intersect, and apparently two line segments do not intersect.

2) Cross-li test
 If the two line segments intersect, then the two line segments bound to each other cross-li.
 If the P1P2 cross-li Q1Q2, then the vector (P1 - Q1)
 And (P2 - Q1) in the vector (Q2 - Q1) on both sides,
 Namely, (P1 - Q1) × (Q2 - Q1) * (P2 - Q1) × (Q2 - Q1) <0.
 On the type can be rewritten as (P1 - Q1) × (Q2 - Q1) * (Q2 - Q1) × (P2 - Q1)> 0.
 When the (P1 - Q1) × (Q2 - Q1) = 0, the statement (P1 - Q1) and (Q2 - Q1) of lines,
 However, as it has been a quick rejection test, so a certain line segment Q1Q2 on P1;
 By the same token, (Q2 - Q1) × (P2 - Q1) = 0 shows a certain line segment Q1Q2 on the P2.
 Therefore, to determine P1P2 is based on cross-li Q1Q2:
 (P1 - Q1) × (Q2 - Q1) * (Q2 - Q1) × (P2 - Q1) <= 0.
 By the same token to determine Q1Q2 is based on cross-li P1P2:
 (Q1 - P1) × (P2 - P1) * (P2 - P1) × (Q2 - P1) <= 0.
*******************************************************************************/
template<class T,int dim>
bool CRtVector<T, dim>::IsIntersect( CRtVector& in )
{
 if (dim != 2) return false;

 if ( max(mvo_start[0], mvo_stop[0]) >= min(in.mvo_start[0], in.mvo_stop[0]) &&
   min(mvo_start[0], mvo_stop[0]) <= max(in.mvo_start[0], in.mvo_stop[0]) &&
   max(mvo_start[1], mvo_stop[1]) >= min(in.mvo_start[1], in.mvo_stop[1]) &&
   min(mvo_start[0], mvo_stop[1]) <= max(in.mvo_start[1], in.mvo_stop[1]) )//exclude
 {
  //cross-li
  double val = 0;
  CRtVector<T, dim> vector(in.mvo_stop, mvo_start);  
  val = vector.CrossProductMode(in.mvo_start);
  vector.SetData(in.mvo_stop, mvo_stop);
  val *= vector.CrossProductMode(in.mvo_start);
  
  if ( val <= 0 )
  {
   vector.SetData(in.mvo_start, mvo_stop);
   val = vector.CrossProductMode(mvo_start);
   vector.SetData(in.mvo_stop, mvo_stop);
   val *= vector.CrossProductMode(mvo_start);
   
   return val <= 0 ? true : false;
  }
 }

 return false;
}

//the area of a point and line segments
template<class T,int dim>
double CRtVector<T, dim>::AreaTriangle( CRtPoint<T>& in)
{
    return fabs( ( mvo_start[0] * mvo_stop[1] + mvo_stop[0] * in[1] + in[0] * mvo_start[1] )
         - ( mvo_start[0] * in[1] + mvo_stop[0] * mvo_start[1] + in[0] * mvo_stop[1] ) / 2.0 );  
}

#endif // !defined(AFX_RTVECTOR_H__150FF33B_0C72_4C7D_97FE_2D9ADA5AA897__INCLUDED_)