1 /*
2 * Authors:
3 * Nathan Hurst <njh@njhurst.com
4 *
5 * Copyright 2009 authors
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it either under the terms of the GNU Lesser General Public
9 * License version 2.1 as published by the Free Software Foundation
10 * (the "LGPL") or, at your option, under the terms of the Mozilla
11 * Public License Version 1.1 (the "MPL"). If you do not alter this
12 * notice, a recipient may use your version of this file under either
13 * the MPL or the LGPL.
14 *
15 * You should have received a copy of the LGPL along with this library
16 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 * You should have received a copy of the MPL along with this library
19 * in the file COPYING-MPL-1.1
20 *
21 * The contents of this file are subject to the Mozilla Public License
22 * Version 1.1 (the "License"); you may not use this file except in
23 * compliance with the License. You may obtain a copy of the License at
24 * http://www.mozilla.org/MPL/
25 *
26 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
27 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
28 * the specific language governing rights and limitations.
29 */
30
31
32 #include <2geom/conicsec.h>
33 #include <2geom/conic_section_clipper.h>
34 #include <2geom/numeric/fitting-tool.h>
35 #include <2geom/numeric/fitting-model.h>
36
37
38 // File: convert.h
39 #include <utility>
40 #include <sstream>
41 #include <stdexcept>
42 #include <optional>
43
44 namespace Geom
45 {
46
intersection(Line l,Rect r)47 LineSegment intersection(Line l, Rect r) {
48 std::optional<LineSegment> seg = l.clip(r);
49 if (seg) {
50 return *seg;
51 } else {
52 return LineSegment(Point(0,0), Point(0,0));
53 }
54 }
55
det(Point a,Point b)56 static double det(Point a, Point b) {
57 return a[0]*b[1] - a[1]*b[0];
58 }
59
60 template <typename T>
det(T a,T b,T c,T d)61 static T det(T a, T b, T c, T d) {
62 return a*d - b*c;
63 }
64
65 template <typename T>
det(T M[2][2])66 static T det(T M[2][2]) {
67 return M[0][0]*M[1][1] - M[1][0]*M[0][1];
68 }
69
70 template <typename T>
det3(T M[3][3])71 static T det3(T M[3][3]) {
72 return ( M[0][0] * det(M[1][1], M[1][2],
73 M[2][1], M[2][2])
74 -M[1][0] * det(M[0][1], M[0][2],
75 M[2][1], M[2][2])
76 +M[2][0] * det(M[0][1], M[0][2],
77 M[1][1], M[1][2]));
78 }
79
boxprod(Point a,Point b,Point c)80 static double boxprod(Point a, Point b, Point c) {
81 return det(a,b) - det(a,c) + det(b,c);
82 }
83
84 class BadConversion : public std::runtime_error {
85 public:
BadConversion(const std::string & s)86 BadConversion(const std::string& s)
87 : std::runtime_error(s)
88 { }
89 };
90
91 template <typename T>
stringify(T x)92 inline std::string stringify(T x)
93 {
94 std::ostringstream o;
95 if (!(o << x))
96 throw BadConversion("stringify(T)");
97 return o.str();
98 }
99
100 /* A G4 continuous cubic parametric approximation for rational quadratics.
101 See
102 An analysis of cubic approximation schemes for conic sections
103 Michael Floater
104 SINTEF
105
106 This is less accurate overall than some of his other schemes, but
107 produces very smooth joins and is still optimally h^-6
108 convergent.
109 */
110
lambda() const111 double RatQuad::lambda() const {
112 return 2*(6*w*w +1 -std::sqrt(3*w*w+1))/(12*w*w+3);
113 }
114
fromPointsTangents(Point P0,Point dP0,Point P,Point P2,Point dP2)115 RatQuad RatQuad::fromPointsTangents(Point P0, Point dP0,
116 Point P,
117 Point P2, Point dP2) {
118 Line Line0 = Line::from_origin_and_vector(P0, dP0);
119 Line Line2 = Line::from_origin_and_vector(P2, dP2);
120 try {
121 OptCrossing oc = intersection(Line0, Line2);
122 if(!oc) // what to do?
123 return RatQuad(Point(), Point(), Point(), 0); // need opt really
124 //assert(0);
125 Point P1 = Line0.pointAt((*oc).ta);
126 double triarea = boxprod(P0, P1, P2);
127 // std::cout << "RatQuad::fromPointsTangents: triarea = " << triarea << std::endl;
128 if (triarea == 0)
129 {
130 return RatQuad(P0, 0.5*(P0+P2), P2, 1);
131 }
132 double tau0 = boxprod(P, P1, P2)/triarea;
133 double tau1 = boxprod(P0, P, P2)/triarea;
134 double tau2 = boxprod(P0, P1, P)/triarea;
135 if (tau0 == 0 || tau1 == 0 || tau2 == 0)
136 {
137 return RatQuad(P0, 0.5*(P0+P2), P2, 1);
138 }
139 double w = tau1/(2*std::sqrt(tau0*tau2));
140 // std::cout << "RatQuad::fromPointsTangents: tau0 = " << tau0 << std::endl;
141 // std::cout << "RatQuad::fromPointsTangents: tau1 = " << tau1 << std::endl;
142 // std::cout << "RatQuad::fromPointsTangents: tau2 = " << tau2 << std::endl;
143 // std::cout << "RatQuad::fromPointsTangents: w = " << w << std::endl;
144 return RatQuad(P0, P1, P2, w);
145 } catch(Geom::InfiniteSolutions const&) {
146 return RatQuad(P0, 0.5*(P0+P2), P2, 1);
147 }
148 return RatQuad(Point(), Point(), Point(), 0); // need opt really
149 }
150
circularArc(Point P0,Point P1,Point P2)151 RatQuad RatQuad::circularArc(Point P0, Point P1, Point P2) {
152 return RatQuad(P0, P1, P2, dot(unit_vector(P0 - P1), unit_vector(P0 - P2)));
153 }
154
155
toCubic() const156 CubicBezier RatQuad::toCubic() const {
157 return toCubic(lambda());
158 }
159
toCubic(double lamb) const160 CubicBezier RatQuad::toCubic(double lamb) const {
161 return CubicBezier(P[0],
162 (1-lamb)*P[0] + lamb*P[1],
163 (1-lamb)*P[2] + lamb*P[1],
164 P[2]);
165 }
166
pointAt(double t) const167 Point RatQuad::pointAt(double t) const {
168 Bezier xt(P[0][0], P[1][0]*w, P[2][0]);
169 Bezier yt(P[0][1], P[1][1]*w, P[2][1]);
170 double wt = Bezier(1, w, 1).valueAt(t);
171 return Point(xt.valueAt(t)/wt,
172 yt.valueAt(t)/wt);
173 }
174
split(RatQuad & a,RatQuad & b) const175 void RatQuad::split(RatQuad &a, RatQuad &b) const {
176 a.P[0] = P[0];
177 b.P[2] = P[2];
178 a.P[1] = (P[0]+w*P[1])/(1+w);
179 b.P[1] = (w*P[1]+P[2])/(1+w);
180 a.w = b.w = std::sqrt((1+w)/2);
181 a.P[2] = b.P[0] = (0.5*a.P[1]+0.5*b.P[1]);
182 }
183
184
hermite() const185 D2<SBasis> RatQuad::hermite() const {
186 SBasis t = Linear(0, 1);
187 SBasis omt = Linear(1, 0);
188
189 D2<SBasis> out(omt*omt*P[0][0]+2*omt*t*P[1][0]*w+t*t*P[2][0],
190 omt*omt*P[0][1]+2*omt*t*P[1][1]*w+t*t*P[2][1]);
191 for(int dim = 0; dim < 2; dim++) {
192 out[dim] = divide(out[dim], (omt*omt+2*omt*t*w+t*t), 2);
193 }
194 return out;
195 }
196
homogeneous() const197 std::vector<SBasis> RatQuad::homogeneous() const {
198 std::vector<SBasis> res(3, SBasis());
199 Bezier xt(P[0][0], P[1][0]*w, P[2][0]);
200 bezier_to_sbasis(res[0],xt);
201 Bezier yt(P[0][1], P[1][1]*w, P[2][1]);
202 bezier_to_sbasis(res[1],yt);
203 Bezier wt(1, w, 1);
204 bezier_to_sbasis(res[2],wt);
205 return res;
206 }
207
208 #if 0
209 std::string xAx::categorise() const {
210 double M[3][3] = {{c[0], c[1], c[3]},
211 {c[1], c[2], c[4]},
212 {c[3], c[4], c[5]}};
213 double D = det3(M);
214 if (c[0] == 0 && c[1] == 0 && c[2] == 0)
215 return "line";
216 std::string res = stringify(D);
217 double descr = c[1]*c[1] - c[0]*c[2];
218 if (descr < 0) {
219 if (c[0] == c[2] && c[1] == 0)
220 return res + "circle";
221 return res + "ellipse";
222 } else if (descr == 0) {
223 return res + "parabola";
224 } else if (descr > 0) {
225 if (c[0] + c[2] == 0) {
226 if (D == 0)
227 return res + "two lines";
228 return res + "rectangular hyperbola";
229 }
230 return res + "hyperbola";
231
232 }
233 return "no idea!";
234 }
235 #endif
236
237
decompose_degenerate(xAx const & C1,xAx const & C2,xAx const & xC0)238 std::vector<Point> decompose_degenerate(xAx const & C1, xAx const & C2, xAx const & xC0) {
239 std::vector<Point> res;
240 double A[2][2] = {{2*xC0.c[0], xC0.c[1]},
241 {xC0.c[1], 2*xC0.c[2]}};
242 //Point B0 = xC0.bottom();
243 double const determ = det(A);
244 //std::cout << determ << "\n";
245 if (fabs(determ) >= 1e-20) { // hopeful, I know
246 Geom::Coord const ideterm = 1.0 / determ;
247
248 double b[2] = {-xC0.c[3], -xC0.c[4]};
249 Point B0((A[1][1]*b[0] -A[0][1]*b[1]),
250 (-A[1][0]*b[0] + A[0][0]*b[1]));
251 B0 *= ideterm;
252 Point n0, n1;
253 // Are these just the eigenvectors of A11?
254 if(xC0.c[0] == xC0.c[2]) {
255 double b = 0.5*xC0.c[1]/xC0.c[0];
256 double c = xC0.c[2]/xC0.c[0];
257 //assert(fabs(b*b-c) > 1e-10);
258 double d = std::sqrt(b*b-c);
259 //assert(fabs(b-d) > 1e-10);
260 n0 = Point(1, b+d);
261 n1 = Point(1, b-d);
262 } else if(fabs(xC0.c[0]) > fabs(xC0.c[2])) {
263 double b = 0.5*xC0.c[1]/xC0.c[0];
264 double c = xC0.c[2]/xC0.c[0];
265 //assert(fabs(b*b-c) > 1e-10);
266 double d = std::sqrt(b*b-c);
267 //assert(fabs(b-d) > 1e-10);
268 n0 = Point(1, b+d);
269 n1 = Point(1, b-d);
270 } else {
271 double b = 0.5*xC0.c[1]/xC0.c[2];
272 double c = xC0.c[0]/xC0.c[2];
273 //assert(fabs(b*b-c) > 1e-10);
274 double d = std::sqrt(b*b-c);
275 //assert(fabs(b-d) > 1e-10);
276 n0 = Point(b+d, 1);
277 n1 = Point(b-d, 1);
278 }
279
280 Line L0 = Line::from_origin_and_vector(B0, rot90(n0));
281 Line L1 = Line::from_origin_and_vector(B0, rot90(n1));
282
283 std::vector<double> rts = C1.roots(L0);
284 for(double rt : rts) {
285 Point P = L0.pointAt(rt);
286 res.push_back(P);
287 }
288 rts = C1.roots(L1);
289 for(double rt : rts) {
290 Point P = L1.pointAt(rt);
291 res.push_back(P);
292 }
293 } else {
294 // single or double line
295 // check for completely zero case (what to do?)
296 assert(xC0.c[0] || xC0.c[1] ||
297 xC0.c[2] || xC0.c[3] ||
298 xC0.c[4] || xC0.c[5]);
299 Point trial_pt(0,0);
300 Point g = xC0.gradient(trial_pt);
301 if(L2sq(g) == 0) {
302 trial_pt[0] += 1;
303 g = xC0.gradient(trial_pt);
304 if(L2sq(g) == 0) {
305 trial_pt[1] += 1;
306 g = xC0.gradient(trial_pt);
307 if(L2sq(g) == 0) {
308 trial_pt[0] += 1;
309 g = xC0.gradient(trial_pt);
310 if(L2sq(g) == 0) {
311 trial_pt = Point(1.5,0.5);
312 g = xC0.gradient(trial_pt);
313 }
314 }
315 }
316 }
317 //std::cout << trial_pt << ", " << g << "\n";
318 /**
319 * At this point we have tried up to 4 points: 0,0, 1,0, 1,1, 2,1, 1.5,1.5
320 *
321 * No degenerate conic can pass through these points, so we can assume
322 * that we've found a perpendicular to the double line.
323 * Proof:
324 * any degenerate must consist of at most 2 lines. 1.5,0.5 is not on any pair of lines
325 * passing through the previous 4 trials.
326 *
327 * alternatively, there may be a way to determine this directly from xC0
328 */
329 assert(L2sq(g) != 0);
330
331 Line Lx = Line::from_origin_and_vector(trial_pt, g); // a line along the gradient
332 std::vector<double> rts = xC0.roots(Lx);
333 for(double rt : rts) {
334 Point P0 = Lx.pointAt(rt);
335 //std::cout << P0 << "\n";
336 Line L = Line::from_origin_and_vector(P0, rot90(g));
337 std::vector<double> cnrts;
338 // It's very likely that at least one of the conics is degenerate, this will hopefully pick the more generate of the two.
339 if(fabs(C1.hessian().det()) > fabs(C2.hessian().det()))
340 cnrts = C1.roots(L);
341 else
342 cnrts = C2.roots(L);
343 for(double cnrt : cnrts) {
344 Point P = L.pointAt(cnrt);
345 res.push_back(P);
346 }
347 }
348 }
349 return res;
350 }
351
xAx_descr(xAx const & C)352 double xAx_descr(xAx const & C) {
353 double mC[3][3] = {{C.c[0], (C.c[1])/2, (C.c[3])/2},
354 {(C.c[1])/2, C.c[2], (C.c[4])/2},
355 {(C.c[3])/2, (C.c[4])/2, C.c[5]}};
356
357 return det3(mC);
358 }
359
360
intersect(xAx const & C1,xAx const & C2)361 std::vector<Point> intersect(xAx const & C1, xAx const & C2) {
362 // You know, if either of the inputs are degenerate we should use them first!
363 if(xAx_descr(C1) == 0) {
364 return decompose_degenerate(C1, C2, C1);
365 }
366 if(xAx_descr(C2) == 0) {
367 return decompose_degenerate(C1, C2, C2);
368 }
369 std::vector<Point> res;
370 SBasis T(Linear(-1,1));
371 SBasis S(Linear(1,1));
372 SBasis C[3][3] = {{T*C1.c[0]+S*C2.c[0], (T*C1.c[1]+S*C2.c[1])/2, (T*C1.c[3]+S*C2.c[3])/2},
373 {(T*C1.c[1]+S*C2.c[1])/2, T*C1.c[2]+S*C2.c[2], (T*C1.c[4]+S*C2.c[4])/2},
374 {(T*C1.c[3]+S*C2.c[3])/2, (T*C1.c[4]+S*C2.c[4])/2, T*C1.c[5]+S*C2.c[5]}};
375
376 SBasis D = det3(C);
377 std::vector<double> rts = Geom::roots(D);
378 if(rts.empty()) {
379 T = Linear(1,1);
380 S = Linear(-1,1);
381 SBasis C[3][3] = {{T*C1.c[0]+S*C2.c[0], (T*C1.c[1]+S*C2.c[1])/2, (T*C1.c[3]+S*C2.c[3])/2},
382 {(T*C1.c[1]+S*C2.c[1])/2, T*C1.c[2]+S*C2.c[2], (T*C1.c[4]+S*C2.c[4])/2},
383 {(T*C1.c[3]+S*C2.c[3])/2, (T*C1.c[4]+S*C2.c[4])/2, T*C1.c[5]+S*C2.c[5]}};
384
385 D = det3(C);
386 rts = Geom::roots(D);
387 }
388 // at this point we have a T and S and perhaps some roots that represent our degenerate conic
389 // Let's just pick one randomly (can we do better?)
390 //for(unsigned i = 0; i < rts.size(); i++) {
391 if(!rts.empty()) {
392 unsigned i = 0;
393 double t = T.valueAt(rts[i]);
394 double s = S.valueAt(rts[i]);
395 xAx xC0 = C1*t + C2*s;
396 //::draw(cr, xC0, screen_rect); // degen
397
398 return decompose_degenerate(C1, C2, xC0);
399
400
401 } else {
402 std::cout << "What?" << std::endl;
403 ;//std::cout << D << "\n";
404 }
405 return res;
406 }
407
408
fromPoint(Point p)409 xAx xAx::fromPoint(Point p) {
410 return xAx(1., 0, 1., -2*p[0], -2*p[1], dot(p,p));
411 }
412
fromDistPoint(Point,double)413 xAx xAx::fromDistPoint(Point /*p*/, double /*d*/) {
414 return xAx();//1., 0, 1., -2*(1+d)*p[0], -2*(1+d)*p[1], dot(p,p)+d*d);
415 }
416
fromLine(Point n,double d)417 xAx xAx::fromLine(Point n, double d) {
418 return xAx(n[0]*n[0], 2*n[0]*n[1], n[1]*n[1], 2*d*n[0], 2*d*n[1], d*d);
419 }
420
fromLine(Line l)421 xAx xAx::fromLine(Line l) {
422 double dist;
423 Point norm = l.normalAndDist(dist);
424
425 return fromLine(norm, dist);
426 }
427
fromPoints(std::vector<Geom::Point> const & pt)428 xAx xAx::fromPoints(std::vector<Geom::Point> const &pt) {
429 Geom::NL::Vector V(pt.size(), -1.0);
430 Geom::NL::Matrix M(pt.size(), 5);
431 for(unsigned i = 0; i < pt.size(); i++) {
432 Geom::Point P = pt[i];
433 Geom::NL::VectorView vv = M.row_view(i);
434 vv[0] = P[0]*P[0];
435 vv[1] = P[0]*P[1];
436 vv[2] = P[1]*P[1];
437 vv[3] = P[0];
438 vv[4] = P[1];
439 }
440
441 Geom::NL::LinearSystem ls(M, V);
442
443 Geom::NL::Vector x = ls.SV_solve();
444 return Geom::xAx(x[0], x[1], x[2], x[3], x[4], 1);
445
446 }
447
448
449
valueAt(Point P) const450 double xAx::valueAt(Point P) const {
451 return evaluate_at(P[0], P[1]);
452 }
453
scale(double sx,double sy) const454 xAx xAx::scale(double sx, double sy) const {
455 return xAx(c[0]*sx*sx, c[1]*sx*sy, c[2]*sy*sy,
456 c[3]*sx, c[4]*sy, c[5]);
457 }
458
gradient(Point p) const459 Point xAx::gradient(Point p) const{
460 double x = p[0];
461 double y = p[1];
462 return Point(2*c[0]*x + c[1]*y + c[3],
463 c[1]*x + 2*c[2]*y + c[4]);
464 }
465
operator -(xAx const & b) const466 xAx xAx::operator-(xAx const &b) const {
467 xAx res;
468 for(int i = 0; i < 6; i++) {
469 res.c[i] = c[i] - b.c[i];
470 }
471 return res;
472 }
operator +(xAx const & b) const473 xAx xAx::operator+(xAx const &b) const {
474 xAx res;
475 for(int i = 0; i < 6; i++) {
476 res.c[i] = c[i] + b.c[i];
477 }
478 return res;
479 }
operator +(double const & b) const480 xAx xAx::operator+(double const &b) const {
481 xAx res;
482 for(int i = 0; i < 5; i++) {
483 res.c[i] = c[i];
484 }
485 res.c[5] = c[5] + b;
486 return res;
487 }
488
operator *(double const & b) const489 xAx xAx::operator*(double const &b) const {
490 xAx res;
491 for(int i = 0; i < 6; i++) {
492 res.c[i] = c[i] * b;
493 }
494 return res;
495 }
496
crossings(Rect r) const497 std::vector<Point> xAx::crossings(Rect r) const {
498 std::vector<Point> res;
499 for(int ei = 0; ei < 4; ei++) {
500 Geom::LineSegment ls(r.corner(ei), r.corner(ei+1));
501 D2<SBasis> lssb = ls.toSBasis();
502 SBasis edge_curve = evaluate_at(lssb[0], lssb[1]);
503 std::vector<double> rts = Geom::roots(edge_curve);
504 for(double rt : rts) {
505 res.push_back(lssb.valueAt(rt));
506 }
507 }
508 return res;
509 }
510
toCurve(Rect const & bnd) const511 std::optional<RatQuad> xAx::toCurve(Rect const & bnd) const {
512 std::vector<Point> crs = crossings(bnd);
513 if(crs.size() == 1) {
514 Point A = crs[0];
515 Point dA = rot90(gradient(A));
516 if(L2sq(dA) <= 1e-10) { // perhaps a single point?
517 return std::optional<RatQuad> ();
518 }
519 LineSegment ls = intersection(Line::from_origin_and_vector(A, dA), bnd);
520 return RatQuad::fromPointsTangents(A, dA, ls.pointAt(0.5), ls[1], dA);
521 }
522 else if(crs.size() >= 2 && crs.size() < 4) {
523 Point A = crs[0];
524 Point C = crs[1];
525 if(crs.size() == 3) {
526 if(distance(A, crs[2]) > distance(A, C))
527 C = crs[2];
528 else if(distance(C, crs[2]) > distance(A, C))
529 A = crs[2];
530 }
531 Line bisector = make_bisector_line(LineSegment(A, C));
532 std::vector<double> bisect_rts = this->roots(bisector);
533 if(!bisect_rts.empty()) {
534 int besti = -1;
535 for(unsigned i =0; i < bisect_rts.size(); i++) {
536 Point p = bisector.pointAt(bisect_rts[i]);
537 if(bnd.contains(p)) {
538 besti = i;
539 }
540 }
541 if(besti >= 0) {
542 Point B = bisector.pointAt(bisect_rts[besti]);
543
544 Point dA = gradient(A);
545 Point dC = gradient(C);
546 if(L2sq(dA) <= 1e-10 || L2sq(dC) <= 1e-10) {
547 return RatQuad::fromPointsTangents(A, C-A, B, C, A-C);
548 }
549
550 RatQuad rq = RatQuad::fromPointsTangents(A, rot90(dA),
551 B, C, rot90(dC));
552 return rq;
553 //std::vector<SBasis> hrq = rq.homogeneous();
554 /*SBasis vertex_poly = evaluate_at(hrq[0], hrq[1], hrq[2]);
555 std::vector<double> rts = roots(vertex_poly);
556 for(unsigned i = 0; i < rts.size(); i++) {
557 //draw_circ(cr, Point(rq.pointAt(rts[i])));
558 }*/
559 }
560 }
561 }
562 return std::optional<RatQuad>();
563 }
564
roots(Point d,Point o) const565 std::vector<double> xAx::roots(Point d, Point o) const {
566 // Find the roots on line l
567 // form the quadratic Q(t) = 0 by composing l with xAx
568 double q2 = c[0]*d[0]*d[0] + c[1]*d[0]*d[1] + c[2]*d[1]*d[1];
569 double q1 = (2*c[0]*d[0]*o[0] +
570 c[1]*(d[0]*o[1]+d[1]*o[0]) +
571 2*c[2]*d[1]*o[1] +
572 c[3]*d[0] + c[4]*d[1]);
573 double q0 = c[0]*o[0]*o[0] + c[1]*o[0]*o[1] + c[2]*o[1]*o[1] + c[3]*o[0] + c[4]*o[1] + c[5];
574 std::vector<double> r;
575 if(q2 == 0) {
576 if(q1 == 0) {
577 return r;
578 }
579 r.push_back(-q0/q1);
580 } else {
581 double desc = q1*q1 - 4*q2*q0;
582 /*std::cout << q2 << ", "
583 << q1 << ", "
584 << q0 << "; "
585 << desc << "\n";*/
586 if (desc < 0)
587 return r;
588 else if (desc == 0)
589 r.push_back(-q1/(2*q2));
590 else {
591 desc = std::sqrt(desc);
592 double t;
593 if (q1 == 0)
594 {
595 t = -0.5 * desc;
596 }
597 else
598 {
599 t = -0.5 * (q1 + sgn(q1) * desc);
600 }
601 r.push_back(t/q2);
602 r.push_back(q0/t);
603 }
604 }
605 return r;
606 }
607
roots(Line const & l) const608 std::vector<double> xAx::roots(Line const &l) const {
609 return roots(l.versor(), l.origin());
610 }
611
quad_ex(double a,double b,double c,Interval ivl)612 Interval xAx::quad_ex(double a, double b, double c, Interval ivl) {
613 double cx = -b*0.5/a;
614 Interval bnds((a*ivl.min()+b)*ivl.min()+c, (a*ivl.max()+b)*ivl.max()+c);
615 if(ivl.contains(cx))
616 bnds.expandTo((a*cx+b)*cx+c);
617 return bnds;
618 }
619
hessian() const620 Geom::Affine xAx::hessian() const {
621 Geom::Affine m(2*c[0], c[1],
622 c[1], 2*c[2],
623 0, 0);
624 return m;
625 }
626
627
solve(double A[2][2],double b[2])628 std::optional<Point> solve(double A[2][2], double b[2]) {
629 double const determ = det(A);
630 if (determ != 0.0) { // hopeful, I know
631 Geom::Coord const ideterm = 1.0 / determ;
632
633 return Point ((A[1][1]*b[0] -A[0][1]*b[1]),
634 (-A[1][0]*b[0] + A[0][0]*b[1]))* ideterm;
635 } else {
636 return std::optional<Point>();
637 }
638 }
639
bottom() const640 std::optional<Point> xAx::bottom() const {
641 double A[2][2] = {{2*c[0], c[1]},
642 {c[1], 2*c[2]}};
643 double b[2] = {-c[3], -c[4]};
644 return solve(A, b);
645 //return Point(-c[3], -c[4])*hessian().inverse();
646 }
647
extrema(Rect r) const648 Interval xAx::extrema(Rect r) const {
649 if (c[0] == 0 && c[1] == 0 && c[2] == 0) {
650 Interval ext(valueAt(r.corner(0)));
651 for(int i = 1; i < 4; i++)
652 ext |= Interval(valueAt(r.corner(i)));
653 return ext;
654 }
655 double k = r[X].min();
656 Interval ext = quad_ex(c[2], c[1]*k+c[4], (c[0]*k + c[3])*k + c[5], r[Y]);
657 k = r[X].max();
658 ext |= quad_ex(c[2], c[1]*k+c[4], (c[0]*k + c[3])*k + c[5], r[Y]);
659 k = r[Y].min();
660 ext |= quad_ex(c[0], c[1]*k+c[3], (c[2]*k + c[4])*k + c[5], r[X]);
661 k = r[Y].max();
662 ext |= quad_ex(c[0], c[1]*k+c[3], (c[2]*k + c[4])*k + c[5], r[X]);
663 std::optional<Point> B0 = bottom();
664 if (B0 && r.contains(*B0))
665 ext.expandTo(0);
666 return ext;
667 }
668
669
670
671
672
673
674
675
676
677 /*
678 * helper functions
679 */
680
at_infinity(Point const & p)681 bool at_infinity (Point const& p)
682 {
683 if (p[X] == infinity() || p[X] == -infinity()
684 || p[Y] == infinity() || p[Y] == -infinity())
685 {
686 return true;
687 }
688 return false;
689 }
690
691 inline
signed_triangle_area(Point const & p1,Point const & p2,Point const & p3)692 double signed_triangle_area (Point const& p1, Point const& p2, Point const& p3)
693 {
694 return (cross(p2, p3) - cross(p1, p3) + cross(p1, p2));
695 }
696
697
698
699 /*
700 * Define a conic section by computing the one that fits better with
701 * N points.
702 *
703 * points: points to fit
704 *
705 * precondition: there must be at least 5 non-overlapping points
706 */
set(std::vector<Point> const & points)707 void xAx::set(std::vector<Point> const& points)
708 {
709 size_t sz = points.size();
710 if (sz < 5)
711 {
712 THROW_RANGEERROR("fitting error: too few points passed");
713 }
714 NL::LFMConicSection model;
715 NL::least_squeares_fitter<NL::LFMConicSection> fitter(model, sz);
716
717 for (size_t i = 0; i < sz; ++i)
718 {
719 fitter.append(points[i]);
720 }
721 fitter.update();
722
723 NL::Vector z(sz, 0.0);
724 model.instance(*this, fitter.result(z));
725 }
726
727 /*
728 * Define a section conic by providing the coordinates of one of its vertex,
729 * the major axis inclination angle and the coordinates of its foci
730 * with respect to the unidimensional system defined by the major axis with
731 * origin set at the provided vertex.
732 *
733 * _vertex : section conic vertex V
734 * _angle : section conic major axis angle
735 * _dist1: +/-distance btw V and nearest focus
736 * _dist2: +/-distance btw V and farest focus
737 *
738 * prerequisite: _dist1 <= _dist2
739 */
set(const Point & _vertex,double _angle,double _dist1,double _dist2)740 void xAx::set (const Point& _vertex, double _angle, double _dist1, double _dist2)
741 {
742 using std::swap;
743
744 if (_dist2 == infinity() || _dist2 == -infinity()) // parabola
745 {
746 if (_dist1 == infinity()) // degenerate to a line
747 {
748 Line l(_vertex, _angle);
749 std::vector<double> lcoeff = l.coefficients();
750 coeff(3) = lcoeff[0];
751 coeff(4) = lcoeff[1];
752 coeff(5) = lcoeff[2];
753 return;
754 }
755
756 // y^2 - 4px == 0
757 double cD = -4 * _dist1;
758
759 double cosa = std::cos (_angle);
760 double sina = std::sin (_angle);
761 double cca = cosa * cosa;
762 double ssa = sina * sina;
763 double csa = cosa * sina;
764
765 coeff(0) = ssa;
766 coeff(1) = -2 * csa;
767 coeff(2) = cca;
768 coeff(3) = cD * cosa;
769 coeff(4) = cD * sina;
770
771 double VxVx = _vertex[X] * _vertex[X];
772 double VxVy = _vertex[X] * _vertex[Y];
773 double VyVy = _vertex[Y] * _vertex[Y];
774
775 coeff(5) = coeff(0) * VxVx + coeff(1) * VxVy + coeff(2) * VyVy
776 - coeff(3) * _vertex[X] - coeff(4) * _vertex[Y];
777 coeff(3) -= (2 * coeff(0) * _vertex[X] + coeff(1) * _vertex[Y]);
778 coeff(4) -= (2 * coeff(2) * _vertex[Y] + coeff(1) * _vertex[X]);
779
780 return;
781 }
782
783 if (std::fabs(_dist1) > std::fabs(_dist2))
784 {
785 swap (_dist1, _dist2);
786 }
787 if (_dist1 < 0)
788 {
789 _angle -= M_PI;
790 _dist1 = -_dist1;
791 _dist2 = -_dist2;
792 }
793
794 // ellipse and hyperbola
795 double lin_ecc = (_dist2 - _dist1) / 2;
796 double rx = (_dist2 + _dist1) / 2;
797
798 double cA = rx * rx - lin_ecc * lin_ecc;
799 double cC = rx * rx;
800 double cF = - cA * cC;
801 // std::cout << "cA: " << cA << std::endl;
802 // std::cout << "cC: " << cC << std::endl;
803 // std::cout << "cF: " << cF << std::endl;
804
805 double cosa = std::cos (_angle);
806 double sina = std::sin (_angle);
807 double cca = cosa * cosa;
808 double ssa = sina * sina;
809 double csa = cosa * sina;
810
811 coeff(0) = cca * cA + ssa * cC;
812 coeff(2) = ssa * cA + cca * cC;
813 coeff(1) = 2 * csa * (cA - cC);
814
815 Point C (rx * cosa + _vertex[X], rx * sina + _vertex[Y]);
816 double CxCx = C[X] * C[X];
817 double CxCy = C[X] * C[Y];
818 double CyCy = C[Y] * C[Y];
819
820 coeff(3) = -2 * coeff(0) * C[X] - coeff(1) * C[Y];
821 coeff(4) = -2 * coeff(2) * C[Y] - coeff(1) * C[X];
822 coeff(5) = cF + coeff(0) * CxCx + coeff(1) * CxCy + coeff(2) * CyCy;
823 }
824
825 /*
826 * Define a conic section by providing one of its vertex and its foci.
827 *
828 * _vertex: section conic vertex
829 * _focus1: section conic focus
830 * _focus2: section conic focus
831 */
set(const Point & _vertex,const Point & _focus1,const Point & _focus2)832 void xAx::set (const Point& _vertex, const Point& _focus1, const Point& _focus2)
833 {
834 if (at_infinity(_vertex))
835 {
836 THROW_RANGEERROR("case not handled: vertex at infinity");
837 }
838 if (at_infinity(_focus2))
839 {
840 if (at_infinity(_focus1))
841 {
842 THROW_RANGEERROR("case not handled: both focus at infinity");
843 }
844 Point VF = _focus1 - _vertex;
845 double dist1 = L2(VF);
846 double angle = atan2(VF);
847 set(_vertex, angle, dist1, infinity());
848 return;
849 }
850 else if (at_infinity(_focus1))
851 {
852 Point VF = _focus2 - _vertex;
853 double dist1 = L2(VF);
854 double angle = atan2(VF);
855 set(_vertex, angle, dist1, infinity());
856 return;
857 }
858 assert (are_collinear (_vertex, _focus1, _focus2));
859 if (!are_near(_vertex, _focus1))
860 {
861 Point VF = _focus1 - _vertex;
862 Line axis(_vertex, _focus1);
863 double angle = atan2(VF);
864 double dist1 = L2(VF);
865 double dist2 = distance (_vertex, _focus2);
866 double t = axis.timeAt(_focus2);
867 if (t < 0) dist2 = -dist2;
868 // std::cout << "t = " << t << std::endl;
869 // std::cout << "dist2 = " << dist2 << std::endl;
870 set (_vertex, angle, dist1, dist2);
871 }
872 else if (!are_near(_vertex, _focus2))
873 {
874 Point VF = _focus2 - _vertex;
875 double angle = atan2(VF);
876 double dist1 = 0;
877 double dist2 = L2(VF);
878 set (_vertex, angle, dist1, dist2);
879 }
880 else
881 {
882 coeff(0) = coeff(2) = 1;
883 coeff(1) = coeff(3) = coeff(4) = coeff(5) = 0;
884 }
885 }
886
887 /*
888 * Define a conic section by passing a focus, the related directrix,
889 * and the eccentricity (e)
890 * (e < 1 -> ellipse; e = 1 -> parabola; e > 1 -> hyperbola)
891 *
892 * _focus: a focus of the conic section
893 * _directrix: the directrix related to the given focus
894 * _eccentricity: the eccentricity parameter of the conic section
895 */
set(const Point & _focus,const Line & _directrix,double _eccentricity)896 void xAx::set (const Point & _focus, const Line & _directrix, double _eccentricity)
897 {
898 Point O = _directrix.pointAt (_directrix.timeAtProjection (_focus));
899 //std::cout << "O = " << O << std::endl;
900 Point OF = _focus - O;
901 double p = L2(OF);
902
903 coeff(0) = 1 - _eccentricity * _eccentricity;
904 coeff(1) = 0;
905 coeff(2) = 1;
906 coeff(3) = -2 * p;
907 coeff(4) = 0;
908 coeff(5) = p * p;
909
910 double angle = atan2 (OF);
911
912 (*this) = rotate (angle);
913 //std::cout << "O = " << O << std::endl;
914 (*this) = translate (O);
915 }
916
917 /*
918 * Made up a degenerate conic section as a pair of lines
919 *
920 * l1, l2: lines that made up the conic section
921 */
set(const Line & l1,const Line & l2)922 void xAx::set (const Line& l1, const Line& l2)
923 {
924 std::vector<double> cl1 = l1.coefficients();
925 std::vector<double> cl2 = l2.coefficients();
926
927 coeff(0) = cl1[0] * cl2[0];
928 coeff(2) = cl1[1] * cl2[1];
929 coeff(5) = cl1[2] * cl2[2];
930 coeff(1) = cl1[0] * cl2[1] + cl1[1] * cl2[0];
931 coeff(3) = cl1[0] * cl2[2] + cl1[2] * cl2[0];
932 coeff(4) = cl1[1] * cl2[2] + cl1[2] * cl2[1];
933 }
934
935
936
937 /*
938 * Return the section conic kind
939 */
kind() const940 xAx::kind_t xAx::kind () const
941 {
942
943 xAx conic(*this);
944 NL::SymmetricMatrix<3> C = conic.get_matrix();
945 NL::ConstSymmetricMatrixView<2> A = C.main_minor_const_view();
946
947 double t1 = trace(A);
948 double t2 = det(A);
949 //double T3 = det(C);
950 int st1 = trace_sgn(A);
951 int st2 = det_sgn(A);
952 int sT3 = det_sgn(C);
953
954 //std::cout << "T3 = " << T3 << std::endl;
955 //std::cout << "sT3 = " << sT3 << std::endl;
956 //std::cout << "t2 = " << t2 << std::endl;
957 //std::cout << "t1 = " << t1 << std::endl;
958 //std::cout << "st2 = " << st2 << std::endl;
959
960 if (sT3 != 0)
961 {
962 if (st2 == 0)
963 {
964 return PARABOLA;
965 }
966 else if (st2 == 1)
967 {
968
969 if (sT3 * st1 < 0)
970 {
971 NL::SymmetricMatrix<2> discr;
972 discr(0,0) = 4; discr(1,1) = t2; discr(1,0) = t1;
973 int discr_sgn = - det_sgn (discr);
974 //std::cout << "t1 * t1 - 4 * t2 = "
975 // << (t1 * t1 - 4 * t2) << std::endl;
976 //std::cout << "discr_sgn = " << discr_sgn << std::endl;
977 if (discr_sgn == 0)
978 {
979 return CIRCLE;
980 }
981 else
982 {
983 return REAL_ELLIPSE;
984 }
985 }
986 else // sT3 * st1 > 0
987 {
988 return IMAGINARY_ELLIPSE;
989 }
990 }
991 else // t2 < 0
992 {
993 if (st1 == 0)
994 {
995 return RECTANGULAR_HYPERBOLA;
996 }
997 else
998 {
999 return HYPERBOLA;
1000 }
1001 }
1002 }
1003 else // T3 == 0
1004 {
1005 if (st2 == 0)
1006 {
1007 //double T2 = NL::trace<2>(C);
1008 int sT2 = NL::trace_sgn<2>(C);
1009 //std::cout << "T2 = " << T2 << std::endl;
1010 //std::cout << "sT2 = " << sT2 << std::endl;
1011
1012 if (sT2 == 0)
1013 {
1014 return DOUBLE_LINE;
1015 }
1016 if (sT2 == -1)
1017 {
1018 return TWO_REAL_PARALLEL_LINES;
1019 }
1020 else // T2 > 0
1021 {
1022 return TWO_IMAGINARY_PARALLEL_LINES;
1023 }
1024 }
1025 else if (st2 == -1)
1026 {
1027 return TWO_REAL_CROSSING_LINES;
1028 }
1029 else // t2 > 0
1030 {
1031 return TWO_IMAGINARY_CROSSING_LINES;
1032 }
1033 }
1034 return UNKNOWN;
1035 }
1036
1037 /*
1038 * Return a string representing the conic section kind
1039 */
categorise() const1040 std::string xAx::categorise() const
1041 {
1042 kind_t KIND = kind();
1043
1044 switch (KIND)
1045 {
1046 case PARABOLA :
1047 return "parabola";
1048 case CIRCLE :
1049 return "circle";
1050 case REAL_ELLIPSE :
1051 return "real ellispe";
1052 case IMAGINARY_ELLIPSE :
1053 return "imaginary ellispe";
1054 case RECTANGULAR_HYPERBOLA :
1055 return "rectangular hyperbola";
1056 case HYPERBOLA :
1057 return "hyperbola";
1058 case DOUBLE_LINE :
1059 return "double line";
1060 case TWO_REAL_PARALLEL_LINES :
1061 return "two real parallel lines";
1062 case TWO_IMAGINARY_PARALLEL_LINES :
1063 return "two imaginary parallel lines";
1064 case TWO_REAL_CROSSING_LINES :
1065 return "two real crossing lines";
1066 case TWO_IMAGINARY_CROSSING_LINES :
1067 return "two imaginary crossing lines";
1068 default :
1069 return "unknown";
1070 }
1071 }
1072
1073 /*
1074 * Compute the solutions of the conic section algebraic equation with respect to
1075 * one coordinate after substituting to the other coordinate the passed value
1076 *
1077 * sol: the computed solutions
1078 * v: the provided value
1079 * d: the index of the coordinate the passed value have to be substituted to
1080 */
roots(std::vector<double> & sol,Coord v,Dim2 d) const1081 void xAx::roots (std::vector<double>& sol, Coord v, Dim2 d) const
1082 {
1083 sol.clear();
1084 if (d < 0 || d > Y)
1085 {
1086 THROW_RANGEERROR("dimension parameter out of range");
1087 }
1088
1089 // p*t^2 + q*t + r = 0;
1090 double p, q, r;
1091
1092 if (d == X)
1093 {
1094 p = coeff(2);
1095 q = coeff(4) + coeff(1) * v;
1096 r = coeff(5) + (coeff(0) * v + coeff(3)) * v;
1097 }
1098 else
1099 {
1100 p = coeff(0);
1101 q = coeff(3) + coeff(1) * v;
1102 r = coeff(5) + (coeff(2) * v + coeff(4)) * v;
1103 }
1104
1105 if (p == 0)
1106 {
1107 if (q == 0) return;
1108 double t = -r/q;
1109 sol.push_back(t);
1110 return;
1111 }
1112
1113 if (q == 0)
1114 {
1115 if ((p > 0 && r > 0) || (p < 0 && r < 0)) return;
1116 double t = -r / p;
1117 t = std::sqrt (t);
1118 sol.push_back(-t);
1119 sol.push_back(t);
1120 return;
1121 }
1122
1123 if (r == 0)
1124 {
1125 double t = -q/p;
1126 sol.push_back(0);
1127 sol.push_back(t);
1128 return;
1129 }
1130
1131
1132 //std::cout << "p = " << p << ", q = " << q << ", r = " << r << std::endl;
1133 double delta = q * q - 4 * p * r;
1134 if (delta < 0) return;
1135 if (delta == 0)
1136 {
1137 double t = -q / (2 * p);
1138 sol.push_back(t);
1139 return;
1140 }
1141 // else
1142 double srd = std::sqrt(delta);
1143 double t = - (q + sgn(q) * srd) / 2;
1144 sol.push_back (t/p);
1145 sol.push_back (r/t);
1146
1147 }
1148
1149 /*
1150 * Return the inclination angle of the major axis of the conic section
1151 */
axis_angle() const1152 double xAx::axis_angle() const
1153 {
1154 if (coeff(0) == 0 && coeff(1) == 0 && coeff(2) == 0)
1155 {
1156 Line l (coeff(3), coeff(4), coeff(5));
1157 return l.angle();
1158 }
1159 if (coeff(1) == 0 && (coeff(0) == coeff(2))) return 0;
1160
1161 double angle;
1162
1163 int sgn_discr = det_sgn (get_matrix().main_minor_const_view());
1164 if (sgn_discr == 0)
1165 {
1166 //std::cout << "rotation_angle: sgn_discr = "
1167 // << sgn_discr << std::endl;
1168 angle = std::atan2 (-coeff(1), 2 * coeff(2));
1169 if (angle < 0) angle += 2*M_PI;
1170 if (angle >= M_PI) angle -= M_PI;
1171
1172 }
1173 else
1174 {
1175 angle = std::atan2 (coeff(1), coeff(0) - coeff(2));
1176 if (angle < 0) angle += 2*M_PI;
1177 angle -= M_PI;
1178 if (angle < 0) angle += 2*M_PI;
1179 angle /= 2;
1180 if (angle >= M_PI) angle -= M_PI;
1181 }
1182 //std::cout << "rotation_angle : angle = " << angle << std::endl;
1183 return angle;
1184 }
1185
1186 /*
1187 * Translate the conic section by the given vector offset
1188 *
1189 * _offset: represent the vector offset
1190 */
translate(const Point & _offset) const1191 xAx xAx::translate (const Point & _offset) const
1192 {
1193 double B = coeff(1) / 2;
1194 double D = coeff(3) / 2;
1195 double E = coeff(4) / 2;
1196
1197 Point T = - _offset;
1198
1199 xAx cs;
1200 cs.coeff(0) = coeff(0);
1201 cs.coeff(1) = coeff(1);
1202 cs.coeff(2) = coeff(2);
1203
1204 Point DE;
1205 DE[0] = coeff(0) * T[0] + B * T[1];
1206 DE[1] = B * T[0] + coeff(2) * T[1];
1207
1208 cs.coeff(3) = (DE[0] + D) * 2;
1209 cs.coeff(4) = (DE[1] + E) * 2;
1210
1211 cs.coeff(5) = dot (T, DE) + 2 * (T[0] * D + T[1] * E) + coeff(5);
1212
1213 return cs;
1214 }
1215
1216
1217 /*
1218 * Rotate the conic section by the given angle wrt the point (0,0)
1219 *
1220 * angle: represent the rotation angle
1221 */
rotate(double angle) const1222 xAx xAx::rotate (double angle) const
1223 {
1224 double c = std::cos(-angle);
1225 double s = std::sin(-angle);
1226 double cc = c * c;
1227 double ss = s * s;
1228 double cs = c * s;
1229
1230 xAx result;
1231 result.coeff(5) = coeff(5);
1232
1233 // quadratic terms
1234 double Bcs = coeff(1) * cs;
1235
1236 result.coeff(0) = coeff(0) * cc + Bcs + coeff(2) * ss;
1237 result.coeff(2) = coeff(0) * ss - Bcs + coeff(2) * cc;
1238 result.coeff(1) = coeff(1) * (cc - ss) + 2 * (coeff(2) - coeff(0)) * cs;
1239
1240 // linear terms
1241 result.coeff(3) = coeff(3) * c + coeff(4) * s;
1242 result.coeff(4) = coeff(4) * c - coeff(3) * s;
1243
1244 return result;
1245 }
1246
1247
1248 /*
1249 * Decompose a degenerate conic in two lines the conic section is made by.
1250 * Return true if the decomposition is successful, else if it fails.
1251 *
1252 * l1, l2: out parameters where the decomposed conic section is returned
1253 */
decompose(Line & l1,Line & l2) const1254 bool xAx::decompose (Line& l1, Line& l2) const
1255 {
1256 NL::SymmetricMatrix<3> C = get_matrix();
1257 if (!is_quadratic() || !isDegenerate())
1258 {
1259 return false;
1260 }
1261 NL::Matrix M(C);
1262 NL::SymmetricMatrix<3> D = -adj(C);
1263
1264 if (!D.is_zero()) // D == 0 <=> rank(C) < 2
1265 {
1266
1267 //if (D.get<0,0>() < 0 || D.get<1,1>() < 0 || D.get<2,2>() < 0)
1268 //{
1269 //std::cout << "C: \n" << C << std::endl;
1270 //std::cout << "D: \n" << D << std::endl;
1271
1272 /*
1273 * This case should be impossible because any diagonal element
1274 * of D is a square, but due to non exact aritmethic computation
1275 * it can actually happen; however the algorithm seems to work
1276 * correctly even if some diagonal term is negative, the only
1277 * difference is that we should compute the absolute value of
1278 * diagonal elements. So until we elaborate a better degenerate
1279 * test it's better not rising exception when we have a negative
1280 * diagonal element.
1281 */
1282 //}
1283
1284 NL::Vector d(3);
1285 d[0] = std::fabs (D.get<0,0>());
1286 d[1] = std::fabs (D.get<1,1>());
1287 d[2] = std::fabs (D.get<2,2>());
1288
1289 size_t idx = d.max_index();
1290 if (d[idx] == 0)
1291 {
1292 THROW_LOGICALERROR ("xAx::decompose: "
1293 "rank 2 but adjoint with null diagonal");
1294 }
1295 d[0] = D(idx,0); d[1] = D(idx,1); d[2] = D(idx,2);
1296 d.scale (1 / std::sqrt (std::fabs (D(idx,idx))));
1297 M(1,2) += d[0]; M(2,1) -= d[0];
1298 M(0,2) -= d[1]; M(2,0) += d[1];
1299 M(0,1) += d[2]; M(1,0) -= d[2];
1300
1301 //std::cout << "C: \n" << C << std::endl;
1302 //std::cout << "D: \n" << D << std::endl;
1303 //std::cout << "d = " << d << std::endl;
1304 //std::cout << "M = " << M << std::endl;
1305 }
1306
1307 std::pair<size_t, size_t> max_ij = M.max_index();
1308 std::pair<size_t, size_t> min_ij = M.min_index();
1309 double abs_max = std::fabs (M(max_ij.first, max_ij.second));
1310 double abs_min = std::fabs (M(min_ij.first, min_ij.second));
1311 size_t i_max, j_max;
1312 if (abs_max > abs_min)
1313 {
1314 i_max = max_ij.first;
1315 j_max = max_ij.second;
1316 }
1317 else
1318 {
1319 i_max = min_ij.first;
1320 j_max = min_ij.second;
1321 }
1322 l1.setCoefficients (M(i_max,0), M(i_max,1), M(i_max,2));
1323 l2.setCoefficients (M(0, j_max), M(1,j_max), M(2,j_max));
1324
1325 return true;
1326 }
1327
1328
1329 /*
1330 * Return the rectangle that bound the conic section arc characterized by
1331 * the passed points.
1332 *
1333 * P1: the initial point of the arc
1334 * Q: the inner point of the arc
1335 * P2: the final point of the arc
1336 *
1337 * prerequisite: the passed points must lie on the conic
1338 */
arc_bound(const Point & P1,const Point & Q,const Point & P2) const1339 Rect xAx::arc_bound (const Point & P1, const Point & Q, const Point & P2) const
1340 {
1341 using std::swap;
1342 //std::cout << "BOUND: P1 = " << P1 << std::endl;
1343 //std::cout << "BOUND: Q = " << Q << std::endl;
1344 //std::cout << "BOUND: P2 = " << P2 << std::endl;
1345
1346 Rect B(P1, P2);
1347 double Qside = signed_triangle_area (P1, Q, P2);
1348 //std::cout << "BOUND: Qside = " << Qside << std::endl;
1349
1350 Line gl[2];
1351 bool empty[2] = {false, false};
1352
1353 try // if the passed coefficients lead to an equation 0x + 0y + c == 0,
1354 { // with c != 0 the setCoefficients rise an exception
1355 gl[0].setCoefficients (coeff(1), 2 * coeff(2), coeff(4));
1356 }
1357 catch(Geom::LogicalError const &e)
1358 {
1359 empty[0] = true;
1360 }
1361
1362 try
1363 {
1364 gl[1].setCoefficients (2 * coeff(0), coeff(1), coeff(3));
1365 }
1366 catch(Geom::LogicalError const &e)
1367 {
1368 empty[1] = true;
1369 }
1370
1371 std::vector<double> rts;
1372 std::vector<Point> M;
1373 for (size_t dim = 0; dim < 2; ++dim)
1374 {
1375 if (empty[dim]) continue;
1376 rts = roots (gl[dim]);
1377 M.clear();
1378 for (double rt : rts)
1379 M.push_back (gl[dim].pointAt (rt));
1380 if (M.size() == 1)
1381 {
1382 double Mside = signed_triangle_area (P1, M[0], P2);
1383 if (sgn(Mside) == sgn(Qside))
1384 {
1385 //std::cout << "BOUND: M.size() == 1" << std::endl;
1386 B[dim].expandTo(M[0][dim]);
1387 }
1388 }
1389 else if (M.size() == 2)
1390 {
1391 //std::cout << "BOUND: M.size() == 2" << std::endl;
1392 if (M[0][dim] > M[1][dim])
1393 swap (M[0], M[1]);
1394
1395 if (M[0][dim] > B[dim].max())
1396 {
1397 double Mside = signed_triangle_area (P1, M[0], P2);
1398 if (sgn(Mside) == sgn(Qside))
1399 B[dim].setMax(M[0][dim]);
1400 }
1401 else if (M[1][dim] < B[dim].min())
1402 {
1403 double Mside = signed_triangle_area (P1, M[1], P2);
1404 if (sgn(Mside) == sgn(Qside))
1405 B[dim].setMin(M[1][dim]);
1406 }
1407 else
1408 {
1409 double Mside = signed_triangle_area (P1, M[0], P2);
1410 if (sgn(Mside) == sgn(Qside))
1411 B[dim].setMin(M[0][dim]);
1412 Mside = signed_triangle_area (P1, M[1], P2);
1413 if (sgn(Mside) == sgn(Qside))
1414 B[dim].setMax(M[1][dim]);
1415 }
1416 }
1417 }
1418
1419 return B;
1420 }
1421
1422 /*
1423 * Return all points on the conic section nearest to the passed point "P".
1424 *
1425 * P: the point to compute the nearest one
1426 */
allNearestTimes(const Point & P) const1427 std::vector<Point> xAx::allNearestTimes (const Point &P) const
1428 {
1429 // TODO: manage the circle - centre case
1430 std::vector<Point> points;
1431
1432 // named C the conic we look for points (x,y) on C such that
1433 // dot (grad (C(x,y)), rot90 (P -(x,y))) == 0; the set of points satisfying
1434 // this equation is still a conic G, so the wanted points can be found by
1435 // intersecting C with G
1436 xAx G (-coeff(1),
1437 2 * (coeff(0) - coeff(2)),
1438 coeff(1),
1439 -coeff(4) + coeff(1) * P[X] - 2 * coeff(0) * P[Y],
1440 coeff(3) - coeff(1) * P[Y] + 2 * coeff(2) * P[X],
1441 -coeff(3) * P[Y] + coeff(4) * P[X]);
1442
1443 std::vector<Point> crs = intersect (*this, G);
1444
1445 //std::cout << "NEAREST POINT: crs.size = " << crs.size() << std::endl;
1446 if (crs.empty()) return points;
1447
1448 size_t idx = 0;
1449 double mindist = distanceSq (crs[0], P);
1450 std::vector<double> dist;
1451 dist.push_back (mindist);
1452
1453 for (size_t i = 1; i < crs.size(); ++i)
1454 {
1455 dist.push_back (distanceSq (crs[i], P));
1456 if (mindist > dist.back())
1457 {
1458 idx = i;
1459 mindist = dist.back();
1460 }
1461 }
1462
1463 points.push_back (crs[idx]);
1464 for (size_t i = 0; i < crs.size(); ++i)
1465 {
1466 if (i == idx) continue;
1467 if (dist[i] == mindist)
1468 points.push_back (crs[i]);
1469 }
1470
1471 return points;
1472 }
1473
1474
1475
clip(std::vector<RatQuad> & rq,const xAx & cs,const Rect & R)1476 bool clip (std::vector<RatQuad> & rq, const xAx & cs, const Rect & R)
1477 {
1478 clipper aclipper (cs, R);
1479 return aclipper.clip (rq);
1480 }
1481
1482
1483 } // end namespace Geom
1484
1485
1486
1487
1488 /*
1489 Local Variables:
1490 mode:c++
1491 c-file-style:"stroustrup"
1492 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
1493 indent-tabs-mode:nil
1494 fill-column:99
1495 End:
1496 */
1497 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
1498