1 // Circle.cpp
2 // Copyright 2011, Dan Heeks
3 // This program is released under the BSD license. See the file COPYING for details.
4
5 // Circle.cpp
6 // Copyright 2011, Dan Heeks
7 // This program is released under the BSD license. See the file COPYING for details.
8
9 #include "Circle.h"
10
CircleOrLine(const Point & p0,const Point & p1,const Point & p2)11 CircleOrLine::CircleOrLine(const Point& p0, const Point& p1, const Point& p2)
12 {
13 // from TangentCircles in http://code.google.com/p/heekscad/source/browse/trunk/src/Geom.cpp
14
15 // set default values, in case this fails
16 m_radius = 0.0;
17 m_c = Point(0, 0);
18
19 if (p0 == p2)
20 {
21 if (p1 == p0)
22 {
23 m_is_a_line = true;
24 m_p0 = p0;
25 m_p1 = p0;
26 }
27 else
28 {
29 m_is_a_line = false;
30 m_c = (p0 + p1) * 0.5;
31 m_radius = p0.dist(m_c);
32 }
33 return;
34 }
35
36 geoff_geometry::CLine line(geoff_geometry::Point(p0.x, p0.y), geoff_geometry::Point(p2.x, p2.y));
37 if(fabs(line.Dist(geoff_geometry::Point(p1.x, p1.y))) <= Point::tolerance)
38 {
39 m_is_a_line = true;
40 m_p0 = p0;
41 m_p1 = p2;
42 return;
43 }
44
45 m_is_a_line = false;
46
47 double x1 = p0.x;
48 double y1 = p0.y;
49 double x2 = p1.x;
50 double y2 = p1.y;
51 double x3 = p2.x;
52 double y3 = p2.y;
53
54 double a = 2 * (x1 - x2);
55 double b = 2 * (y1 - y2);
56 double d = (x1 * x1 + y1 * y1) - (x2 * x2 + y2 * y2);
57
58 double A = 2 * (x1 - x3);
59 double B = 2 * (y1 - y3);
60 double D = (x1 * x1 + y1 * y1) - (x3 * x3 + y3 * y3);
61
62 double aBmbA = (a*B - b * A); // aB - bA
63
64 // x = k + Kr where
65 double k = (B*d - b * D) / aBmbA;
66
67 // y = l + Lr where
68 double l = (-A * d + a * D) / aBmbA;
69
70 double qa = -1;
71 double qb = 0.0;
72 double qc = k * k + x1 * x1 - 2 * k*x1 + l * l + y1 * y1 - 2 * l*y1;
73
74 // solve the quadratic equation, r = (-b +- sqrt(b*b - 4*a*c))/(2 * a)
75 for (int qs = 0; qs < 2; qs++) {
76 double bb = qb * qb;
77 double ac4 = 4 * qa*qc;
78 if (ac4 <= bb) {
79 double r = (-qb + ((qs == 0) ? 1 : -1) * sqrt(bb - ac4)) / (2 * qa);
80 double x = k;
81 double y = l;
82
83 // set the circle
84 if (r >= 0.0) {
85 m_c = Point(x, y);
86 m_radius = r;
87 }
88 }
89 }
90 }
91
PointIsOn(const Point & p,double accuracy)92 bool CircleOrLine::PointIsOn(const Point& p, double accuracy)
93 {
94 if (m_is_a_line)
95 {
96 geoff_geometry::CLine line(geoff_geometry::Point(m_p0.x, m_p0.y), geoff_geometry::Point(m_p1.x, m_p1.y));
97 return (fabs(line.Dist(geoff_geometry::Point(p.x, p.y))) <= accuracy);
98 }
99 double rp = p.dist(m_c);
100 bool on = fabs(m_radius - rp) < accuracy;
101 return on;
102 }
103
LineIsOn(const Point & p0,const Point & p1,double accuracy)104 bool CircleOrLine::LineIsOn(const Point& p0, const Point& p1, double accuracy)
105 {
106 // checks the points are on the arc, to the given accuracy, and the mid point of the line.
107
108 if (this->m_is_a_line)
109 {
110 if (!PointIsOn(p0, accuracy))return false;
111 if (!PointIsOn(p1, accuracy))return false;
112
113 Point this_dir = m_p1 - m_p0;
114 this_dir.normalize();
115 Point dir = p1 - p0;
116 dir.normalize();
117 return ((dir * this_dir) >= -0.0000000001); // they are going in the same direction
118 }
119
120 if (!PointIsOn(p0, accuracy))return false;
121 if (!PointIsOn(p1, accuracy))return false;
122
123 Point mid = Point((p0 + p1) / 2);
124 if (!PointIsOn(mid, accuracy))return false;
125
126 return true;
127 }