| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /** @file | ||
| 2 | * @brief Basic intersection routines | ||
| 3 | *//* | ||
| 4 | * Authors: | ||
| 5 | * Nathan Hurst <njh@njhurst.com> | ||
| 6 | * Marco Cecchetti <mrcekets at gmail.com> | ||
| 7 | * Jean-François Barraud <jf.barraud@gmail.com> | ||
| 8 | * | ||
| 9 | * Copyright 2008-2009 Authors | ||
| 10 | * | ||
| 11 | * This library is free software; you can redistribute it and/or | ||
| 12 | * modify it either under the terms of the GNU Lesser General Public | ||
| 13 | * License version 2.1 as published by the Free Software Foundation | ||
| 14 | * (the "LGPL") or, at your option, under the terms of the Mozilla | ||
| 15 | * Public License Version 1.1 (the "MPL"). If you do not alter this | ||
| 16 | * notice, a recipient may use your version of this file under either | ||
| 17 | * the MPL or the LGPL. | ||
| 18 | * | ||
| 19 | * You should have received a copy of the LGPL along with this library | ||
| 20 | * in the file COPYING-LGPL-2.1; if not, write to the Free Software | ||
| 21 | * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | ||
| 22 | * You should have received a copy of the MPL along with this library | ||
| 23 | * in the file COPYING-MPL-1.1 | ||
| 24 | * | ||
| 25 | * The contents of this file are subject to the Mozilla Public License | ||
| 26 | * Version 1.1 (the "License"); you may not use this file except in | ||
| 27 | * compliance with the License. You may obtain a copy of the License at | ||
| 28 | * http://www.mozilla.org/MPL/ | ||
| 29 | * | ||
| 30 | * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY | ||
| 31 | * OF ANY KIND, either express or implied. See the LGPL or the MPL for | ||
| 32 | * the specific language governing rights and limitations. | ||
| 33 | * | ||
| 34 | */ | ||
| 35 | |||
| 36 | #include <2geom/basic-intersection.h> | ||
| 37 | #include <2geom/sbasis-to-bezier.h> | ||
| 38 | #include <2geom/exception.h> | ||
| 39 | |||
| 40 | #ifdef HAVE_GSL | ||
| 41 | #include <gsl/gsl_vector.h> | ||
| 42 | #include <gsl/gsl_multiroots.h> | ||
| 43 | #endif | ||
| 44 | |||
| 45 | using std::vector; | ||
| 46 | namespace Geom { | ||
| 47 | |||
| 48 | //#ifdef USE_RECURSIVE_INTERSECTOR | ||
| 49 | |||
| 50 | // void find_intersections(std::vector<std::pair<double, double> > &xs, | ||
| 51 | // D2<SBasis> const & A, | ||
| 52 | // D2<SBasis> const & B) { | ||
| 53 | // vector<Point> BezA, BezB; | ||
| 54 | // sbasis_to_bezier(BezA, A); | ||
| 55 | // sbasis_to_bezier(BezB, B); | ||
| 56 | |||
| 57 | // xs.clear(); | ||
| 58 | |||
| 59 | // find_intersections_bezier_recursive(xs, BezA, BezB); | ||
| 60 | // } | ||
| 61 | // void find_intersections(std::vector< std::pair<double, double> > & xs, | ||
| 62 | // std::vector<Point> const& A, | ||
| 63 | // std::vector<Point> const& B, | ||
| 64 | // double precision){ | ||
| 65 | // find_intersections_bezier_recursive(xs, A, B, precision); | ||
| 66 | // } | ||
| 67 | |||
| 68 | //#else | ||
| 69 | |||
| 70 | namespace detail{ namespace bezier_clipping { | ||
| 71 | void portion(std::vector<Point> &B, Interval const &I); | ||
| 72 | void derivative(std::vector<Point> &D, std::vector<Point> const &B); | ||
| 73 | }; }; | ||
| 74 | |||
| 75 | 58 | void find_intersections(std::vector<std::pair<double, double> > &xs, | |
| 76 | D2<Bezier> const & A, | ||
| 77 | D2<Bezier> const & B, | ||
| 78 | double precision) | ||
| 79 | { | ||
| 80 |
3/6✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 58 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 58 times.
✗ Branch 10 not taken.
|
58 | find_intersections_bezier_clipping(xs, bezier_points(A), bezier_points(B), precision); |
| 81 | 58 | } | |
| 82 | |||
| 83 | ✗ | void find_intersections(std::vector<std::pair<double, double> > &xs, | |
| 84 | D2<SBasis> const & A, | ||
| 85 | D2<SBasis> const & B, | ||
| 86 | double precision) | ||
| 87 | { | ||
| 88 | ✗ | vector<Point> BezA, BezB; | |
| 89 | ✗ | sbasis_to_bezier(BezA, A); | |
| 90 | ✗ | sbasis_to_bezier(BezB, B); | |
| 91 | |||
| 92 | ✗ | find_intersections_bezier_clipping(xs, BezA, BezB, precision); | |
| 93 | ✗ | } | |
| 94 | |||
| 95 | ✗ | void find_intersections(std::vector< std::pair<double, double> > & xs, | |
| 96 | std::vector<Point> const& A, | ||
| 97 | std::vector<Point> const& B, | ||
| 98 | double precision) | ||
| 99 | { | ||
| 100 | ✗ | find_intersections_bezier_clipping(xs, A, B, precision); | |
| 101 | ✗ | } | |
| 102 | |||
| 103 | //#endif | ||
| 104 | |||
| 105 | /* | ||
| 106 | * split the curve at the midpoint, returning an array with the two parts | ||
| 107 | * Temporary storage is minimized by using part of the storage for the result | ||
| 108 | * to hold an intermediate value until it is no longer needed. | ||
| 109 | */ | ||
| 110 | // TODO replace with Bezier method | ||
| 111 | ✗ | void split(vector<Point> const &p, double t, | |
| 112 | vector<Point> &left, vector<Point> &right) { | ||
| 113 | ✗ | const unsigned sz = p.size(); | |
| 114 | //Geom::Point Vtemp[sz][sz]; | ||
| 115 | ✗ | vector<vector<Point> > Vtemp(sz); | |
| 116 | ✗ | for ( size_t i = 0; i < sz; ++i ) | |
| 117 | ✗ | Vtemp[i].reserve(sz); | |
| 118 | |||
| 119 | /* Copy control points */ | ||
| 120 | ✗ | std::copy(p.begin(), p.end(), Vtemp[0].begin()); | |
| 121 | |||
| 122 | /* Triangle computation */ | ||
| 123 | ✗ | for (unsigned i = 1; i < sz; i++) { | |
| 124 | ✗ | for (unsigned j = 0; j < sz - i; j++) { | |
| 125 | ✗ | Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]); | |
| 126 | } | ||
| 127 | } | ||
| 128 | |||
| 129 | ✗ | left.resize(sz); | |
| 130 | ✗ | right.resize(sz); | |
| 131 | ✗ | for (unsigned j = 0; j < sz; j++) | |
| 132 | ✗ | left[j] = Vtemp[j][0]; | |
| 133 | ✗ | for (unsigned j = 0; j < sz; j++) | |
| 134 | ✗ | right[j] = Vtemp[sz-1-j][j]; | |
| 135 | ✗ | } | |
| 136 | |||
| 137 | |||
| 138 | |||
| 139 | ✗ | void find_self_intersections(std::vector<std::pair<double, double> > &xs, | |
| 140 | D2<Bezier> const &A, | ||
| 141 | double precision) | ||
| 142 | { | ||
| 143 | ✗ | std::vector<double> dr = derivative(A[X]).roots(); | |
| 144 | { | ||
| 145 | ✗ | std::vector<double> dyr = derivative(A[Y]).roots(); | |
| 146 | ✗ | dr.insert(dr.begin(), dyr.begin(), dyr.end()); | |
| 147 | ✗ | } | |
| 148 | ✗ | dr.push_back(0); | |
| 149 | ✗ | dr.push_back(1); | |
| 150 | // We want to be sure that we have no empty segments | ||
| 151 | ✗ | std::sort(dr.begin(), dr.end()); | |
| 152 | ✗ | std::vector<double>::iterator new_end = std::unique(dr.begin(), dr.end()); | |
| 153 | ✗ | dr.resize( new_end - dr.begin() ); | |
| 154 | |||
| 155 | ✗ | std::vector< D2<Bezier> > pieces; | |
| 156 | ✗ | for (unsigned i = 0; i < dr.size() - 1; ++i) { | |
| 157 | ✗ | pieces.push_back(portion(A, dr[i], dr[i+1])); | |
| 158 | } | ||
| 159 | /*{ | ||
| 160 | vector<Point> l, r, in = A; | ||
| 161 | for(unsigned i = 0; i < dr.size()-1; i++) { | ||
| 162 | split(in, (dr[i+1]-dr[i]) / (1 - dr[i]), l, r); | ||
| 163 | pieces.push_back(l); | ||
| 164 | in = r; | ||
| 165 | } | ||
| 166 | }*/ | ||
| 167 | |||
| 168 | ✗ | for(unsigned i = 0; i < dr.size()-1; i++) { | |
| 169 | ✗ | for(unsigned j = i+1; j < dr.size()-1; j++) { | |
| 170 | ✗ | std::vector<std::pair<double, double> > section; | |
| 171 | |||
| 172 | ✗ | find_intersections(section, pieces[i], pieces[j], precision); | |
| 173 | ✗ | for(auto & k : section) { | |
| 174 | ✗ | double l = k.first; | |
| 175 | ✗ | double r = k.second; | |
| 176 | // XXX: This condition will prune out false positives, but it might create some false negatives. Todo: Confirm it is correct. | ||
| 177 | ✗ | if(j == i+1) | |
| 178 | //if((l == 1) && (r == 0)) | ||
| 179 | ✗ | if( ( l > precision ) && (r < precision) )//FIXME: what precision should be used here??? | |
| 180 | ✗ | continue; | |
| 181 | ✗ | xs.emplace_back((1-l)*dr[i] + l*dr[i+1], | |
| 182 | ✗ | (1-r)*dr[j] + r*dr[j+1]); | |
| 183 | } | ||
| 184 | ✗ | } | |
| 185 | } | ||
| 186 | |||
| 187 | // Because i is in order, xs should be roughly already in order? | ||
| 188 | //sort(xs.begin(), xs.end()); | ||
| 189 | //unique(xs.begin(), xs.end()); | ||
| 190 | ✗ | } | |
| 191 | |||
| 192 | ✗ | void find_self_intersections(std::vector<std::pair<double, double> > &xs, | |
| 193 | D2<SBasis> const &A, | ||
| 194 | double precision) | ||
| 195 | { | ||
| 196 | ✗ | D2<Bezier> in; | |
| 197 | ✗ | sbasis_to_bezier(in, A); | |
| 198 | ✗ | find_self_intersections(xs, in, precision); | |
| 199 | ✗ | } | |
| 200 | |||
| 201 | |||
| 202 | ✗ | void subdivide(D2<Bezier> const &a, | |
| 203 | D2<Bezier> const &b, | ||
| 204 | std::vector< std::pair<double, double> > const &xs, | ||
| 205 | std::vector< D2<Bezier> > &av, | ||
| 206 | std::vector< D2<Bezier> > &bv) | ||
| 207 | { | ||
| 208 | ✗ | if (xs.empty()) { | |
| 209 | ✗ | av.push_back(a); | |
| 210 | ✗ | bv.push_back(b); | |
| 211 | ✗ | return; | |
| 212 | } | ||
| 213 | |||
| 214 | ✗ | std::pair<double, double> prev = std::make_pair(0., 0.); | |
| 215 | ✗ | for (const auto & x : xs) { | |
| 216 | ✗ | av.push_back(portion(a, prev.first, x.first)); | |
| 217 | ✗ | bv.push_back(portion(b, prev.second, x.second)); | |
| 218 | ✗ | av.back()[X].at0() = bv.back()[X].at0() = lerp(0.5, av.back()[X].at0(), bv.back()[X].at0()); | |
| 219 | ✗ | av.back()[X].at1() = bv.back()[X].at1() = lerp(0.5, av.back()[X].at1(), bv.back()[X].at1()); | |
| 220 | ✗ | av.back()[Y].at0() = bv.back()[Y].at0() = lerp(0.5, av.back()[Y].at0(), bv.back()[Y].at0()); | |
| 221 | ✗ | av.back()[Y].at1() = bv.back()[Y].at1() = lerp(0.5, av.back()[Y].at1(), bv.back()[Y].at1()); | |
| 222 | ✗ | prev = x; | |
| 223 | } | ||
| 224 | ✗ | av.push_back(portion(a, prev.first, 1)); | |
| 225 | ✗ | bv.push_back(portion(b, prev.second, 1)); | |
| 226 | ✗ | av.back()[X].at0() = bv.back()[X].at0() = lerp(0.5, av.back()[X].at0(), bv.back()[X].at0()); | |
| 227 | ✗ | av.back()[X].at1() = bv.back()[X].at1() = lerp(0.5, av.back()[X].at1(), bv.back()[X].at1()); | |
| 228 | ✗ | av.back()[Y].at0() = bv.back()[Y].at0() = lerp(0.5, av.back()[Y].at0(), bv.back()[Y].at0()); | |
| 229 | ✗ | av.back()[Y].at1() = bv.back()[Y].at1() = lerp(0.5, av.back()[Y].at1(), bv.back()[Y].at1()); | |
| 230 | } | ||
| 231 | |||
| 232 | #ifdef HAVE_GSL | ||
| 233 | #include <gsl/gsl_multiroots.h> | ||
| 234 | |||
| 235 | struct rparams | ||
| 236 | { | ||
| 237 | D2<SBasis> const &A; | ||
| 238 | D2<SBasis> const &B; | ||
| 239 | }; | ||
| 240 | |||
| 241 | static int | ||
| 242 | ✗ | intersect_polish_f (const gsl_vector * x, void *params, | |
| 243 | gsl_vector * f) | ||
| 244 | { | ||
| 245 | ✗ | const double x0 = gsl_vector_get (x, 0); | |
| 246 | ✗ | const double x1 = gsl_vector_get (x, 1); | |
| 247 | |||
| 248 | ✗ | Geom::Point dx = ((struct rparams *) params)->A(x0) - | |
| 249 | ✗ | ((struct rparams *) params)->B(x1); | |
| 250 | |||
| 251 | ✗ | gsl_vector_set (f, 0, dx[0]); | |
| 252 | ✗ | gsl_vector_set (f, 1, dx[1]); | |
| 253 | |||
| 254 | ✗ | return GSL_SUCCESS; | |
| 255 | } | ||
| 256 | #endif | ||
| 257 | |||
| 258 | union dbl_64{ | ||
| 259 | long long i64; | ||
| 260 | double d64; | ||
| 261 | }; | ||
| 262 | |||
| 263 | ✗ | static double EpsilonBy(double value, int eps) | |
| 264 | { | ||
| 265 | dbl_64 s; | ||
| 266 | ✗ | s.d64 = value; | |
| 267 | ✗ | s.i64 += eps; | |
| 268 | ✗ | return s.d64; | |
| 269 | } | ||
| 270 | |||
| 271 | |||
| 272 | ✗ | static void intersect_polish_root (D2<SBasis> const &A, double &s, | |
| 273 | D2<SBasis> const &B, double &t) { | ||
| 274 | #ifdef HAVE_GSL | ||
| 275 | const gsl_multiroot_fsolver_type *T; | ||
| 276 | gsl_multiroot_fsolver *sol; | ||
| 277 | |||
| 278 | int status; | ||
| 279 | ✗ | size_t iter = 0; | |
| 280 | #endif | ||
| 281 | ✗ | std::vector<Point> as, bs; | |
| 282 | ✗ | as = A.valueAndDerivatives(s, 2); | |
| 283 | ✗ | bs = B.valueAndDerivatives(t, 2); | |
| 284 | ✗ | Point F = as[0] - bs[0]; | |
| 285 | ✗ | double best = dot(F, F); | |
| 286 | |||
| 287 | ✗ | for(int i = 0; i < 4; i++) { | |
| 288 | |||
| 289 | /** | ||
| 290 | we want to solve | ||
| 291 | J*(x1 - x0) = f(x0) | ||
| 292 | |||
| 293 | |dA(s)[0] -dB(t)[0]| (X1 - X0) = A(s) - B(t) | ||
| 294 | |dA(s)[1] -dB(t)[1]| | ||
| 295 | **/ | ||
| 296 | |||
| 297 | // We're using the standard transformation matricies, which is numerically rather poor. Much better to solve the equation using elimination. | ||
| 298 | |||
| 299 | ✗ | Affine jack(as[1][0], as[1][1], | |
| 300 | ✗ | -bs[1][0], -bs[1][1], | |
| 301 | ✗ | 0, 0); | |
| 302 | ✗ | Point soln = (F)*jack.inverse(); | |
| 303 | ✗ | double ns = s - soln[0]; | |
| 304 | ✗ | double nt = t - soln[1]; | |
| 305 | |||
| 306 | ✗ | as = A.valueAndDerivatives(ns, 2); | |
| 307 | ✗ | bs = B.valueAndDerivatives(nt, 2); | |
| 308 | ✗ | F = as[0] - bs[0]; | |
| 309 | ✗ | double trial = dot(F, F); | |
| 310 | ✗ | if (trial > best*0.1) {// we have standards, you know | |
| 311 | // At this point we could do a line search | ||
| 312 | ✗ | break; | |
| 313 | } | ||
| 314 | ✗ | best = trial; | |
| 315 | ✗ | s = ns; | |
| 316 | ✗ | t = nt; | |
| 317 | } | ||
| 318 | |||
| 319 | #ifdef HAVE_GSL | ||
| 320 | ✗ | const size_t n = 2; | |
| 321 | ✗ | struct rparams p = {A, B}; | |
| 322 | ✗ | gsl_multiroot_function f = {&intersect_polish_f, n, &p}; | |
| 323 | |||
| 324 | ✗ | double x_init[2] = {s, t}; | |
| 325 | ✗ | gsl_vector *x = gsl_vector_alloc (n); | |
| 326 | |||
| 327 | ✗ | gsl_vector_set (x, 0, x_init[0]); | |
| 328 | ✗ | gsl_vector_set (x, 1, x_init[1]); | |
| 329 | |||
| 330 | ✗ | T = gsl_multiroot_fsolver_hybrids; | |
| 331 | ✗ | sol = gsl_multiroot_fsolver_alloc (T, 2); | |
| 332 | ✗ | gsl_multiroot_fsolver_set (sol, &f, x); | |
| 333 | |||
| 334 | do | ||
| 335 | { | ||
| 336 | ✗ | iter++; | |
| 337 | ✗ | status = gsl_multiroot_fsolver_iterate (sol); | |
| 338 | |||
| 339 | ✗ | if (status) /* check if solver is stuck */ | |
| 340 | ✗ | break; | |
| 341 | |||
| 342 | status = | ||
| 343 | ✗ | gsl_multiroot_test_residual (sol->f, 1e-12); | |
| 344 | } | ||
| 345 | ✗ | while (status == GSL_CONTINUE && iter < 1000); | |
| 346 | |||
| 347 | ✗ | s = gsl_vector_get (sol->x, 0); | |
| 348 | ✗ | t = gsl_vector_get (sol->x, 1); | |
| 349 | |||
| 350 | ✗ | gsl_multiroot_fsolver_free (sol); | |
| 351 | ✗ | gsl_vector_free (x); | |
| 352 | #endif | ||
| 353 | |||
| 354 | { | ||
| 355 | // This code does a neighbourhood search for minor improvements. | ||
| 356 | ✗ | double best_v = L1(A(s) - B(t)); | |
| 357 | //std::cout << "------\n" << best_v << std::endl; | ||
| 358 | ✗ | Point best(s,t); | |
| 359 | while (true) { | ||
| 360 | ✗ | Point trial = best; | |
| 361 | ✗ | double trial_v = best_v; | |
| 362 | ✗ | for(int nsi = -1; nsi < 2; nsi++) { | |
| 363 | ✗ | for(int nti = -1; nti < 2; nti++) { | |
| 364 | ✗ | Point n(EpsilonBy(best[0], nsi), | |
| 365 | ✗ | EpsilonBy(best[1], nti)); | |
| 366 | ✗ | double c = L1(A(n[0]) - B(n[1])); | |
| 367 | //std::cout << c << "; "; | ||
| 368 | ✗ | if (c < trial_v) { | |
| 369 | ✗ | trial = n; | |
| 370 | ✗ | trial_v = c; | |
| 371 | } | ||
| 372 | } | ||
| 373 | } | ||
| 374 | ✗ | if(trial == best) { | |
| 375 | //std::cout << "\n" << s << " -> " << s - best[0] << std::endl; | ||
| 376 | //std::cout << t << " -> " << t - best[1] << std::endl; | ||
| 377 | //std::cout << best_v << std::endl; | ||
| 378 | ✗ | s = best[0]; | |
| 379 | ✗ | t = best[1]; | |
| 380 | ✗ | return; | |
| 381 | } else { | ||
| 382 | ✗ | best = trial; | |
| 383 | ✗ | best_v = trial_v; | |
| 384 | } | ||
| 385 | ✗ | } | |
| 386 | } | ||
| 387 | ✗ | } | |
| 388 | |||
| 389 | |||
| 390 | ✗ | void polish_intersections(std::vector<std::pair<double, double> > &xs, | |
| 391 | D2<SBasis> const &A, D2<SBasis> const &B) | ||
| 392 | { | ||
| 393 | ✗ | for(auto & x : xs) | |
| 394 | ✗ | intersect_polish_root(A, x.first, | |
| 395 | ✗ | B, x.second); | |
| 396 | ✗ | } | |
| 397 | |||
| 398 | /** | ||
| 399 | * Compute the Hausdorf distance from A to B only. | ||
| 400 | */ | ||
| 401 | ✗ | double hausdorfl(D2<SBasis>& A, D2<SBasis> const& B, | |
| 402 | double m_precision, | ||
| 403 | double *a_t, double* b_t) { | ||
| 404 | ✗ | std::vector< std::pair<double, double> > xs; | |
| 405 | ✗ | std::vector<Point> Az, Bz; | |
| 406 | ✗ | sbasis_to_bezier (Az, A); | |
| 407 | ✗ | sbasis_to_bezier (Bz, B); | |
| 408 | ✗ | find_collinear_normal(xs, Az, Bz, m_precision); | |
| 409 | ✗ | double h_dist = 0, h_a_t = 0, h_b_t = 0; | |
| 410 | ✗ | double dist = 0; | |
| 411 | ✗ | Point Ax = A.at0(); | |
| 412 | ✗ | double t = Geom::nearest_time(Ax, B); | |
| 413 | ✗ | dist = Geom::distance(Ax, B(t)); | |
| 414 | ✗ | if (dist > h_dist) { | |
| 415 | ✗ | h_a_t = 0; | |
| 416 | ✗ | h_b_t = t; | |
| 417 | ✗ | h_dist = dist; | |
| 418 | } | ||
| 419 | ✗ | Ax = A.at1(); | |
| 420 | ✗ | t = Geom::nearest_time(Ax, B); | |
| 421 | ✗ | dist = Geom::distance(Ax, B(t)); | |
| 422 | ✗ | if (dist > h_dist) { | |
| 423 | ✗ | h_a_t = 1; | |
| 424 | ✗ | h_b_t = t; | |
| 425 | ✗ | h_dist = dist; | |
| 426 | } | ||
| 427 | ✗ | for (auto & x : xs) | |
| 428 | { | ||
| 429 | ✗ | Point At = A(x.first); | |
| 430 | ✗ | Point Bu = B(x.second); | |
| 431 | ✗ | double distAtBu = Geom::distance(At, Bu); | |
| 432 | ✗ | t = Geom::nearest_time(At, B); | |
| 433 | ✗ | dist = Geom::distance(At, B(t)); | |
| 434 | //FIXME: we might miss it due to floating point precision... | ||
| 435 | ✗ | if (dist >= distAtBu-.1 && distAtBu > h_dist) { | |
| 436 | ✗ | h_a_t = x.first; | |
| 437 | ✗ | h_b_t = x.second; | |
| 438 | ✗ | h_dist = distAtBu; | |
| 439 | } | ||
| 440 | |||
| 441 | } | ||
| 442 | ✗ | if(a_t) *a_t = h_a_t; | |
| 443 | ✗ | if(b_t) *b_t = h_b_t; | |
| 444 | |||
| 445 | ✗ | return h_dist; | |
| 446 | ✗ | } | |
| 447 | |||
| 448 | /** | ||
| 449 | * Compute the symmetric Hausdorf distance. | ||
| 450 | */ | ||
| 451 | ✗ | double hausdorf(D2<SBasis>& A, D2<SBasis> const& B, | |
| 452 | double m_precision, | ||
| 453 | double *a_t, double* b_t) { | ||
| 454 | ✗ | double h_dist = hausdorfl(A, B, m_precision, a_t, b_t); | |
| 455 | |||
| 456 | ✗ | double dist = 0; | |
| 457 | ✗ | Point Bx = B.at0(); | |
| 458 | ✗ | double t = Geom::nearest_time(Bx, A); | |
| 459 | ✗ | dist = Geom::distance(Bx, A(t)); | |
| 460 | ✗ | if (dist > h_dist) { | |
| 461 | ✗ | if(a_t) *a_t = t; | |
| 462 | ✗ | if(b_t) *b_t = 0; | |
| 463 | ✗ | h_dist = dist; | |
| 464 | } | ||
| 465 | ✗ | Bx = B.at1(); | |
| 466 | ✗ | t = Geom::nearest_time(Bx, A); | |
| 467 | ✗ | dist = Geom::distance(Bx, A(t)); | |
| 468 | ✗ | if (dist > h_dist) { | |
| 469 | ✗ | if(a_t) *a_t = t; | |
| 470 | ✗ | if(b_t) *b_t = 1; | |
| 471 | ✗ | h_dist = dist; | |
| 472 | } | ||
| 473 | |||
| 474 | ✗ | return h_dist; | |
| 475 | } | ||
| 476 | |||
| 477 | 27 | bool non_collinear_segments_intersect(const Point &A, const Point &B, const Point &C, const Point &D) | |
| 478 | { | ||
| 479 |
6/10✓ Branch 10 taken 15 times.
✓ Branch 11 taken 12 times.
✓ Branch 16 taken 27 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 27 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 27 times.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 27 times.
|
141 | return cross(D - C, A - C) * cross(D - C, B - C) < 0 && // |
| 480 |
10/10✓ Branch 6 taken 11 times.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 15 times.
✓ Branch 9 taken 12 times.
✓ Branch 11 taken 15 times.
✓ Branch 12 taken 12 times.
✓ Branch 14 taken 15 times.
✓ Branch 15 taken 12 times.
✓ Branch 17 taken 15 times.
✓ Branch 18 taken 12 times.
|
87 | cross(B - A, C - A) * cross(B - A, D - A) < 0; |
| 481 | } | ||
| 482 | }; | ||
| 483 | |||
| 484 | /* | ||
| 485 | Local Variables: | ||
| 486 | mode:c++ | ||
| 487 | c-file-style:"stroustrup" | ||
| 488 | c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) | ||
| 489 | indent-tabs-mode:nil | ||
| 490 | fill-column:99 | ||
| 491 | End: | ||
| 492 | */ | ||
| 493 | // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : | ||
| 494 |