代码改变世界

三点定面的算法实现

2011-05-24 10:05  爱车龟速兔  阅读(1033)  评论(0编辑  收藏  举报
/Files/dunnice/GeometryHelper.txt
在做Gis的时候, 遇到一个通过三点定面的问题,在求解面公式的时候, 没有找到简单的第三方控件.就自己写了一个简单的方法.

三点定面公式:
Ax + By + Cz + D =0
x: x轴值
y: y轴值
z: z轴值

以下是数值版的逻辑推理版本:
已知三点P(1,-1,1) Q(2,4,6) S(1,0,3)

求出两个向量PQ(1,5,5)和PS(0,1,2)两条线

三点确定一个平面
设平面上任意一点T的坐标为(x,y,z)
则PT可以表示成aPQ+bPS  其中a,b为任意常数
即(x-1,y+1,z-1)=a(1,5,5)+b(0,1,2)=(a,5a+b,5a+2b)
所以
x-1=a
y+1=5a+b
z-1=5a+2b
把第一个式子代入第二个式子得
y+1=5(x-1)+b
b=y+1-5x+5=y-5x+6
代入第三个式子
z-1=5(x-1)+2(y-5x+6)
整理得
5x-2y+z-8=0

该平面的公式为
5x-2y+z-8=0




以下是纯逻辑推理版本:
三个点:
(x1, y1, z1) (x2, y2, z2) (x2, y3, z3)

向量坐标的基准点为(x1, y1, z1):
PQ(x2-x1, y2-y1, z2-z1)  =>  PQ(x4, y4, z4)
PS(x3-x1, y3-y1, z3-z1)  =>  PS(x5, y5, z5)

对于平面中人和一个点(x, y, z):
(x-x1, y-y1, z-z1) = a*PQ + b*PS => a*(x4, y4, z4) + b*(x5, y5, z5)
                   => (a*x4+b*x5, a*y4+b*y5, a*z4+b*z5)

所以:
x-x1 = a*x4 + b*x5
y-y1 = a*y4 + b*y5
z-z1 = a*z4 + b*z5

有第一个公式得出:
a = (x - x1 - b*x5)/x4
  => x/x4 - x1/x4 - b*x5/x4

将a的转换公式代入第二个公式:
y-y1 = (x/x4 - x1/x4 - b*x5/x4)*y4 + b*y5
     => x*y4/x4 - x1*y4/x4 - b*x5*y4/x4 + b*y5
     => x*y4/x4 - x1*y4/x4 + b*(y5 - x5*y4/x4)
得出:
b = (y - y1 - x*y4/x4 + x1*y4/x4)/(y5 - x5*y4/x4)

将a的转换公式代入第三个公式:
z-z1 = (x/x4 - x1/x4 - b*x5/x4)*z4 + b*z5
     => x*z4/x4 - x1*z4/x4 - b*x5*z4/x4 + b*z5
     => x*z4/x4 - x1*z4/x4 + b*(z5 - x5*z4/x4)

将b的转换公式代入第三个公式得到:
z-z1 = x*z4/x4 - x1*z4/x4  + ((y - y1 - x*y4/x4 + x1*y4/x4)/(y5 - x5*y4/x4))*(z5 - x5*z4/x4)
     => x*z4/x4 - x1*z4/x4  + ((y - y1 - x*y4/x4 + x1*y4/x4)/(y5 - x5*y4/x4))*(z5 - x5*z4/x4)
     => x*z4/x4 - x1*z4/x4  + ((y - y1 - x*y4/x4 + x1*y4/x4)/(y5 - x5*y4/x4))*(z5 - x5*z4/x4)

设(y5 - x5*y4/x4) 为 k
设(z5 - x5*z4/x4) 为 j
设(x1*y4/x4- y1)  为 l

z-z1 => x*z4/x4 - x1*z4/x4  + ((y  - x*y4/x4 + l)/k)*j  
     => x*z4/x4 - x1*z4/x4  + y*j/k - x*y4*j/(x4*k) + l*j/k

所以平面公式为:
-x*z4/x4 + x1*z4/x4  - y*j/k + x*y4*j/(x4*k) - l*j/k + z - z1 = 0

将常量放到公式末尾:
x*(y4*j/(x4*k) - *z4/x4) - y*j/k + z + (x1*z4/x4 - l*j/k - z1) = 0

所以
A = j*y4/(x4*k) - z4/x4
B = -j/k
C = 1
D = x1*z4/x4 - l*j/k - z1

对应的程序就非常简单了. 自己写就行了.
 
 
后来才发现, 上面的方法, 在求任何两个点的x坐标和y坐标都不相同的时候, 是可以的, 但是如果有相同的情况, 那就处理不了了.
所以还是需要采用线性代数的方式进行求解.
 
Matrix类型如下:
 
具体求解的代码如下:
 
这个确实好使.