| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /** @file | ||
| 2 | * @brief Circle shape | ||
| 3 | *//* | ||
| 4 | * Authors: | ||
| 5 | * Marco Cecchetti <mrcekets at gmail.com> | ||
| 6 | * Krzysztof Kosiński <tweenk.pl@gmail.com> | ||
| 7 | * | ||
| 8 | * Copyright 2008-2014 Authors | ||
| 9 | * | ||
| 10 | * This library is free software; you can redistribute it and/or | ||
| 11 | * modify it either under the terms of the GNU Lesser General Public | ||
| 12 | * License version 2.1 as published by the Free Software Foundation | ||
| 13 | * (the "LGPL") or, at your option, under the terms of the Mozilla | ||
| 14 | * Public License Version 1.1 (the "MPL"). If you do not alter this | ||
| 15 | * notice, a recipient may use your version of this file under either | ||
| 16 | * the MPL or the LGPL. | ||
| 17 | * | ||
| 18 | * You should have received a copy of the LGPL along with this library | ||
| 19 | * in the file COPYING-LGPL-2.1; if not, write to the Free Software | ||
| 20 | * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | ||
| 21 | * You should have received a copy of the MPL along with this library | ||
| 22 | * in the file COPYING-MPL-1.1 | ||
| 23 | * | ||
| 24 | * The contents of this file are subject to the Mozilla Public License | ||
| 25 | * Version 1.1 (the "License"); you may not use this file except in | ||
| 26 | * compliance with the License. You may obtain a copy of the License at | ||
| 27 | * http://www.mozilla.org/MPL/ | ||
| 28 | * | ||
| 29 | * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY | ||
| 30 | * OF ANY KIND, either express or implied. See the LGPL or the MPL for | ||
| 31 | * the specific language governing rights and limitations. | ||
| 32 | */ | ||
| 33 | |||
| 34 | #include <2geom/circle.h> | ||
| 35 | #include <2geom/ellipse.h> | ||
| 36 | #include <2geom/elliptical-arc.h> | ||
| 37 | #include <2geom/numeric/fitting-tool.h> | ||
| 38 | #include <2geom/numeric/fitting-model.h> | ||
| 39 | |||
| 40 | namespace Geom { | ||
| 41 | |||
| 42 | ✗ | Rect Circle::boundsFast() const | |
| 43 | { | ||
| 44 | ✗ | Point rr(_radius, _radius); | |
| 45 | ✗ | Rect bbox(_center - rr, _center + rr); | |
| 46 | ✗ | return bbox; | |
| 47 | } | ||
| 48 | |||
| 49 | 1 | void Circle::setCoefficients(Coord A, Coord B, Coord C, Coord D) | |
| 50 | { | ||
| 51 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (A == 0) { |
| 52 | ✗ | THROW_RANGEERROR("square term coefficient == 0"); | |
| 53 | } | ||
| 54 | |||
| 55 | //std::cerr << "B = " << B << " C = " << C << " D = " << D << std::endl; | ||
| 56 | |||
| 57 | 1 | Coord b = B / A; | |
| 58 | 1 | Coord c = C / A; | |
| 59 | 1 | Coord d = D / A; | |
| 60 | |||
| 61 | 1 | _center[X] = -b/2; | |
| 62 | 1 | _center[Y] = -c/2; | |
| 63 | 1 | Coord r2 = _center[X] * _center[X] + _center[Y] * _center[Y] - d; | |
| 64 | |||
| 65 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (r2 < 0) { |
| 66 | ✗ | THROW_RANGEERROR("ray^2 < 0"); | |
| 67 | } | ||
| 68 | |||
| 69 | 1 | _radius = std::sqrt(r2); | |
| 70 | 1 | } | |
| 71 | |||
| 72 | 1 | void Circle::coefficients(Coord &A, Coord &B, Coord &C, Coord &D) const | |
| 73 | { | ||
| 74 | 1 | A = 1; | |
| 75 | 1 | B = -2 * _center[X]; | |
| 76 | 1 | C = -2 * _center[Y]; | |
| 77 | 1 | D = _center[X] * _center[X] + _center[Y] * _center[Y] - _radius * _radius; | |
| 78 | 1 | } | |
| 79 | |||
| 80 | ✗ | std::vector<Coord> Circle::coefficients() const | |
| 81 | { | ||
| 82 | ✗ | std::vector<Coord> c(4); | |
| 83 | ✗ | coefficients(c[0], c[1], c[2], c[3]); | |
| 84 | ✗ | return c; | |
| 85 | ✗ | } | |
| 86 | |||
| 87 | |||
| 88 | 1 | Zoom Circle::unitCircleTransform() const | |
| 89 | { | ||
| 90 | 1 | Zoom ret(_radius, _center / _radius); | |
| 91 | 1 | return ret; | |
| 92 | } | ||
| 93 | |||
| 94 | 1 | Zoom Circle::inverseUnitCircleTransform() const | |
| 95 | { | ||
| 96 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (_radius == 0) { |
| 97 | ✗ | THROW_RANGEERROR("degenerate circle does not have an inverse unit circle transform"); | |
| 98 | } | ||
| 99 | |||
| 100 | 1 | Zoom ret(1/_radius, Translate(-_center)); | |
| 101 | 1 | return ret; | |
| 102 | } | ||
| 103 | |||
| 104 | ✗ | Point Circle::initialPoint() const | |
| 105 | { | ||
| 106 | ✗ | Point p(_center); | |
| 107 | ✗ | p[X] += _radius; | |
| 108 | ✗ | return p; | |
| 109 | } | ||
| 110 | |||
| 111 | 114 | Point Circle::pointAt(Coord t) const { | |
| 112 |
1/2✓ Branch 3 taken 114 times.
✗ Branch 4 not taken.
|
114 | return _center + Point::polar(t) * _radius; |
| 113 | } | ||
| 114 | |||
| 115 | ✗ | Coord Circle::valueAt(Coord t, Dim2 d) const { | |
| 116 | ✗ | Coord delta = (d == X ? std::cos(t) : std::sin(t)); | |
| 117 | ✗ | return _center[d] + delta * _radius; | |
| 118 | } | ||
| 119 | |||
| 120 | 13 | Coord Circle::timeAt(Point const &p) const { | |
| 121 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 13 times.
|
13 | if (_center == p) return 0; |
| 122 |
1/2✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
|
13 | return atan2(p - _center); |
| 123 | } | ||
| 124 | |||
| 125 | ✗ | Coord Circle::nearestTime(Point const &p) const { | |
| 126 | ✗ | return timeAt(p); | |
| 127 | } | ||
| 128 | |||
| 129 | ✗ | bool Circle::contains(Rect const &r) const | |
| 130 | { | ||
| 131 | ✗ | for (unsigned i = 0; i < 4; ++i) { | |
| 132 | ✗ | if (!contains(r.corner(i))) return false; | |
| 133 | } | ||
| 134 | ✗ | return true; | |
| 135 | } | ||
| 136 | |||
| 137 | 6 | bool Circle::contains(Circle const &other) const | |
| 138 | { | ||
| 139 | 6 | Coord cdist = distance(_center, other._center); | |
| 140 | 6 | Coord rdist = fabs(_radius - other._radius); | |
| 141 | 6 | return cdist <= rdist; | |
| 142 | } | ||
| 143 | |||
| 144 | 3 | bool Circle::intersects(Line const &l) const | |
| 145 | { | ||
| 146 | // http://mathworld.wolfram.com/Circle-LineIntersection.html | ||
| 147 |
1/2✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
3 | Coord dr = l.vector().length(); |
| 148 | 3 | Coord r = _radius; | |
| 149 | 3 | Coord D = cross(l.initialPoint(), l.finalPoint()); | |
| 150 | 3 | Coord delta = r*r * dr*dr - D*D; | |
| 151 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
|
3 | if (delta >= 0) return true; |
| 152 | 1 | return false; | |
| 153 | } | ||
| 154 | |||
| 155 | 6 | bool Circle::intersects(Circle const &other) const | |
| 156 | { | ||
| 157 | 6 | Coord cdist = distance(_center, other._center); | |
| 158 | 6 | Coord rsum = _radius + other._radius; | |
| 159 | 6 | return cdist <= rsum; | |
| 160 | } | ||
| 161 | |||
| 162 | |||
| 163 | 3 | std::vector<ShapeIntersection> Circle::intersect(Line const &l) const | |
| 164 | { | ||
| 165 | // http://mathworld.wolfram.com/Circle-LineIntersection.html | ||
| 166 |
1/2✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
3 | Coord dr = l.vector().length(); |
| 167 |
1/2✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
3 | Coord dx = l.vector().x(); |
| 168 |
1/2✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
3 | Coord dy = l.vector().y(); |
| 169 | 3 | Coord D = cross(l.initialPoint() - _center, l.finalPoint() - _center); | |
| 170 | 3 | Coord delta = _radius*_radius * dr*dr - D*D; | |
| 171 | |||
| 172 | 3 | std::vector<ShapeIntersection> result; | |
| 173 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
|
3 | if (delta < 0) return result; |
| 174 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
|
2 | if (delta == 0) { |
| 175 | 1 | Coord ix = (D*dy) / (dr*dr); | |
| 176 | 1 | Coord iy = (-D*dx) / (dr*dr); | |
| 177 | 1 | Point ip(ix, iy); ip += _center; | |
| 178 |
3/6✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
|
1 | result.emplace_back(timeAt(ip), l.timeAt(ip), ip); |
| 179 | 1 | return result; | |
| 180 | } | ||
| 181 | |||
| 182 | 1 | Coord sqrt_delta = std::sqrt(delta); | |
| 183 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | Coord signmod = dy < 0 ? -1 : 1; |
| 184 | |||
| 185 | 1 | Coord i1x = (D*dy + signmod * dx * sqrt_delta) / (dr*dr); | |
| 186 | 1 | Coord i1y = (-D*dx + fabs(dy) * sqrt_delta) / (dr*dr); | |
| 187 | 1 | Point i1p(i1x, i1y); i1p += _center; | |
| 188 | |||
| 189 | 1 | Coord i2x = (D*dy - signmod * dx * sqrt_delta) / (dr*dr); | |
| 190 | 1 | Coord i2y = (-D*dx - fabs(dy) * sqrt_delta) / (dr*dr); | |
| 191 | 1 | Point i2p(i2x, i2y); i2p += _center; | |
| 192 | |||
| 193 |
3/6✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
|
1 | result.emplace_back(timeAt(i1p), l.timeAt(i1p), i1p); |
| 194 |
3/6✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
|
1 | result.emplace_back(timeAt(i2p), l.timeAt(i2p), i2p); |
| 195 | 1 | return result; | |
| 196 | ✗ | } | |
| 197 | |||
| 198 | ✗ | std::vector<ShapeIntersection> Circle::intersect(LineSegment const &l) const | |
| 199 | { | ||
| 200 | ✗ | std::vector<ShapeIntersection> result = intersect(Line(l)); | |
| 201 | ✗ | filter_line_segment_intersections(result); | |
| 202 | ✗ | return result; | |
| 203 | ✗ | } | |
| 204 | |||
| 205 | 6 | std::vector<ShapeIntersection> Circle::intersect(Circle const &other) const | |
| 206 | { | ||
| 207 | 6 | std::vector<ShapeIntersection> result; | |
| 208 | |||
| 209 |
2/4✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
|
6 | if (*this == other) { |
| 210 | ✗ | THROW_INFINITESOLUTIONS(); | |
| 211 | } | ||
| 212 |
2/4✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
|
6 | if (contains(other)) return result; |
| 213 |
3/4✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 3 times.
|
6 | if (!intersects(other)) return result; |
| 214 | |||
| 215 | // See e.g. http://mathworld.wolfram.com/Circle-CircleIntersection.html | ||
| 216 | // Basically, we figure out where is the third point of a triangle | ||
| 217 | // with two points in the centers and with edge lengths equal to radii | ||
| 218 | 3 | Point cv = other._center - _center; | |
| 219 | 3 | Coord d = cv.length(); | |
| 220 | 3 | Coord R = radius(), r = other.radius(); | |
| 221 | |||
| 222 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
|
3 | if (d == R + r) { |
| 223 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
1 | Point px = lerp(R / d, _center, other._center); |
| 224 |
2/4✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
1 | Coord T = timeAt(px), t = other.timeAt(px); |
| 225 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | result.emplace_back(T, t, px); |
| 226 | 1 | return result; | |
| 227 | } | ||
| 228 | |||
| 229 | // q is the distance along the line between centers to the perpendicular line | ||
| 230 | // that goes through both intersections. | ||
| 231 | 2 | Coord q = (d*d - r*r + R*R) / (2*d); | |
| 232 |
1/2✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
2 | Point qp = lerp(q/d, _center, other._center); |
| 233 | |||
| 234 | // The triangle given by the points: | ||
| 235 | // _center, qp, intersection | ||
| 236 | // is a right triangle. Determine the distance between qp and intersection | ||
| 237 | // using the Pythagorean theorem. | ||
| 238 | 2 | Coord h = std::sqrt(R*R - q*q); | |
| 239 | 2 | Point qd = (h/d) * cv.cw(); | |
| 240 | |||
| 241 | // now compute the intersection points | ||
| 242 | 2 | Point x1 = qp + qd; | |
| 243 | 2 | Point x2 = qp - qd; | |
| 244 | |||
| 245 |
3/6✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
|
2 | result.emplace_back(timeAt(x1), other.timeAt(x1), x1); |
| 246 |
3/6✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
|
2 | result.emplace_back(timeAt(x2), other.timeAt(x2), x2); |
| 247 | 2 | return result; | |
| 248 | ✗ | } | |
| 249 | |||
| 250 | /** | ||
| 251 | @param inner a point whose angle with the circle center is inside the angle that the arc spans | ||
| 252 | */ | ||
| 253 | EllipticalArc * | ||
| 254 | ✗ | Circle::arc(Point const& initial, Point const& inner, Point const& final) const | |
| 255 | { | ||
| 256 | // TODO native implementation! | ||
| 257 | ✗ | Ellipse e(_center[X], _center[Y], _radius, _radius, 0); | |
| 258 | ✗ | return e.arc(initial, inner, final); | |
| 259 | } | ||
| 260 | |||
| 261 | 9 | bool Circle::operator==(Circle const &other) const | |
| 262 | { | ||
| 263 |
2/2✓ Branch 1 taken 8 times.
✓ Branch 2 taken 1 times.
|
9 | if (_center != other._center) return false; |
| 264 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (_radius != other._radius) return false; |
| 265 | 1 | return true; | |
| 266 | } | ||
| 267 | |||
| 268 | ✗ | D2<SBasis> Circle::toSBasis() const | |
| 269 | { | ||
| 270 | ✗ | D2<SBasis> B; | |
| 271 | ✗ | Linear bo = Linear(0, 2 * M_PI); | |
| 272 | |||
| 273 | ✗ | B[0] = cos(bo,4); | |
| 274 | ✗ | B[1] = sin(bo,4); | |
| 275 | |||
| 276 | ✗ | B = B * _radius + _center; | |
| 277 | |||
| 278 | ✗ | return B; | |
| 279 | ✗ | } | |
| 280 | |||
| 281 | |||
| 282 | ✗ | void Circle::fit(std::vector<Point> const& points) | |
| 283 | { | ||
| 284 | ✗ | size_t sz = points.size(); | |
| 285 | ✗ | if (sz < 2) { | |
| 286 | ✗ | THROW_RANGEERROR("fitting error: too few points passed"); | |
| 287 | } | ||
| 288 | ✗ | if (sz == 2) { | |
| 289 | ✗ | _center = points[0] * 0.5 + points[1] * 0.5; | |
| 290 | ✗ | _radius = distance(points[0], points[1]) / 2; | |
| 291 | ✗ | return; | |
| 292 | } | ||
| 293 | |||
| 294 | ✗ | NL::LFMCircle model; | |
| 295 | ✗ | NL::least_squeares_fitter<NL::LFMCircle> fitter(model, sz); | |
| 296 | |||
| 297 | ✗ | for (size_t i = 0; i < sz; ++i) { | |
| 298 | ✗ | fitter.append(points[i]); | |
| 299 | } | ||
| 300 | ✗ | fitter.update(); | |
| 301 | |||
| 302 | ✗ | NL::Vector z(sz, 0.0); | |
| 303 | ✗ | model.instance(*this, fitter.result(z)); | |
| 304 | ✗ | } | |
| 305 | |||
| 306 | |||
| 307 | 8 | bool are_near(Circle const &a, Circle const &b, Coord eps) | |
| 308 | { | ||
| 309 | // to check whether no point on a is further than eps from b, | ||
| 310 | // we check two things: | ||
| 311 | // 1. if radii differ by more than eps, there is definitely a point that fails | ||
| 312 | // 2. if they differ by less, we check the centers. They have to be closer | ||
| 313 | // together if the radius differs, since the maximum distance will be | ||
| 314 | // equal to sum of radius difference and distance between centers. | ||
| 315 |
2/2✓ Branch 3 taken 3 times.
✓ Branch 4 taken 5 times.
|
8 | if (!are_near(a.radius(), b.radius(), eps)) return false; |
| 316 | 5 | Coord adjusted_eps = eps - fabs(a.radius() - b.radius()); | |
| 317 |
1/2✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
|
5 | return are_near(a.center(), b.center(), adjusted_eps); |
| 318 | } | ||
| 319 | |||
| 320 | ✗ | std::ostream &operator<<(std::ostream &out, Circle const &c) | |
| 321 | { | ||
| 322 | ✗ | out << "Circle(" << c.center() << ", " << format_coord_nice(c.radius()) << ")"; | |
| 323 | ✗ | return out; | |
| 324 | } | ||
| 325 | |||
| 326 | } // end namespace Geom | ||
| 327 | |||
| 328 | /* | ||
| 329 | Local Variables: | ||
| 330 | mode:c++ | ||
| 331 | c-file-style:"stroustrup" | ||
| 332 | c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) | ||
| 333 | indent-tabs-mode:nil | ||
| 334 | fill-column:99 | ||
| 335 | End: | ||
| 336 | */ | ||
| 337 | // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : | ||
| 338 |