三点求圆弧中心arx

admin

1. 
   
2. 
   
3. 已知空间三个点，解算外接圆圆心坐标，C++编程实现
   
4. struct PT3
   
5. {
   
6.         double x, y, z;
   
7. };
   
8. 
   
9. int solveCenterPointOfCircle(std::vector<PT3> pd, double centerpoint[])
   
10. {
    
11.         double a1, b1, c1, d1;
    
12.         double a2, b2, c2, d2;
    
13.         double a3, b3, c3, d3;
    
14. 
    
15.         double x1 = pd[0].x, y1 = pd[0].y, z1 = pd[0].z;
    
16.         double x2 = pd[1].x, y2 = pd[1].y, z2 = pd[1].z;
    
17.         double x3 = pd[2].x, y3 = pd[2].y, z3 = pd[2].z;
    
18. 
    
19.         a1 = (y1*z2 - y2*z1 - y1*z3 + y3*z1 + y2*z3 - y3*z2);
    
20.         b1 = -(x1*z2 - x2*z1 - x1*z3 + x3*z1 + x2*z3 - x3*z2);
    
21.         c1 = (x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2);
    
22.         d1 = -(x1*y2*z3 - x1*y3*z2 - x2*y1*z3 + x2*y3*z1 + x3*y1*z2 - x3*y2*z1);
    
23. 
    
24.         a2 = 2 * (x2 - x1);
    
25.         b2 = 2 * (y2 - y1);
    
26.         c2 = 2 * (z2 - z1);
    
27.         d2 = x1 * x1 + y1 * y1 + z1 * z1 - x2 * x2 - y2 * y2 - z2 * z2;
    
28. 
    
29.         a3 = 2 * (x3 - x1);
    
30.         b3 = 2 * (y3 - y1);
    
31.         c3 = 2 * (z3 - z1);
    
32.         d3 = x1 * x1 + y1 * y1 + z1 * z1 - x3 * x3 - y3 * y3 - z3 * z3;
    
33. 
    
34.         centerpoint[0] = -(b1*c2*d3 - b1*c3*d2 - b2*c1*d3 + b2*c3*d1 + b3*c1*d2 - b3*c2*d1)
    
35.                 /(a1*b2*c3 - a1*b3*c2 - a2*b1*c3 + a2*b3*c1 + a3*b1*c2 - a3*b2*c1);
    
36.         centerpoint[1] =  (a1*c2*d3 - a1*c3*d2 - a2*c1*d3 + a2*c3*d1 + a3*c1*d2 - a3*c2*d1)
    
37.                 /(a1*b2*c3 - a1*b3*c2 - a2*b1*c3 + a2*b3*c1 + a3*b1*c2 - a3*b2*c1);
    
38.         centerpoint[2] = -(a1*b2*d3 - a1*b3*d2 - a2*b1*d3 + a2*b3*d1 + a3*b1*d2 - a3*b2*d1)
    
39.         centerpoint[2] = -(a1*b2*d3 - a1*b3*d2 - a2*b1*d3 + a2*b3*d1 + a3*b1*d2 - a3*b2*d1)
    
40.                 /(a1*b2*c3 - a1*b3*c2 - a2*b1*c3 + a2*b3*c1 + a3*b1*c2 - a3*b2*c1);
    
41. 
    
42.         return 0;
    
43. }
    
44. 
    
45. 
    
46. 
    

复制代码

admin

1. 根据三点创建圆弧
   
2. 
   
3. 
   
4. private Point3d GetArcCenter(Point3d pt1,Point3d pt2,Point3d pt3,ref double radius)
   
5.         {
   
6.             double xysm, xyse, xy;
   
7.             double[] m_ArcCenter = new double[3];
   
8.             xy = Math.Pow(pt1.X, 2) + Math.Pow(pt1.Y, 2);
   
9.             xyse = xy - Math.Pow(pt3.X, 2) - Math.Pow(pt3.Y, 2);
   
10.             xysm = xy - Math.Pow(pt2.X, 2) - Math.Pow(pt2.Y, 2);
    
11.             xy = (pt1.X - pt2.X) * (pt1.Y - pt3.Y) - (pt1.X - pt3.X) * (pt1.Y - pt2.Y);
    
12.             // 判断参数有效性
    
13.             if (Math.Abs(xy)<0.000001)
    
14.             {
    
15.                 MessageBox.Show("所输入的参数无法创建圆形!");
    
16.                 return Point3d.Origin;
    
17.             }
    
18. 
    
19.             // 获得圆心和半径
    
20.             m_ArcCenter[0] = ((xysm * (pt1.Y - pt3.Y)) - (xyse * (pt1.Y - pt2.Y))) / (2 * xy);
    
21.             m_ArcCenter[1] = ((xyse * (pt1.X - pt2.X)) - (xysm * (pt1.X - pt3.X))) / (2 * xy);
    
22.             m_ArcCenter[2] = 0;
    
23.             radius = Math.Sqrt((pt1.X - m_ArcCenter[0]) * (pt1.X - m_ArcCenter[0]) + (pt1.Y - m_ArcCenter[1]) * (pt1.Y - m_ArcCenter[1]));
    
24.             if (radius<0.000001)
    
25.             {
    
26.                 MessageBox.Show("半径过小");
    
27.                 return Point3d.Origin;
    
28.             }
    
29. 
    
30.             // 函数返回圆心的位置，而半径在参数中通过引用方式返回
    
31.             Point3d arcCenter = new Point3d(m_ArcCenter[0], m_ArcCenter[1], m_ArcCenter[2]);
    
32.             return arcCenter;
    
33.         }
    
34. 
    
35.         public Arc AddArcFromPt(Point3d ptSt,Point3d ptSc,Point3d ptEn,ObjectId objid,string featurename)
    
36.         {
    
37.             // 三点法创建圆弧
    
38.             Point3d ptCenter;
    
39.             double radius = 0.0;
    
40.             Arc objArc;
    
41.             ptCenter = GetArcCenter(ptSt, ptSc, ptEn, ref radius);
    
42.             if (isClockWise(ptCenter, ptSt, ptSc, ptEn))
    
43.                 objArc = AddArcCSEP(ptCenter, ptSt, ptEn, radius,objid,featurename);
    
44.             else
    
45.                 objArc = AddArcCSEP(ptCenter, ptEn, ptSt, radius,objid,featurename);
    
46.             return objArc;
    
47.         }
    
48. 
    
49.         private Arc AddArcCSEP(Point3d ptCenter,Point3d ptSt,Point3d ptEn,double radius,ObjectId objid,string featurename)
    
50.         {
    
51.             double radiusent;
    
52.             double stAng, enAng;
    
53.             // 计算半径
    
54.             radiusent = radius;
    
55.             // 计算起点角度和终点角度
    
56.             stAng = CalLineAngle(ptCenter, ptSt);
    
57.             enAng = CalLineAngle(ptCenter, ptEn);
    
58.             Arc objArc = new Arc(ptCenter, radius, stAng, enAng);
    
59.             objArc.LinetypeId = objid;
    
60.             objArc.Layer = featurename;
    
61.             objArc.SetDatabaseDefaults();
    
62.             return objArc;
    
63.         }
    
64. 
    
65.         private bool isClockWise(Point3d ptCenter,Point3d ptSt,Point3d ptSc,Point3d ptEn)
    
66.         {
    
67.             double a1, a2, a3;
    
68.             a1 = CalLineAngle(ptSt, ptCenter);
    
69.             a2 = CalLineAngle(ptSc, ptCenter);
    
70.             a3 = CalLineAngle(ptEn, ptCenter);
    
71.             bool isClock = (a1 < a2) | (a2 < a3) | (a1 < a3);
    
72.             return isClock;
    
73.         }
    
74. 
    
75.         private double CalLineAngle(Point3d vStartPt,Point3d vEndPt)
    
76.         {
    
77.             double X, Y;
    
78.             double CalculateLineAngle = 0.0;
    
79.             double dbJudgeAngle = 0.001;
    
80.             X = vEndPt.X - vStartPt.X;
    
81.             Y = vEndPt.Y - vStartPt.Y;
    
82.             if ((X > 0 && Y > 0)||(X > 0 && Y < 0))
    
83.             {
    
84.                 CalculateLineAngle = Math.Atan(Y / X);
    
85.             }
    
86.             else if ((X < 0 && Y > 0) || (X < 0 && Y < 0))
    
87.             {
    
88.                 CalculateLineAngle = Math.PI + Math.Atan(Y / X);
    
89.             }
    
90.             else if(X == 0)
    
91.             {
    
92.                 if (Y < 0)
    
93.                 {
    
94.                     CalculateLineAngle = -Math.PI / 2;
    
95.                 }
    
96.                 else if (Y > 0)
    
97.                 {
    
98.                     CalculateLineAngle = Math.PI / 2;
    
99.                 }
    
100.             }
     
101.             else if (Y == 0)
     
102.             {
     
103.                 if (X < 0)
     
104.                 {
     
105.                     CalculateLineAngle = Math.PI;
     
106.                 }
     
107.                 else if (X > 0)
     
108.                 {
     
109.                     CalculateLineAngle = 0;
     
110.                 }
     
111.             }
     
112.             else
     
113.             {
     
114.                 ;
     
115.             }
     
116. 
     
117.             if (Math.Abs(CalculateLineAngle - 0) < dbJudgeAngle)
     
118.             {
     
119.                 CalculateLineAngle = 0;
     
120.                 return CalculateLineAngle;
     
121.             }
     
122.             else if (CalculateLineAngle < 0)
     
123.             {
     
124.                 CalculateLineAngle = 2 * Math.PI + CalculateLineAngle;
     
125.                 return CalculateLineAngle;
     
126.             }
     
127.             else
     
128.                 return CalculateLineAngle;
     
129.         }
     
130. 
     
131. #autoca

复制代码

admin

已知圆弧的一些信息,求得圆弧上中间点的坐标的方法(C++语言描述)

1. // 从圆弧一些信息得到圆弧上一个特殊的点(计算的是二维图形的情况)
   
2. Point3d DBOPERATION::GetArcTangencyPoint(Point3d pStartPoint, Point3d pEndPoint,
   
3.                                          Point3d pCenterPoint,  double  dRadius,
   
4.                                           double  dStartAngle,  double  dSweepAngle)
   
5. {
   
6.     Point3d pResultPoint;
   
7.     pResultPoint.X  =   0 ;
   
8.     pResultPoint.Y  =   0 ;
   
9.     pResultPoint.Z  =   0 ;
   
10. 
    
11.      // 因为会得到两个点,哪个点在弧上需要再进行判断
    
12.      double  dRx1  =   0 ;
    
13.      double  dRx2  =   0 ;
    
14.      double  dRy1  =   0 ;
    
15.      double  dRy2  =   0 ;
    
16. 
    
17.      const   double  PI  =   3.1415926535897932 ;
    
18. 
    
19.     Point3d pMiddlePoint;
    
20.     pMiddlePoint.X  =  (pStartPoint.X  +  pEndPoint.X)  /   2 ;
    
21.     pMiddlePoint.Y  =  (pStartPoint.Y  +  pEndPoint.Y)  /   2 ;
    
22.     pMiddlePoint.Z  =   0.0 ;
    
23. 
    
24.      // 扫角的角度值
    
25.      double  dArcSweepAngle  =   180   *  dSweepAngle  /  PI;
    
26.      // 中点到圆心的距离
    
27.      double  ddy  =   0 ;
    
28.      double  ddx  =   0 ;
    
29.      if  (pMiddlePoint.X  -  pCenterPoint.X  >=   0.1e-6 )
    
30.     {
    
31.         ddx  =  pMiddlePoint.X  -  pCenterPoint.X;
    
32.     }
    
33.      if  (pMiddlePoint.Y  -  pCenterPoint.Y  >=   0.1e-6 )
    
34.     {
    
35.         ddx  =  pMiddlePoint.Y  -  pCenterPoint.Y;
    
36.     }
    
37.      double  dDistance  =  sqrt(ddx  *  ddx  +  ddy  *  ddy);
    
38. 
    
39.      if  ((fabs(pMiddlePoint.Y  -  pCenterPoint.Y)  <   0.1e-6
    
40.          &&  fabs(pMiddlePoint.X  -  pCenterPoint.X)  <   0.1e-6 )
    
41.          &&  (fabs(pStartPoint.X  -  pEndPoint.X)  <   0.1e-6
    
42.          ||  fabs(pStartPoint.Y  -  pEndPoint.Y)  < 0.1e-6  ))
    
43.     {
    
44.          // 半圆且斜率=0或者斜率为无穷大的情况
    
45.          if  (fabs(pStartPoint.X  -  pEndPoint.X)  <   0.1e-6 )
    
46.         {
    
47.              // 结果的两个X的值
    
48.             dRx1  =  pMiddlePoint.X  +  (dRadius  -  dDistance);
    
49.             dRx2  =  pMiddlePoint.X  -  (dRadius  +  dDistance);
    
50.              // 对应的两个Y的值
    
51.             dRy1  =  pMiddlePoint.Y;
    
52.             dRy2  =  pMiddlePoint.Y;
    
53.         }
    
54.          if  (fabs(pStartPoint.Y  -  pEndPoint.Y)  < 0.1e-6 )
    
55.         {
    
56.              // 结果的两个X的值
    
57.             dRx1  =  pMiddlePoint.X;
    
58.             dRx2  =  pMiddlePoint.X;
    
59.              // 对应的两个Y的值
    
60.             dRy1  =  pMiddlePoint.Y  +  (dRadius  -  dDistance);
    
61.             dRy2  =  pMiddlePoint.Y  -  (dRadius  +  dDistance);
    
62.         }
    
63.     }
    
64.      else   if  (fabs(pMiddlePoint.Y  -  pCenterPoint.Y)  <   0.1e-6
    
65.          &&  fabs(pMiddlePoint.X  -  pCenterPoint.X)  >=   0.1e-6 )
    
66.     {
    
67.          // 圆心与端点中心点在水平线上的情况,斜率=0
    
68.          // 结果的两个X的值
    
69.         dRx1  =  pMiddlePoint.X  +  (dRadius  -  dDistance);
    
70.         dRx2  =  pMiddlePoint.X  -  (dRadius  +  dDistance);
    
71.          // 对应的两个Y的值
    
72.         dRy1  =  pMiddlePoint.Y;
    
73.         dRy2  =  pMiddlePoint.Y;
    
74.     }
    
75.      else   if  (fabs(pMiddlePoint.X  -  pCenterPoint.X)  <   0.1e-6
    
76.          &&  fabs(pMiddlePoint.Y  -  pCenterPoint.Y)  >=   0.1e-6 )
    
77.     {
    
78.          // 圆心与端点中心点在垂直线上的情况,斜率=无穷大
    
79.          // 结果的两个X的值
    
80.         dRx1  =  pMiddlePoint.X;
    
81.         dRx2  =  pMiddlePoint.X;
    
82.          // 对应的两个Y的值
    
83.         dRy1  =  pMiddlePoint.Y  +  (dRadius  -  dDistance);
    
84.         dRy2  =  pMiddlePoint.Y  -  (dRadius  +  dDistance);
    
85.     }
    
86.      else
    
87.     {
    
88.          // 圆弧两个端点成线的斜率,已经排除了等于0或者无穷大的可能性
    
89.          double  dTheKSAE  =  (pEndPoint.Y  -  pStartPoint.Y)  /  (pEndPoint.X  -  pStartPoint.X);
    
90.          // 与端点组成的线垂直的线的斜率(0度,180度,360度弧的dk需要特殊处理)
    
91.          double  dK  =  ( 1   /  dTheKSAE)  *  ( - 1 );
    
92. 
    
93.          if  (fabs(fabs(dArcSweepAngle)  -   360 )  <   0.1e-6 )
    
94.         {
    
95.             dK  =  (pMiddlePoint.Y  -  pCenterPoint.Y)  /  (pMiddlePoint.X  -  pCenterPoint.X);
    
96.         }
    
97.          if  (fabs(fabs(dArcSweepAngle)  -   180 )  <   0.1e-6 )
    
98.         {
    
99.             dK  =  (pStartPoint.Y  -  pCenterPoint.Y)  /  (pStartPoint.X  -  pCenterPoint.X);
    
100.             dK  =  ( 1   /  dK)  *  ( - 1 );
     
101.         }
     
102.          if  (fabs(fabs(dArcSweepAngle)  -   0 )  <   0.1e-6 )
     
103.         {
     
104.             dK  =  (pMiddlePoint.Y  -  pCenterPoint.Y)  /  (pMiddlePoint.X  -  pCenterPoint.X);
     
105.         }
     
106.          /
     
107.          // 这是经过两端点中点并与两端点组成的直线垂直的直线的方程
     
108.          // k为斜率,M(x),M(y)分别为中心点坐标值
     
109.          // y=k * x - k * M(x) + M(y)
     
110.          /
     
111. 
     
112.          // 求出- k * M(x) + M(y)部分的值
     
113.          double  dDif  =  ( - 1 )  *  dK  *  pMiddlePoint.X  +  pMiddlePoint.Y;
     
114. 
     
115.          /
     
116.          //  [x - C(x)]^2 + [y - C(y)]^2 = r^2
     
117.          //  代入上面的y,其中- k * M(x) + M(y)以dDif代替
     
118.          //  [x - C(x)]^2 + [kx + dDif - C(y)]^2 = r^2
     
119.          //  下面以dBetween代替+ dDif - C(y)
     
120.          /
     
121.          double  dBetween  =  dDif  -  pCenterPoint.Y;
     
122. 
     
123.          /
     
124.          //  x^2 + 2*C(x)*x + [C(x)]^2
     
125.          //  +
     
126.          //  (kx)^2 + 2*k*dBetween*x + (dBetween)^2
     
127.          //  = r^2
     
128.          /
     
129. 
     
130.          // 再求中间值
     
131.          double  dNx  =  ( 2   *  dK  *  dBetween  -   2   *  pCenterPoint.X)  /  ( 1   +  dK  *  dK);
     
132.          double  dLa  =  (dRadius  *  dRadius  -  pCenterPoint.X  *  pCenterPoint.X
     
133.              -  dBetween  *  dBetween)  /  ( 1   +  dK  *  dK);
     
134.          double  dAd  =  dLa  +  dNx  *  dNx  /   4 ;
     
135. 
     
136.          // 结果的两个X的值
     
137.         dRx1  =  sqrt(dAd)  -  dNx  /   2 ;
     
138.         dRx2  =  ( - 1 )  *  sqrt(dAd)  -  dNx  /   2 ;
     
139.          // 对应的两个Y的值
     
140.         dRy1  =  dK  *  dRx1  -  dK  *  pMiddlePoint.X  +  pMiddlePoint.Y;
     
141.         dRy2  =  dK  *  dRx2  -  dK  *  pMiddlePoint.X  +  pMiddlePoint.Y;
     
142.     }
     
143. 
     
144.      // 对得到的两个点进行判断,找到在弧上的那个
     
145.     pResultPoint.X  =  dRx1;
     
146.     pResultPoint.Y  =  dRy1;
     
147.      if  (fabs(fabs(dArcSweepAngle)  -   180.0 )  <   0.1e-4 )
     
148.     {
     
149.          // 半圆情况处理
     
150.          double  p1x,p1y,p2x,p2y,qx,qy,dx1,dy1,dx2,dy2,det;
     
151.         p2x  =  pStartPoint.X;
     
152.         p2y  =  pStartPoint.Y;
     
153.         p1x  =  pCenterPoint.X;
     
154.         p1y  =  pCenterPoint.Y;
     
155.         qx  =  dRx2;
     
156.         qy  =  dRy2;
     
157. 
     
158.         dx1  =  p2x  -  p1x;
     
159.         dy1  =  p2y  -  p1y;
     
160.         dx2  =  qx -  p2x;
     
161.         dy2  =  qy  -  p2y;
     
162.         det  =  dx1 * dy2  -  dx2 * dy1;
     
163.         BOOL bBlockWise  =  TRUE;
     
164.          if  (det  >   0.0 )
     
165.         {
     
166.              // 逆时针方向
     
167.             bBlockWise  =  FALSE;
     
168.         }
     
169.          if  ((( ! bBlockWise)  &&  (dArcSweepAngle  >   0.0 ))
     
170.              ||  ((bBlockWise)  &&  (dArcSweepAngle  <   0.0 )))
     
171.         {
     
172.             pResultPoint.X  =  dRx2;
     
173.             pResultPoint.Y  =  dRy2;
     
174.         }
     
175.     }
     
176.      else
     
177.     {
     
178.          // 非半圆的情况
     
179.          double  dd1  =  sqrt(pow(dRx1  -  pMiddlePoint.X,  2 )
     
180.              +  pow(dRy1  -  pMiddlePoint.Y,  2 ));
     
181.          double  dd2  =  sqrt(pow(dRx2  -  pMiddlePoint.X,  2 )
     
182.              +  pow(dRy2  -  pMiddlePoint.Y,  2 ));
     
183.          if  (((dd1  <  dd2)  &&  (fabs(dArcSweepAngle)  >   180.0 ))
     
184.              ||  ((dd1  >  dd2)  &&  (fabs(dArcSweepAngle)  <   180.0 )))
     
185.         {
     
186.             pResultPoint.X  =  dRx2;
     
187.             pResultPoint.Y  =  dRy2;
     
188.         }
     
189.     }
     
190. 
     
191.      return  pResultPoint;
     
192. }

复制代码