線性代數之平面

一，理論

平面方程：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.<br /> Plane(const Vec3<Type> & a_normal, const Vec3<Type> & a_pt)<br /> {<br /> m_normal = a_normal;<br /> m_normal.normalize();<br /> m_distance = m_normal.dot(a_pt);<br /> } 

2，使用平面上點序列構造平面

點序列按照逆時針順序：如果這些點共面則計算其重心，任選相鄰三點，兩個向量叉乘計演算法向轉化為問題1；如果這些點不共面，需要計算出一個 最適合的平面，使用：Newell's Method for Computing the Plane Equation of a Polygon

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

3，判斷某個點是否位於平面上，不在平面上確定位於平面哪邊、

 將該點和法向點積結果為改點到平面的距離，比較該距離和原點到平面的距離。

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

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

 

三，應用

對位於一條射線上的若干點，按照射線方向排序，相鄰頂點比較：只需要構造過一點該射線的法平面，判斷另一點位於這個平面左或者右邊即可：

struct my_greater<br />: public std::binary_function<Point, Point, bool><br />{<br />typedefPoint_Type;<br />my_greater(const Vector& d)<br />: _dir(d)<br />{}<br />// 比較兩點在射線方向的前後順序<br />// true : far_pt 確實比 near_pt離射線起點遠<br />// false : far_pt 不比 near_pt離射線起點遠<br />// functor for operator>=<br />bool operator()(const _Type& far_pt, const _Type& near_pt) const<br />{// apply operator>= to operands</p><p>// 通過原點到 "過lfs且垂直dir的平面" 的距離比較<br />// 由於Point不能直接跟Vector求點積，出於效率直接計算<br />float origin_to_near_plane_dist = _dir.x*near_pt.x + _dir.y*near_pt.y + _dir.z*near_pt.z;<br />float origin_to_far_plane_dist = far_pt.x*_dir.x + far_pt.y*_dir.y + far_pt.z*_dir.z;<br />// 距離有正負之分：<br />//+ 當原點位於平面下方時，距離為正<br />//+ 當原點位於平面上方時，距離為負<br />// 三種可能：<br />//+ 原點同時位於兩個平面的下方，此時: origin_to_far_plane_dist >= origin_to_near_plane_dist > 0.<br />//+ 原點同時位於兩個平面的上方，此時：0 > origin_to_far_plane_dist >= origin_to_near_plane_dist.<br />//+ 原點位於兩個平面之間，此時：origin_to_far_plane_dist >= 0 > origin_to_near_plane_dist.<br />return (origin_to_far_plane_dist >= origin_to_near_plane_dist);<br />}<br />const Vector&_dir;<br />private:<br />my_greater();<br />}; 

比較的過程只是一些乘積、求和和比較。

 

gtl在之前一篇blog給出連結：http://blog.csdn.net/dizuo/archive/2008/05/28/2491400.aspx