GCC Code Coverage Report

Line	Branch	Exec	Source
1			/** @file
2			* @brief Conic Section
3			*//*
4			* Authors:
5			* Nathan Hurst <njh@njhurst.com>
6			*
7			* Copyright 2009 authors
8			*
9			* This library is free software; you can redistribute it and/or
10			* modify it either under the terms of the GNU Lesser General Public
11			* License version 2.1 as published by the Free Software Foundation
12			* (the "LGPL") or, at your option, under the terms of the Mozilla
13			* Public License Version 1.1 (the "MPL"). If you do not alter this
14			* notice, a recipient may use your version of this file under either
15			* the MPL or the LGPL.
16			*
17			* You should have received a copy of the LGPL along with this library
18			* in the file COPYING-LGPL-2.1; if not, write to the Free Software
19			* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20			* You should have received a copy of the MPL along with this library
21			* in the file COPYING-MPL-1.1
22			*
23			* The contents of this file are subject to the Mozilla Public License
24			* Version 1.1 (the "License"); you may not use this file except in
25			* compliance with the License. You may obtain a copy of the License at
26			* http://www.mozilla.org/MPL/
27			*
28			* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
29			* OF ANY KIND, either express or implied. See the LGPL or the MPL for
30			* the specific language governing rights and limitations.
31			*/
32
33
34			#ifndef LIB2GEOM_SEEN_CONICSEC_H
35			#define LIB2GEOM_SEEN_CONICSEC_H
36
37			#include <2geom/exception.h>
38			#include <2geom/angle.h>
39			#include <2geom/rect.h>
40			#include <2geom/affine.h>
41			#include <2geom/point.h>
42			#include <2geom/line.h>
43			#include <2geom/bezier-curve.h>
44			#include <2geom/numeric/linear_system.h>
45			#include <2geom/numeric/symmetric-matrix-fs.h>
46			#include <2geom/numeric/symmetric-matrix-fs-operation.h>
47
48			#include <optional>
49
50			#include <string>
51			#include <vector>
52			#include <ostream>
53
54
55
56
57			namespace Geom
58			{
59
60			class RatQuad{
61			/**
62			* A curve of the form B02*A + B12*B*w + B22*C/(B02 + B12*w + B22)
63			* These curves can exactly represent a piece conic section less than a certain angle (find out)
64			*
65			**/
66
67			public:
68			Point P[3];
69			double w;
70		✗	RatQuad() {}
71		✗	RatQuad(Point a, Point b, Point c, double w) : w(w) {
72		✗	P[0] = a;
73		✗	P[1] = b;
74		✗	P[2] = c;
75		✗	}
76			double lambda() const;
77
78			static RatQuad fromPointsTangents(Point P0, Point dP0,
79			Point P,
80			Point P2, Point dP2);
81			static RatQuad circularArc(Point P0, Point P1, Point P2);
82
83			CubicBezier toCubic() const;
84			CubicBezier toCubic(double lam) const;
85
86			Point pointAt(double t) const;
87		✗	Point at0() const {return P[0];}
88		✗	Point at1() const {return P[2];}
89
90			void split(RatQuad &a, RatQuad &b) const;
91
92			D2<SBasis> hermite() const;
93			std::vector<SBasis> homogeneous() const;
94			};
95
96
97
98
99			class xAx{
100			public:
101			double c[6];
102
103			enum kind_t
104			{
105			PARABOLA,
106			CIRCLE,
107			REAL_ELLIPSE,
108			IMAGINARY_ELLIPSE,
109			RECTANGULAR_HYPERBOLA,
110			HYPERBOLA,
111			DOUBLE_LINE,
112			TWO_REAL_PARALLEL_LINES,
113			TWO_IMAGINARY_PARALLEL_LINES,
114			TWO_REAL_CROSSING_LINES,
115			TWO_IMAGINARY_CROSSING_LINES, // crossing at a real point
116			SINGLE_POINT = TWO_IMAGINARY_CROSSING_LINES,
117			UNKNOWN
118			};
119
120
121		✗	xAx() {}
122
123			/*
124			* Define the conic section by its algebraic equation coefficients
125			*
126			* c0, .., c5: equation coefficients
127			*/
128		13	xAx (double c0, double c1, double c2, double c3, double c4, double c5)
129		13	{
130		13	set (c0, c1, c2, c3, c4, c5);
131		13	}
132
133			/*
134			* Define a conic section by its related symmetric matrix
135			*/
136			xAx (const NL::ConstSymmetricMatrixView<3> & C)
137			{
138			set(C);
139			}
140
141			/*
142			* Define a conic section by computing the one that fits better with
143			* N points.
144			*
145			* points: points to fit
146			*
147			* precondition: there must be at least 5 non-overlapping points
148			*/
149			xAx (std::vector<Point> const& points)
150			{
151			set (points);
152			}
153
154			/*
155			* Define a section conic by providing the coordinates of one of its
156			* vertex,the major axis inclination angle and the coordinates of its foci
157			* with respect to the unidimensional system defined by the major axis with
158			* origin set at the provided vertex.
159			*
160			* _vertex : section conic vertex V
161			* _angle : section conic major axis angle
162			* _dist1: +/-distance btw V and nearest focus
163			* _dist2: +/-distance btw V and farest focus
164			*
165			* prerequisite: _dist1 <= _dist2
166			*/
167			xAx (const Point& _vertex, double _angle, double _dist1, double _dist2)
168			{
169			set (_vertex, _angle, _dist1, _dist2);
170			}
171
172			/*
173			* Define a conic section by providing one of its vertex and its foci.
174			*
175			* _vertex: section conic vertex
176			* _focus1: section conic focus
177			* _focus2: section conic focus
178			*/
179			xAx (const Point& _vertex, const Point& _focus1, const Point& _focus2)
180			{
181			set(_vertex, _focus1, _focus2);
182			}
183
184			/*
185			* Define a conic section by passing a focus, the related directrix,
186			* and the eccentricity (e)
187			* (e < 1 -> ellipse; e = 1 -> parabola; e > 1 -> hyperbola)
188			*
189			* _focus: a focus of the conic section
190			* _directrix: the directrix related to the given focus
191			* _eccentricity: the eccentricity parameter of the conic section
192			*/
193			xAx (const Point & _focus, const Line & _directrix, double _eccentricity)
194			{
195			set (_focus, _directrix, _eccentricity);
196			}
197
198			/*
199			* Made up a degenerate conic section as a pair of lines
200			*
201			* l1, l2: lines that made up the conic section
202			*/
203			xAx (const Line& l1, const Line& l2)
204			{
205			set (l1, l2);
206			}
207
208			/*
209			* Define the conic section by its algebraic equation coefficients
210			* c0, ..., c5: equation coefficients
211			*/
212		13	void set (double c0, double c1, double c2, double c3, double c4, double c5)
213			{
214		13	c[0] = c0; c[1] = c1; c[2] = c2; // xx, xy, yy
215		13	c[3] = c3; c[4] = c4; // x, y
216		13	c[5] = c5; // 1
217		13	}
218
219			/*
220			* Define a conic section by its related symmetric matrix
221			*/
222			void set (const NL::ConstSymmetricMatrixView<3> & C)
223			{
224			set(C(0,0), 2*C(1,0), C(1,1), 2*C(2,0), 2*C(2,1), C(2,2));
225			}
226
227			void set (std::vector<Point> const& points);
228
229			void set (const Point& _vertex, double _angle, double _dist1, double _dist2);
230
231			void set (const Point& _vertex, const Point& _focus1, const Point& _focus2);
232
233			void set (const Point & _focus, const Line & _directrix, double _eccentricity);
234
235			void set (const Line& l1, const Line& l2);
236
237
238			static xAx fromPoint(Point p);
239			static xAx fromDistPoint(Point p, double d);
240			static xAx fromLine(Point n, double d);
241			static xAx fromLine(Line l);
242			static xAx fromPoints(std::vector<Point> const &pts);
243
244
245			template<typename T>
246		✗	T evaluate_at(T x, T y) const {
247		✗	return c[0]*x*x + c[1]*x*y + c[2]*y*y + c[3]*x + c[4]*y + c[5];
248			}
249
250			double valueAt(Point P) const;
251
252		✗	std::vector<double> implicit_form_coefficients() const {
253		✗	return std::vector<double>(c, c+6);
254			}
255
256			template<typename T>
257		✗	T evaluate_at(T x, T y, T w) const {
258		✗	return c[0]*x*x + c[1]*x*y + c[2]*y*y + c[3]*x*w + c[4]*y*w + c[5]*w*w;
259			}
260
261			xAx scale(double sx, double sy) const;
262
263			Point gradient(Point p) const;
264
265			xAx operator-(xAx const &b) const;
266			xAx operator+(xAx const &b) const;
267			xAx operator+(double const &b) const;
268			xAx operator*(double const &b) const;
269
270			std::vector<Point> crossings(Rect r) const;
271			std::optional<RatQuad> toCurve(Rect const & bnd) const;
272			std::vector<double> roots(Point d, Point o) const;
273
274			std::vector<double> roots(Line const &l) const;
275
276			static Interval quad_ex(double a, double b, double c, Interval ivl);
277
278			Geom::Affine hessian() const;
279
280			std::optional<Point> bottom() const;
281
282			Interval extrema(Rect r) const;
283
284
285			/*
286			* Return the symmetric matrix related to the conic section.
287			* Modifying the matrix does not modify the conic section
288			*/
289		✗	NL::SymmetricMatrix<3> get_matrix() const
290			{
291		✗	NL::SymmetricMatrix<3> C(c);
292		✗	C(1,0) *= 0.5; C(2,0) *= 0.5; C(2,1) *= 0.5;
293		✗	return C;
294		✗	}
295
296			/*
297			* Return the i-th coefficient of the conic section algebraic equation
298			* Modifying the returned value does not modify the conic section coefficient
299			*/
300		✗	double coeff (size_t i) const
301			{
302		✗	return c[i];
303			}
304
305			/*
306			* Return the i-th coefficient of the conic section algebraic equation
307			* Modifying the returned value modifies the conic section coefficient
308			*/
309		✗	double& coeff (size_t i)
310			{
311		✗	return c[i];
312			}
313
314			kind_t kind () const;
315
316			std::string categorise() const;
317
318			/*
319			* Return true if the equation:
320			* c0*x^2 + c1*xy + c2*y^2 + c3*x + c4*y +c5 == 0
321			* really defines a conic, false otherwise
322			*/
323		✗	bool is_quadratic() const
324			{
325		✗	return (coeff(0) != 0 || coeff(1) != 0 || coeff(2) != 0);
326			}
327
328			/*
329			* Return true if the conic is degenerate, i.e. if the related matrix
330			* determinant is null, false otherwise
331			*/
332		✗	bool isDegenerate() const
333			{
334		✗	return (det_sgn (get_matrix()) == 0);
335			}
336
337			/*
338			* Compute the centre of symmetry of the conic section when it exists,
339			* else it return an uninitialized std::optional<Point> instance.
340			*/
341		✗	std::optional<Point> centre() const
342			{
343			typedef std::optional<Point> opt_point_t;
344
345		✗	double d = coeff(1) * coeff(1) - 4 * coeff(0) * coeff(2);
346		✗	if (are_near (d, 0)) return opt_point_t();
347		✗	NL::Matrix Q(2, 2);
348		✗	Q(0,0) = coeff(0);
349		✗	Q(1,1) = coeff(2);
350		✗	Q(0,1) = Q(1,0) = coeff(1) * 0.5;
351		✗	NL::Vector T(2);
352		✗	T[0] = - coeff(3) * 0.5;
353		✗	T[1] = - coeff(4) * 0.5;
354
355		✗	NL::LinearSystem ls (Q, T);
356		✗	NL::Vector sol = ls.SV_solve();
357		✗	Point C;
358		✗	C[0] = sol[0];
359		✗	C[1] = sol[1];
360
361		✗	return opt_point_t(C);
362		✗	}
363
364			double axis_angle() const;
365
366			void roots (std::vector<double>& sol, Coord v, Dim2 d) const;
367
368			xAx translate (const Point & _offset) const;
369
370			xAx rotate (double angle) const;
371
372			/*
373			* Rotate the conic section by the given angle wrt the provided point.
374			*
375			* _rot_centre: the rotation centre
376			* _angle: the rotation angle
377			*/
378			xAx rotate (const Point & _rot_centre, double _angle) const
379			{
380			xAx result
381			= translate (-_rot_centre).rotate (_angle).translate (_rot_centre);
382			return result;
383			}
384
385			/*
386			* Compute the tangent line of the conic section at the provided point
387			*
388			* _point: the conic section point the tangent line pass through
389			*/
390		✗	Line tangent (const Point & _point) const
391			{
392		✗	NL::Vector pp(3);
393		✗	pp[0] = _point[0]; pp[1] = _point[1]; pp[2] = 1;
394		✗	NL::SymmetricMatrix<3> C = get_matrix();
395		✗	NL::Vector line = C * pp;
396		✗	return Line(line[0], line[1], line[2]);
397		✗	}
398
399			/*
400			* For a non degenerate conic compute the dual conic.
401			* TODO: investigate degenerate case
402			*/
403			xAx dual () const
404			{
405			//assert (! isDegenerate());
406			NL::SymmetricMatrix<3> C = get_matrix();
407			NL::SymmetricMatrix<3> D = adj(C);
408			xAx dc(D);
409			return dc;
410			}
411
412			bool decompose (Line& l1, Line& l2) const;
413
414			/**
415			* @brief Division-free decomposition of a degenerate conic section, without
416			* degeneration test.
417			*
418			* When the conic is degenerate, it consists of 0, 1 or 2 lines in the xy-plane.
419			* This function returns these lines. But it does not check whether the conic
420			* is really degenerate, so calling it on a non-degenerate conic produces a
421			* meaningless result.
422			*
423			* If the number of lines is less than two, the trailing lines in the returned
424			* array will be degenerate. Use Line::isDegenerate() to test for that.
425			*
426			* This version of the decomposition is division-free, which improves numerical
427			* stability compared to decompose().
428			*
429			* @param epsilon The numerical threshold for floating point comparison of discriminants.
430			*/
431			std::array<Line, 2> decompose_df(Coord epsilon = EPSILON) const;
432
433			/*
434			* Generate a RatQuad object from a conic arc.
435			*
436			* p0: the initial point of the arc
437			* p1: the inner point of the arc
438			* p2: the final point of the arc
439			*/
440		✗	RatQuad toRatQuad (const Point & p0,
441			const Point & p1,
442			const Point & p2) const
443			{
444		✗	Point dp0 = gradient (p0);
445		✗	Point dp2 = gradient (p2);
446			return
447		✗	RatQuad::fromPointsTangents (p0, rot90 (dp0), p1, p2, rot90 (dp2));
448			}
449
450			/*
451			* Return the angle related to the normal gradient computed at the passed
452			* point.
453			*
454			* _point: the point at which computes the angle
455			*
456			* prerequisite: the passed point must lie on the conic
457			*/
458		✗	double angle_at (const Point & _point) const
459			{
460		✗	double angle = atan2 (gradient (_point));
461		✗	if (angle < 0) angle += (2*M_PI);
462		✗	return angle;
463			}
464
465			/*
466			* Return true if the given point is contained in the conic arc determined
467			* by the passed points.
468			*
469			* _point: the point to be tested
470			* _initial: the initial point of the arc
471			* _inner: an inner point of the arc
472			* _final: the final point of the arc
473			*
474			* prerequisite: the passed points must lie on the conic, the inner point
475			* has to be strictly contained in the arc, except when the
476			* initial and final points are equal: in such a case if the
477			* inner point is also equal to them, then they define an arc
478			* made up by a single point.
479			*
480			*/
481			bool arc_contains (const Point & _point, const Point & _initial,
482			const Point & _inner, const Point & _final) const
483			{
484			AngleInterval ai(angle_at(_initial), angle_at(_inner), angle_at(_final));
485			return ai.contains(angle_at(_point));
486			}
487
488			Rect arc_bound (const Point & P1, const Point & Q, const Point & P2) const;
489
490			std::vector<Point> allNearestTimes (const Point &P) const;
491
492			/*
493			* Return the point on the conic section nearest to the passed point "P".
494			*
495			* P: the point to compute the nearest one
496			*/
497			Point nearestTime (const Point &P) const
498			{
499			std::vector<Point> points = allNearestTimes (P);
500			if ( !points.empty() )
501			{
502			return points.front();
503			}
504			// else
505			THROW_LOGICALERROR ("nearestTime: no nearest point found");
506			return Point();
507			}
508
509			};
510
511			std::vector<Point> intersect(const xAx & C1, const xAx & C2);
512
513			bool clip (std::vector<RatQuad> & rq, const xAx & cs, const Rect & R);
514
515		✗	inline std::ostream &operator<< (std::ostream &out_file, const xAx &x) {
516		✗	for(double i : x.c) {
517		✗	out_file << i << ", ";
518			}
519		✗	return out_file;
520			}
521
522			};
523
524
525			#endif // LIB2GEOM_SEEN_CONICSEC_H
526
527			/*
528			Local Variables:
529			mode:c++
530			c-file-style:"stroustrup"
531			c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
532			indent-tabs-mode:nil
533			fill-column:99
534			End:
535			*/
536			// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
537
538