GCC Code Coverage Report

Line	Branch	Exec	Source
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
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
56		✗	static double det(Point a, Point b) {
57		✗	return a[0]*b[1] - a[1]*b[0];
58			}
59
60			template <typename T>
61		✗	static T det(T a, T b, T c, T d) {
62		✗	return a*d - b*c;
63			}
64
65			template <typename T>
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>
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
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:
86			BadConversion(const std::string& s)
87			: std::runtime_error(s)
88			{ }
89			};
90
91			template <typename T>
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
111		✗	double RatQuad::lambda() const {
112		✗	return 2*(6*w*w +1 -std::sqrt(3*w*w+1))/(12*w*w+3);
113			}
114
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
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
156		✗	CubicBezier RatQuad::toCubic() const {
157		✗	return toCubic(lambda());
158			}
159
160		✗	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
167		✗	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
175		✗	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
185		✗	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
197		✗	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
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
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
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
409		✗	xAx xAx::fromPoint(Point p) {
410		✗	return xAx(1., 0, 1., -2*p[0], -2*p[1], dot(p,p));
411			}
412
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
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
421		✗	xAx xAx::fromLine(Line l) {
422		✗	double dist;
423		✗	Point norm = l.normalAndDist(dist);
424
425		✗	return fromLine(norm, dist);
426			}
427
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
450		✗	double xAx::valueAt(Point P) const {
451		✗	return evaluate_at(P[0], P[1]);
452			}
453
454		✗	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
459		✗	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
466		✗	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			}
473		✗	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			}
480		✗	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
489		✗	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
497		✗	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
511		✗	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
565		✗	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
608		✗	std::vector<double> xAx::roots(Line const &l) const {
609		✗	return roots(l.versor(), l.origin());
610			}
611
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
620		✗	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
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
640		✗	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
648		✗	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
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
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			*/
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			*/
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			*/
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			*/
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			*/
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			*/
940		✗	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			*/
1040		✗	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			*/
1081		✗	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			*/
1152		✗	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			*/
1191		✗	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			*/
1222		✗	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			*/
1254		✗	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		13	std::array<Line, 2> xAx::decompose_df(Coord epsilon) const
1329			{
1330			// For the classification of degenerate conics, see https://mathworld.wolfram.com/QuadraticCurve.html
1331			using std::sqrt, std::abs;
1332
1333			// Create 2 degenerate lines
1334		13	auto const origin = Point(0, 0);
1335		13	std::array<Line, 2> result = {Line(origin, origin), Line(origin, origin)};
1336
1337		13	double A = c[0];
1338		13	double B = c[1];
1339		13	double C = c[2];
1340		13	double D = c[3];
1341		13	double E = c[4];
1342		13	double F = c[5];
1343		13	Coord discriminant = sqr(B) - 4 * A * C;
1344	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 13 times.	13	if (discriminant < -epsilon) {
1345		✗	return result;
1346			}
1347
1348		13	bool single_line = false; // In the generic case, there will be 2 lines.
1349		13	bool parallel_lines = false;
1350	2/2 ✓ Branch 0 taken 10 times. ✓ Branch 1 taken 3 times.	13	if (discriminant < epsilon) {
1351		10	discriminant = 0;
1352		10	parallel_lines = true;
1353			// Check the secondary discriminant
1354		10	Coord const secondary = sqr(D) + sqr(E) - 4 * F * (A + C);
1355	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 10 times.	10	if (secondary < -epsilon) {
1356		✗	return result;
1357			}
1358		10	single_line = (secondary < epsilon);
1359			}
1360
1361	5/6 ✓ Branch 1 taken 5 times. ✓ Branch 2 taken 8 times. ✗ Branch 4 not taken. ✓ Branch 5 taken 5 times. ✓ Branch 6 taken 8 times. ✓ Branch 7 taken 5 times.	13	if (abs(A) > epsilon || abs(C) > epsilon) {
1362			// This is the typical case: either x² or y² come with a nonzero coefficient.
1363			// To guard against numerical errors, we check which of the coefficients A, C has larger absolute value.
1364
1365		8	bool const swap_xy = abs(C) > abs(A);
1366	2/2 ✓ Branch 0 taken 2 times. ✓ Branch 1 taken 6 times.	8	if (swap_xy) {
1367		2	std::swap(A, C);
1368		2	std::swap(D, E);
1369			}
1370
1371			// From now on, we may assume that A is "reasonably large".
1372	2/2 ✓ Branch 0 taken 5 times. ✓ Branch 1 taken 3 times.	8	if (parallel_lines) {
1373	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5 times.	5	if (single_line) {
1374			// Special case: a single line.
1375		✗	std::array<double, 3> coeffs = {sqrt(abs(A)), sqrt(abs(C)), sqrt(abs(F))};
1376		✗	if (swap_xy) {
1377		✗	std::swap(coeffs[0], coeffs[1]);
1378			}
1379		✗	rescale_homogenous(coeffs);
1380		✗	result[0].setCoefficients(coeffs[0], coeffs[1], coeffs[2]);
1381		✗	return result;
1382			}
1383
1384			// Two parallel lines.
1385		5	Coord const quotient_discriminant = sqr(D) - 4 * A * F;
1386	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5 times.	5	if (quotient_discriminant < 0) {
1387		✗	return result;
1388			}
1389		5	Coord const sqrt_disc = sqrt(quotient_discriminant);
1390		5	double const c1 = 0.5 * (D - sqrt_disc);
1391		5	double const c2 = c1 + sqrt_disc;
1392		5	std::array<double, 3> coeffs = {A, 0.5 * B, c1};
1393	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5 times.	5	if (swap_xy) {
1394		✗	std::swap(coeffs[0], coeffs[1]);
1395			}
1396		5	rescale_homogenous(coeffs);
1397	1/2 ✓ Branch 5 taken 5 times. ✗ Branch 6 not taken.	5	result[0].setCoefficients(coeffs[0], coeffs[1], coeffs[2]);
1398
1399		5	coeffs = {A, 0.5 * B, c2};
1400	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5 times.	5	if (swap_xy) {
1401		✗	std::swap(coeffs[0], coeffs[1]);
1402			}
1403		5	rescale_homogenous(coeffs);
1404	1/2 ✓ Branch 5 taken 5 times. ✗ Branch 6 not taken.	5	result[1].setCoefficients(coeffs[0], coeffs[1], coeffs[2]);
1405		5	return result;
1406			}
1407
1408			// Now for the typical case of 2 non-parallel lines.
1409
1410			// We know that A is further away from 0 than C is.
1411			// The mathematical derivation of the solution is as follows:
1412			// let Δ = B² - 4AC (the discriminant); we know Δ > 0.
1413			// Write δ = sqrt(Δ); we know that this is also positive.
1414			// Then the product AΔ is nonzero, so the equation
1415			// Ax² + Bxy + Cy² + Dx + Ey + F = 0
1416			// is equivalent to
1417			// AΔ (Ax² + Bxy + Cy² + Dx + Ey + F) = 0.
1418			// Consider the two factors
1419			// L_1 = Aδx + 0.5 (Bδ-Δ)y + EA - 0.5 D(B-δ)
1420			// L_2 = Aδx + 0.5 (Bδ+Δ)y - EA + 0.5 D(B+δ)
1421			// With a bit of algebra, you can show that L_1 * L_2 expands
1422			// to AΔ (Ax² + Bxy + Cy² + Dx + Ey + F) (in order to get the
1423			// correct value of F, you have to use the fact that the conic
1424			// is degenerate). Therefore, the factors L_1 and L_2 are in
1425			// fact equations of the two lines to be found.
1426		3	Coord const delta = sqrt(discriminant);
1427		3	std::array<double, 3> coeffs1 = {A * delta, 0.5 * (B * delta - discriminant), E * A - 0.5 * D * (B - delta)};
1428		3	std::array<double, 3> coeffs2 = {coeffs1[0], coeffs1[1] + discriminant, D * delta - coeffs1[2]};
1429	2/2 ✓ Branch 0 taken 2 times. ✓ Branch 1 taken 1 times.	3	if (swap_xy) { // We must unswap the coefficients of x and y
1430		2	std::swap(coeffs1[0], coeffs1[1]);
1431		2	std::swap(coeffs2[0], coeffs2[1]);
1432			}
1433
1434		3	unsigned index = 0;
1435	4/6 ✓ Branch 1 taken 1 times. ✓ Branch 2 taken 2 times. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken. ✓ Branch 6 taken 3 times. ✗ Branch 7 not taken.	3	if (coeffs1[0] != 0 || coeffs1[1] != 0) {
1436		3	rescale_homogenous(coeffs1);
1437	1/2 ✓ Branch 5 taken 3 times. ✗ Branch 6 not taken.	3	result[index++].setCoefficients(coeffs1[0], coeffs1[1], coeffs1[2]);
1438			}
1439	2/6 ✗ Branch 1 not taken. ✓ Branch 2 taken 3 times. ✗ Branch 4 not taken. ✗ Branch 5 not taken. ✓ Branch 6 taken 3 times. ✗ Branch 7 not taken.	3	if (coeffs2[0] != 0 || coeffs2[1] != 0) {
1440		3	rescale_homogenous(coeffs2);
1441	1/2 ✓ Branch 5 taken 3 times. ✗ Branch 6 not taken.	3	result[index].setCoefficients(coeffs2[0], coeffs2[1], coeffs2[2]);
1442			}
1443		3	return result;
1444			}
1445
1446			// If we're here, then A==0 and C==0.
1447	1/2 ✓ Branch 1 taken 5 times. ✗ Branch 2 not taken.	5	if (abs(B) < epsilon) { // A == B == C == 0, so the conic reduces to Dx + Ey + F.
1448	3/4 ✓ Branch 0 taken 1 times. ✓ Branch 1 taken 4 times. ✗ Branch 2 not taken. ✓ Branch 3 taken 1 times.	5	if (D == 0 && E == 0) {
1449		✗	return result;
1450			}
1451		5	std::array<double, 3> coeffs = {D, E, F};
1452		5	rescale_homogenous(coeffs);
1453	1/2 ✓ Branch 5 taken 5 times. ✗ Branch 6 not taken.	5	result[0].setCoefficients(coeffs[0], coeffs[1], coeffs[2]);
1454		5	return result;
1455			}
1456
1457			// OK, so A == C == 0 but B != 0. In other words, the conic has the form
1458			// Bxy + Dx + Ey + F. Since B != 0, the zero set stays the same if we multiply the
1459			// equation by B, which gives us this equation:
1460			// B²xy + BDx + BEy + BF = 0.
1461			// The above factors as (Bx + E)(By + D) = 0.
1462		✗	std::array<double, 2> nonzero_coeffs = {B, E};
1463		✗	rescale_homogenous(nonzero_coeffs);
1464		✗	result[0].setCoefficients(nonzero_coeffs[0], 0, nonzero_coeffs[1]);
1465
1466		✗	nonzero_coeffs = {B, D};
1467		✗	rescale_homogenous(nonzero_coeffs);
1468		✗	result[1].setCoefficients(0, nonzero_coeffs[0], nonzero_coeffs[1]);
1469		✗	return result;
1470			}
1471
1472			/*
1473			* Return the rectangle that bound the conic section arc characterized by
1474			* the passed points.
1475			*
1476			* P1: the initial point of the arc
1477			* Q: the inner point of the arc
1478			* P2: the final point of the arc
1479			*
1480			* prerequisite: the passed points must lie on the conic
1481			*/
1482		✗	Rect xAx::arc_bound (const Point & P1, const Point & Q, const Point & P2) const
1483			{
1484			using std::swap;
1485			//std::cout << "BOUND: P1 = " << P1 << std::endl;
1486			//std::cout << "BOUND: Q = " << Q << std::endl;
1487			//std::cout << "BOUND: P2 = " << P2 << std::endl;
1488
1489		✗	Rect B(P1, P2);
1490		✗	double Qside = signed_triangle_area (P1, Q, P2);
1491			//std::cout << "BOUND: Qside = " << Qside << std::endl;
1492
1493		✗	Line gl[2];
1494		✗	bool empty[2] = {false, false};
1495
1496			try // if the passed coefficients lead to an equation 0x + 0y + c == 0,
1497			{ // with c != 0 the setCoefficients rise an exception
1498		✗	gl[0].setCoefficients (coeff(1), 2 * coeff(2), coeff(4));
1499			}
1500		✗	catch(Geom::LogicalError const &e)
1501			{
1502		✗	empty[0] = true;
1503		✗	}
1504
1505			try
1506			{
1507		✗	gl[1].setCoefficients (2 * coeff(0), coeff(1), coeff(3));
1508			}
1509		✗	catch(Geom::LogicalError const &e)
1510			{
1511		✗	empty[1] = true;
1512		✗	}
1513
1514		✗	std::vector<double> rts;
1515		✗	std::vector<Point> M;
1516		✗	for (size_t dim = 0; dim < 2; ++dim)
1517			{
1518		✗	if (empty[dim]) continue;
1519		✗	rts = roots (gl[dim]);
1520		✗	M.clear();
1521		✗	for (double rt : rts)
1522		✗	M.push_back (gl[dim].pointAt (rt));
1523		✗	if (M.size() == 1)
1524			{
1525		✗	double Mside = signed_triangle_area (P1, M[0], P2);
1526		✗	if (sgn(Mside) == sgn(Qside))
1527			{
1528			//std::cout << "BOUND: M.size() == 1" << std::endl;
1529		✗	B[dim].expandTo(M[0][dim]);
1530			}
1531			}
1532		✗	else if (M.size() == 2)
1533			{
1534			//std::cout << "BOUND: M.size() == 2" << std::endl;
1535		✗	if (M[0][dim] > M[1][dim])
1536		✗	swap (M[0], M[1]);
1537
1538		✗	if (M[0][dim] > B[dim].max())
1539			{
1540		✗	double Mside = signed_triangle_area (P1, M[0], P2);
1541		✗	if (sgn(Mside) == sgn(Qside))
1542		✗	B[dim].setMax(M[0][dim]);
1543			}
1544		✗	else if (M[1][dim] < B[dim].min())
1545			{
1546		✗	double Mside = signed_triangle_area (P1, M[1], P2);
1547		✗	if (sgn(Mside) == sgn(Qside))
1548		✗	B[dim].setMin(M[1][dim]);
1549			}
1550			else
1551			{
1552		✗	double Mside = signed_triangle_area (P1, M[0], P2);
1553		✗	if (sgn(Mside) == sgn(Qside))
1554		✗	B[dim].setMin(M[0][dim]);
1555		✗	Mside = signed_triangle_area (P1, M[1], P2);
1556		✗	if (sgn(Mside) == sgn(Qside))
1557		✗	B[dim].setMax(M[1][dim]);
1558			}
1559			}
1560			}
1561
1562		✗	return B;
1563		✗	}
1564
1565			/*
1566			* Return all points on the conic section nearest to the passed point "P".
1567			*
1568			* P: the point to compute the nearest one
1569			*/
1570		✗	std::vector<Point> xAx::allNearestTimes (const Point &P) const
1571			{
1572			// TODO: manage the circle - centre case
1573		✗	std::vector<Point> points;
1574
1575			// named C the conic we look for points (x,y) on C such that
1576			// dot (grad (C(x,y)), rot90 (P -(x,y))) == 0; the set of points satisfying
1577			// this equation is still a conic G, so the wanted points can be found by
1578			// intersecting C with G
1579		✗	xAx G (-coeff(1),
1580		✗	2 * (coeff(0) - coeff(2)),
1581			coeff(1),
1582		✗	-coeff(4) + coeff(1) * P[X] - 2 * coeff(0) * P[Y],
1583		✗	coeff(3) - coeff(1) * P[Y] + 2 * coeff(2) * P[X],
1584		✗	-coeff(3) * P[Y] + coeff(4) * P[X]);
1585
1586		✗	std::vector<Point> crs = intersect (*this, G);
1587
1588			//std::cout << "NEAREST POINT: crs.size = " << crs.size() << std::endl;
1589		✗	if (crs.empty()) return points;
1590
1591		✗	size_t idx = 0;
1592		✗	double mindist = distanceSq (crs[0], P);
1593		✗	std::vector<double> dist;
1594		✗	dist.push_back (mindist);
1595
1596		✗	for (size_t i = 1; i < crs.size(); ++i)
1597			{
1598		✗	dist.push_back (distanceSq (crs[i], P));
1599		✗	if (mindist > dist.back())
1600			{
1601		✗	idx = i;
1602		✗	mindist = dist.back();
1603			}
1604			}
1605
1606		✗	points.push_back (crs[idx]);
1607		✗	for (size_t i = 0; i < crs.size(); ++i)
1608			{
1609		✗	if (i == idx) continue;
1610		✗	if (dist[i] == mindist)
1611		✗	points.push_back (crs[i]);
1612			}
1613
1614		✗	return points;
1615		✗	}
1616
1617
1618
1619		✗	bool clip (std::vector<RatQuad> & rq, const xAx & cs, const Rect & R)
1620			{
1621		✗	clipper aclipper (cs, R);
1622		✗	return aclipper.clip (rq);
1623		✗	}
1624
1625
1626			} // end namespace Geom
1627
1628
1629
1630
1631			/*
1632			Local Variables:
1633			mode:c++
1634			c-file-style:"stroustrup"
1635			c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
1636			indent-tabs-mode:nil
1637			fill-column:99
1638			End:
1639			*/
1640			// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
1641