1 /*
2 * Ellipse Curve
3 *
4 * Authors:
5 * Marco Cecchetti <mrcekets at gmail.com>
6 *
7 * Copyright 2008 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 */
34 #include <2geom/ellipse.h>
35 #include <2geom/svg-elliptical-arc.h>
36 #include <2geom/numeric/fitting-tool.h>
37 #include <2geom/numeric/fitting-model.h>
40 namespace Geom
41 {
43 void Ellipse::set(double A, double B, double C, double D, double E, double F)
44 {
45 double den = 4*A*C - B*B;
46 if ( den == 0 )
47 {
48 THROW_LOGICALERROR("den == 0, while computing ellipse centre");
49 }
50 m_centre[X] = (B*E - 2*C*D) / den;
51 m_centre[Y] = (B*D - 2*A*E) / den;
53 // evaluate the a coefficient of the ellipse equation in normal form
54 // E(x,y) = a*(x-cx)^2 + b*(x-cx)*(y-cy) + c*(y-cy)^2 = 1
55 // where b = a*B , c = a*C, (cx,cy) == centre
56 double num = A * sqr(m_centre[X])
57 + B * m_centre[X] * m_centre[Y]
58 + C * sqr(m_centre[Y])
59 - F;
62 //evaluate ellipse rotation angle
63 double rot = std::atan2( -B, -(A - C) )/2;
64 // std::cerr << "rot = " << rot << std::endl;
65 bool swap_axes = false;
66 if ( are_near(rot, 0) ) rot = 0;
67 if ( are_near(rot, M_PI/2) || rot < 0 )
68 {
69 swap_axes = true;
70 }
72 // evaluate the length of the ellipse rays
73 double cosrot = std::cos(rot);
74 double sinrot = std::sin(rot);
75 double cos2 = cosrot * cosrot;
76 double sin2 = sinrot * sinrot;
77 double cossin = cosrot * sinrot;
79 den = A * cos2 + B * cossin + C * sin2;
80 if ( den == 0 )
81 {
82 THROW_LOGICALERROR("den == 0, while computing 'rx' coefficient");
83 }
84 double rx2 = num/den;
85 if ( rx2 < 0 )
86 {
87 THROW_LOGICALERROR("rx2 < 0, while computing 'rx' coefficient");
88 }
89 double rx = std::sqrt(rx2);
91 den = C * cos2 - B * cossin + A * sin2;
92 if ( den == 0 )
93 {
94 THROW_LOGICALERROR("den == 0, while computing 'ry' coefficient");
95 }
96 double ry2 = num/den;
97 if ( ry2 < 0 )
98 {
99 THROW_LOGICALERROR("ry2 < 0, while computing 'rx' coefficient");
100 }
101 double ry = std::sqrt(ry2);
103 // the solution is not unique so we choose always the ellipse
104 // with a rotation angle between 0 and PI/2
105 if ( swap_axes ) std::swap(rx, ry);
106 if ( are_near(rot, M_PI/2)
107 || are_near(rot, -M_PI/2)
108 || are_near(rx, ry) )
109 {
110 rot = 0;
111 }
112 else if ( rot < 0 )
113 {
114 rot += M_PI/2;
115 }
117 m_ray[X] = rx;
118 m_ray[Y] = ry;
119 m_angle = rot;
120 }
123 std::vector<double> Ellipse::implicit_form_coefficients() const
124 {
125 if (ray(X) == 0 || ray(Y) == 0)
126 {
127 THROW_LOGICALERROR("a degenerate ellipse doesn't own an implicit form");
128 }
130 std::vector<double> coeff(6);
131 double cosrot = std::cos(rot_angle());
132 double sinrot = std::sin(rot_angle());
133 double cos2 = cosrot * cosrot;
134 double sin2 = sinrot * sinrot;
135 double cossin = cosrot * sinrot;
136 double invrx2 = 1 / (ray(X) * ray(X));
137 double invry2 = 1 / (ray(Y) * ray(Y));
139 coeff[0] = invrx2 * cos2 + invry2 * sin2;
140 coeff[1] = 2 * (invrx2 - invry2) * cossin;
141 coeff[2] = invrx2 * sin2 + invry2 * cos2;
142 coeff[3] = -(2 * coeff[0] * center(X) + coeff[1] * center(Y));
143 coeff[4] = -(2 * coeff[2] * center(Y) + coeff[1] * center(X));
144 coeff[5] = coeff[0] * center(X) * center(X)
145 + coeff[1] * center(X) * center(Y)
146 + coeff[2] * center(Y) * center(Y)
147 - 1;
148 return coeff;
149 }
152 void Ellipse::set(std::vector<Point> const& points)
153 {
154 size_t sz = points.size();
155 if (sz < 5)
156 {
157 THROW_RANGEERROR("fitting error: too few points passed");
158 }
159 NL::LFMEllipse model;
160 NL::least_squeares_fitter<NL::LFMEllipse> fitter(model, sz);
162 for (size_t i = 0; i < sz; ++i)
163 {
164 fitter.append(points[i]);
165 }
166 fitter.update();
168 NL::Vector z(sz, 0.0);
169 model.instance(*this, fitter.result(z));
170 }
173 SVGEllipticalArc
174 Ellipse::arc(Point const& initial, Point const& inner, Point const& final,
175 bool _svg_compliant)
176 {
177 Point sp_cp = initial - center();
178 Point ep_cp = final - center();
179 Point ip_cp = inner - center();
181 double angle1 = angle_between(sp_cp, ep_cp);
182 double angle2 = angle_between(sp_cp, ip_cp);
183 double angle3 = angle_between(ip_cp, ep_cp);
185 bool large_arc_flag = true;
186 bool sweep_flag = true;
188 if ( angle1 > 0 )
189 {
190 if ( angle2 > 0 && angle3 > 0 )
191 {
192 large_arc_flag = false;
193 sweep_flag = true;
194 }
195 else
196 {
197 large_arc_flag = true;
198 sweep_flag = false;
199 }
200 }
201 else
202 {
203 if ( angle2 < 0 && angle3 < 0 )
204 {
205 large_arc_flag = false;
206 sweep_flag = false;
207 }
208 else
209 {
210 large_arc_flag = true;
211 sweep_flag = true;
212 }
213 }
215 SVGEllipticalArc ea( initial, ray(X), ray(Y), rot_angle(),
216 large_arc_flag, sweep_flag, final, _svg_compliant);
217 return ea;
218 }
220 Ellipse Ellipse::transformed(Matrix const& m) const
221 {
222 double cosrot = std::cos(rot_angle());
223 double sinrot = std::sin(rot_angle());
224 Matrix A( ray(X) * cosrot, ray(X) * sinrot,
225 -ray(Y) * sinrot, ray(Y) * cosrot,
226 0, 0 );
227 Point new_center = center() * m;
228 Matrix M = m.without_translation();
229 Matrix AM = A * M;
230 if ( are_near(AM.det(), 0) )
231 {
232 double angle;
233 if (AM[0] != 0)
234 {
235 angle = std::atan2(AM[2], AM[0]);
236 }
237 else if (AM[1] != 0)
238 {
239 angle = std::atan2(AM[3], AM[1]);
240 }
241 else
242 {
243 angle = M_PI/2;
244 }
245 Point V(std::cos(angle), std::sin(angle));
246 V *= AM;
247 double rx = L2(V);
248 angle = atan2(V);
249 return Ellipse(new_center[X], new_center[Y], rx, 0, angle);
250 }
252 std::vector<double> coeff = implicit_form_coefficients();
253 Matrix Q( coeff[0], coeff[1]/2,
254 coeff[1]/2, coeff[2],
255 0, 0 );
257 Matrix invm = M.inverse();
258 Q = invm * Q ;
259 std::swap( invm[1], invm[2] );
260 Q *= invm;
261 Ellipse e(Q[0], 2*Q[1], Q[3], 0, 0, -1);
262 e.m_centre = new_center;
264 return e;
265 }
267 Ellipse::Ellipse(Geom::Circle const &c)
268 {
269 m_centre = c.center();
270 m_ray[X] = m_ray[Y] = c.ray();
271 }
273 } // end namespace Geom
276 /*
277 Local Variables:
278 mode:c++
279 c-file-style:"stroustrup"
280 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
281 indent-tabs-mode:nil
282 fill-column:99
283 End:
284 */
285 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :