线性代数之平面

一，理论

平面方程：Ax + By + Cz + w = 0;

这个方程的意义，可以看成两个向量的点击 N(A,B,C)  V(x,y,z), N.dot(V) = -w; 向量N是一个常量向量，向量V代表一类向量。N可以看成平面的单位法向量，V可以看成平面上的点和原点的连线向量。则此时的w就有了特殊含义：原点到平面的距离。

如下图， N.dot(V) = A*x + B*y + C*z = |N| * |V| * cos(angle) = |V| * cos(angle) = distance; 

角度angle为平面法线和向量V的夹角。

 

所以平面方程可以表示为法向量和距离的二元组。

plane = <Normal, distance>

二，代码实现

根据上面的理论，设计一个平面类只需要两个参数：

Plane

      + distance : float

      + normal : Vector3

 

下面介绍一些常用函数：

1，使用平面的法线和平面上一点构造平面

      根据上图所示，该点和单位化以后的平面法向点积就是原点到平面的距离。

//! Construct a Plane from a normal and a point laying in the plane. /a normal must not be null. Plane(const Vec3<Type> & a_normal, const Vec3<Type> & a_pt) { m_normal = a_normal; m_normal.normalize(); m_distance = m_normal.dot(a_pt); } 

2，使用平面上点序列构造平面

点序列按照逆时针顺序：如果这些点共面则计算其重心，任选相邻三点，两个向量叉乘计算法向转化为问题1；如果这些点不共面，需要计算出一个 最适合的平面，使用：Newell's Method for Computing the Plane Equation of a Polygon

   Plane(const std::vector< Vec3<Type> > & points) { if(points.empty()) return; // "Newell's Method for Computing the Plane Equation of a Polygon", // Graphics Gems III, Academic Press (1992), pp.231-232. m_normal.setValue(0,0,0); Vec3<Type> refpt(0,0,0); // Compute the polygon normal and a reference point on // the plane. Note that the actual reference point is // refpt / nverts for(unsigned int i=0; i<points.size(); i++){ const Vec3<Type> & u = points[i]; const Vec3<Type> & v = points[(i+1)%points.size()]; m_normal[0] += (u[1] - v[1]) * (u[2] + v[2]); m_normal[1] += (u[2] - v[2]) * (u[0] + v[0]); m_normal[2] += (u[0] - v[0]) * (u[1] + v[1]); refpt += u; } // Compute the distance from origin to the plane equation m_distance = -refpt.dot(m_normal) / (m_normal.length() * points.size()); } 

3，判断某个点是否位于平面上，不在平面上确定位于平面哪边、

 将该点和法向点积结果为改点到平面的距离，比较该距离和原点到平面的距离。

//! Check if the given point lies in the halfspace of the plane which the plane normal vector is pointing. bool isInHalfSpace(const Vec3<Type> & a_point) const { return (a_point.dot(m_normal) >= m_distance); } 

4，判断射线Ray和平面P是否相交。

先根据射线方向判断是否和平面相交

angle_cos = Ray::dir.dot( Plane::normal ); 两个单位向量点积，计算射线和平面法线的夹角余弦值。

如果射线跟平面相交，则

d1 = Ray::origin.dot( Plane::normal ) 这句话的含义：计算原点到 经过origin且垂直于normal的平面P1的距离。

d = Plane::distance - d1; 这句话含义：计算平行平面P和P1的距离。

t = d / angle_cos 就是射线方向上起点到平面的距离t, t有正负。

射线上任何一点都可以表示为 pt = origin + t * dir; 射线和平面交点很容易求出。

/*! Ray-plane intersection. Return true if there is an intersection. /a t is the distance from ray origin to plane. */ bool intersect(const Ray<Type> & a_ray, Type & t) const { // Check if the ray is parallel to the plane. Type denom = m_normal.dot(a_ray.getDirection()); if(denom == (Type)0.0) return false; t = (m_distance - m_normal.dot(a_ray.getOrigin())) / denom; return true; } 

 

三，应用

对位于一条射线上的若干点，按照射线方向排序，相邻顶点比较：只需要构造过一点该射线的法平面，判断另一点位于这个平面左或者右边即可：

struct my_greater : public std::binary_function<Point, Point, bool> { typedef Point _Type; my_greater(const Vector& d) : _dir(d) {} // 比较两点在射线方向的前后顺序 // true : far_pt 确实比 near_pt离射线起点远 // false : far_pt 不比 near_pt离射线起点远 // functor for operator>= bool operator()(const _Type& far_pt, const _Type& near_pt) const { // apply operator>= to operands // 通过原点到 "过lfs且垂直dir的平面" 的距离比较 // 由于Point不能直接跟Vector求点积，出于效率直接计算 float origin_to_near_plane_dist = _dir.x*near_pt.x + _dir.y*near_pt.y + _dir.z*near_pt.z; float origin_to_far_plane_dist = far_pt.x*_dir.x + far_pt.y*_dir.y + far_pt.z*_dir.z; // 距离有正负之分： // + 当原点位于平面下方时，距离为正 // + 当原点位于平面上方时，距离为负 // 三种可能： // + 原点同时位于两个平面的下方，此时: origin_to_far_plane_dist >= origin_to_near_plane_dist > 0. // + 原点同时位于两个平面的上方，此时：0 > origin_to_far_plane_dist >= origin_to_near_plane_dist. // + 原点位于两个平面之间，此时：origin_to_far_plane_dist >= 0 > origin_to_near_plane_dist. return (origin_to_far_plane_dist >= origin_to_near_plane_dist); } const Vector& _dir; private: my_greater(); }; 

比较的过程只是一些乘积、求和和比较。

 

gtl在之前一篇blog给出链接：http://blog.youkuaiyun.com/dizuo/archive/2008/05/28/2491400.aspx