空间3点，计算圆心坐标

1. bool CC3PDlg::calCircleCenter(PD3Point p1,PD3Point p2,double &radium)   
2. //input p1[0] p1[1] p1[2]    
3. //output p2[0]  circle center coordinate    
4. //       p2[1]  normal     
5. //       radium      
6. {   
7.     double dPx,dQx,dRx;   
8.     double dPy,dQy,dRy;   
9.     double dPz,dQz,dRz;   
10.     double dX0,dY0,dZ0;   
11.     double dR;   
12.    
13.     dPx=p1[0].x;   
14.     dQx=p1[1].x;   
15.     dRx=p1[2].x;   
16.    
17.     dPy=p1[0].y;   
18.     dQy=p1[1].y;   
19.     dRy=p1[2].y;   
20.    
21.     dPz=p1[0].z;   
22.     dQz=p1[1].z;   
23.     dRz=p1[2].z;   
24.    
25.     //算平面法量    
26.     double pi,pj,pk;   
27.     double x1=dQx-dPx;   
28.     double x2=dRx-dPx;   
29.    
30.     double y1=dQy-dPy;   
31.     double y2=dRy-dPy;   
32.    
33.     double z1=dQz-dPz;   
34.     double z2=dRz-dPz;   
35.        
36.     pi=y1*z2-z1*y2;   
37.     pj=z1*x2-x1*z2;   
38.     pk=x1*y2-y1*x2;   
39.    
40.     if((pi==0)&&(pj==0)&&(pk==0))   
41.     {   
42.         return FALSE;   
43.     }   
44.     //求PQ和PR的中垂线    
45.     //1，过PQ的中点(Mx,My,Mz),    
46.     double dMx,dMy,dMz;   
47.        
48.     dMx=(dPx+dQx)/2;   
49.     dMy=(dPy+dQy)/2;   
50.     dMz=(dPz+dQz)/2;   
51.    
52.     //与（Mi,Mj,Mk）＝（pi，pj，pk）×（x1,y1,z1）垂直    
53.     double dMi,dMj,dMk;   
54.     dMi=pj*z1-pk*y1;   
55.     dMj=pk*x1-pi*z1;   
56.     dMk=pi*y1-pj*x1;   
57.    
58.     //2，过PR的中点(Nx,Ny,Nz),    
59.     double dNx,dNy,dNz;   
60.        
61.     dNx=(dPx+dRx)/2;   
62.     dNy=(dPy+dRy)/2;   
63.     dNz=(dPz+dRz)/2;   
64.    
65.     //与（Ni,Nj,Nk）＝（pi，pj，pk）×（x2,y2,z2）垂直    
66.     double dNi,dNj,dNk,ds;   
67.        
68.     dNi=pj*z2-pk*y2;   
69.     dNj=pk*x2-pi*z2;   
70.     dNk=pi*y2-pj*x2;   
71.    
72.     //归一    
73. //  ds=dNi*dNi+dNj*dNj+dNk*dNk;    
74.     ds=pi*pi+pj*pj+pk*pk;   
75.     ds=sqrt(ds);   
76.    
77.     if(ds==0)    
78.     {   
79.         p2[1].x=0;   
80.         p2[1].y=0;   
81.         p2[1].z=0;   
82.     }   
83.     else   
84.     {   
85.         p2[1].x=pi/ds;   
86.         p2[1].y=pj/ds;   
87.         p2[1].z=pk/ds;   
88.     }   
89.        
90.    
91.    
92.     //得到的两条中垂线为{X=dMx+dMi*tm;Y=dMy+dMj*tm;Z=dMz+dMk*tm;}  或(X-dMx)/dMi=(Y-dMy)/dMj=(Z-dMz)/dMk;    
93.     //                  {X=dNx+dNi*tn;Y=dNy+dNj*tn;Z=dNz+dNk*tn;}  或(X-dNx)/dNi=(Y-dNy)/dNj=(Z-dNz)/dNk;    
94.    
95.     //解两直线交点    
96.     double tm,tn;   
97.     tn=((dMy-dNy)*dMi+dMj*(dNx-dMx))/(dNj*dMi-dMj*dNi);   
98.     tm=(dNx+dNi*tn-dMx)/dMi;   
99.    
100.     dX0=dMx+dMi*tm;   
101.     dY0=dMy+dMj*tm;   
102.     dZ0=dMz+dMk*tm;   
103.    
104.     p2[0].x=dX0;   
105.     p2[0].y=dY0;   
106.     p2[0].z=dZ0;   
107.    
108.    
109.     //得半径    
110.     dR=(dX0-dPx)*(dX0-dPx)+(dY0-dPy)*(dY0-dPy)+(dZ0-dPz)*(dZ0-dPz);   
111.     dR=sqrt(dR);   
112.    
113.     radium=dR;   
114.     return TRUE;   
115.    
116. }