1 /*
2 * Show off crossings between two D2<SBasis> curves.
3 * The intersection points are found by using implicitization tecnique.
4 *
5 * Authors:
6 * Marco Cecchetti <mrcekets at gmail.com>
7 *
8 * Copyright 2008 authors
9 *
10 * This library is free software; you can redistribute it and/or
11 * modify it either under the terms of the GNU Lesser General Public
12 * License version 2.1 as published by the Free Software Foundation
13 * (the "LGPL") or, at your option, under the terms of the Mozilla
14 * Public License Version 1.1 (the "MPL"). If you do not alter this
15 * notice, a recipient may use your version of this file under either
16 * the MPL or the LGPL.
17 *
18 * You should have received a copy of the LGPL along with this library
19 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 * You should have received a copy of the MPL along with this library
22 * in the file COPYING-MPL-1.1
23 *
24 * The contents of this file are subject to the Mozilla Public License
25 * Version 1.1 (the "License"); you may not use this file except in
26 * compliance with the License. You may obtain a copy of the License at
27 * http://www.mozilla.org/MPL/
28 *
29 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
30 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
31 * the specific language governing rights and limitations.
32 */
33
34
35 #include <toys/path-cairo.h>
36 #include <toys/toy-framework-2.h>
37 #include <2geom/d2.h>
38 #include <2geom/sbasis-poly.h>
39 #include <2geom/numeric/linear_system.h>
40 #include <2geom/symbolic/implicit.h>
41
42
43 using namespace Geom;
44
45 /*
46 * helper routines
47 */
poly_to_mvpoly1(SL::MVPoly1 & p,Geom::Poly const & q)48 void poly_to_mvpoly1(SL::MVPoly1 & p, Geom::Poly const& q)
49 {
50 for (size_t i = 0; i < q.size(); ++i)
51 {
52 p.coefficient(i, q[i]);
53 }
54 p.normalize();
55 }
56
mvpoly1_to_poly(Geom::Poly & p,SL::MVPoly1 const & q)57 void mvpoly1_to_poly(Geom::Poly & p, SL::MVPoly1 const& q)
58 {
59 p.resize(q.get_poly().size());
60 for (size_t i = 0; i < q.get_poly().size(); ++i)
61 {
62 p[i] = q[i];
63 }
64 }
65
66
67 /*
68 * intersection_info
69 * structure utilized to store intersection info
70 *
71 * p - the intersection point
72 * t0 - the parameter t value at which the first curve pass through p
73 * t1 - the parameter t value at which the first curve pass through p
74 */
75 struct intersection_info
76 {
intersection_infointersection_info77 intersection_info()
78 {}
79
intersection_infointersection_info80 intersection_info(Point const& _p, Coord _t0, Coord _t1)
81 : p(_p), t0(_t0), t1(_t1)
82 {}
83
84 Point p;
85 Coord t0, t1;
86 };
87
88 typedef std::vector<intersection_info> intersections_info;
89
90
91
92 /*
93 * intersection algorithm
94 */
intersect(intersections_info & xs,D2<SBasis> const & A,D2<SBasis> const & B)95 void intersect(intersections_info& xs, D2<SBasis> const& A, D2<SBasis> const& B)
96 {
97 using std::swap;
98
99 // supposing implicitization the most expensive step
100 // we perform a call to intersect with curve arguments swapped
101 if (A[0].size() > B[0].size())
102 {
103 intersect(xs, B, A);
104 for (size_t i = 0; i < xs.size(); ++i)
105 swap(xs[i].t0, xs[i].t1);
106
107 return;
108 }
109
110 // convert A from symmetric power basis to power basis
111 Geom::Poly A0 = sbasis_to_poly(A[0]);
112 Geom::Poly A1 = sbasis_to_poly(A[1]);
113
114 // convert to MultiPoly type
115 SL::MVPoly1 Af, Ag;
116 poly_to_mvpoly1(Af, A0);
117 poly_to_mvpoly1(Ag, A1);
118
119 // compute a basis of the ideal related to the curve A
120 // in vector form
121 Geom::SL::basis_type b;
122 // if we compute the micro-basis the bezout matrix is made up
123 // by one only entry so we can't do the inversion step.
124 if (A0.size() == 3)
125 {
126 make_initial_basis(b, Af, Ag);
127 }
128 else
129 {
130 microbasis(b, Af, Ag);
131 }
132
133 // we put the basis in of the form of two independent moving line
134 Geom::SL::MVPoly3 p, q;
135 basis_to_poly(p, b[0]);
136 basis_to_poly(q, b[1]);
137
138 // compute the Bezout matrix and the implicit equation of the curve A
139 Geom::SL::Matrix<Geom::SL::MVPoly2> BZ = make_bezout_matrix(p, q);
140 SL::MVPoly2 ic = determinant_minor(BZ);
141 ic.normalize();
142
143
144 // convert B from symmetric power basis to power basis
145 Geom::Poly B0 = sbasis_to_poly(B[0]);
146 Geom::Poly B1 = sbasis_to_poly(B[1]);
147
148 // convert to MultiPoly type
149 SL::MVPoly1 Bf, Bg;
150 poly_to_mvpoly1(Bf, B0);
151 poly_to_mvpoly1(Bg, B1);
152
153 // evaluate the implicit equation of A on B
154 // so we get an s(t) polynomial that give us
155 // the t values for B at which intersection happens
156 SL::MVPoly1 s = ic(Bf, Bg);
157
158 // convert s(t) to Poly type, in order to use the real_solve function
159 Geom::Poly z;
160 mvpoly1_to_poly(z, s);
161
162 // compute t values for the curve B at which intersection happens
163 std::vector<double> sol = solve_reals(z);
164
165 // filter the found solutions wrt the domain interval [0,1] of B
166 // and compute the related point coordinates
167 std::vector<double> pt;
168 pt.reserve(sol.size());
169 std::vector<Point> points;
170 points.reserve(sol.size());
171 for (size_t i = 0; i < sol.size(); ++i)
172 {
173 if (sol[i] >= 0 && sol[i] <= 1)
174 {
175 pt.push_back(sol[i]);
176 points.push_back(B(pt.back()));
177 }
178 }
179
180 // case: A is parametrized by polynomial of degree 1
181 // we compute the t values of A at the intersection points
182 // and filter the results wrt the domain interval [0,1]
183 double t;
184 xs.clear();
185 xs.reserve(pt.size());
186 if (A0.size() == 2)
187 {
188 for (size_t i = 0; i < points.size(); ++i)
189 {
190 t = (points[i][X] - A0[0]) / A0[1];
191 if (t >= 0 && t <= 1)
192 {
193 xs.push_back(intersection_info(points[i], t, pt[i]));
194 }
195 }
196 return;
197 }
198
199 // general case
200 // we compute the value of the parameter t of A at each intersection point
201 // and we filter the final result wrt the domain interval [0,1]
202 // the computation is performed by using the inversion formula for each point
203 // As reference see:
204 // Sederberger - Computer Aided Geometric Design
205 // par 16.5 - Implicitization and Inversion
206 size_t n = BZ.rows();
207 Geom::NL::Matrix BZN(n, n);
208 Geom::NL::MatrixView BZV(BZN, 0, 0, n-1, n-1);
209 Geom::NL::VectorView cv = BZN.column_view(n-1);
210 Geom::NL::VectorView bv(cv, n-1);
211 Geom::NL::LinearSystem ls(BZV, bv);
212 for (size_t i = 0; i < points.size(); ++i)
213 {
214 // evaluate the first main minor of order n-1 at each intersection point
215 polynomial_matrix_evaluate(BZN, BZ, points[i]);
216 // solve the linear system with the powers of t as unknowns
217 ls.SV_solve();
218 // the last element contains the t value
219 t = -ls.solution()[n-2];
220 // filter with respect to the domain of A
221 if (t >= 0 && t <= 1)
222 {
223 xs.push_back(intersection_info(points[i], t, pt[i]));
224 }
225 }
226 }
227
228
229
230 class IntersectImplicit : public Toy
231 {
232
draw(cairo_t * cr,std::ostringstream * notify,int width,int height,bool save,std::ostringstream * timer_stream)233 void draw( cairo_t *cr, std::ostringstream *notify,
234 int width, int height, bool save, std::ostringstream *timer_stream)
235 {
236 cairo_set_line_width (cr, 0.3);
237 cairo_set_source_rgba (cr, 0.8, 0., 0, 1);
238 D2<SBasis> A = pshA.asBezier();
239 cairo_d2_sb(cr, A);
240 cairo_stroke(cr);
241 cairo_set_source_rgba (cr, 0.0, 0., 0, 1);
242 D2<SBasis> B = pshB.asBezier();
243 cairo_d2_sb(cr, B);
244 cairo_stroke(cr);
245
246 intersect(xs, A, B);
247 for (size_t i = 0; i < xs.size(); ++i)
248 {
249 draw_handle(cr, xs[i].p);
250 }
251
252 Toy::draw(cr, notify, width, height, save,timer_stream);
253 }
254
255
256 public:
IntersectImplicit(unsigned int _A_bez_ord,unsigned int _B_bez_ord)257 IntersectImplicit(unsigned int _A_bez_ord, unsigned int _B_bez_ord)
258 : A_bez_ord(_A_bez_ord), B_bez_ord(_B_bez_ord)
259 {
260 handles.push_back(&pshA);
261 for (unsigned int i = 0; i <= A_bez_ord; ++i)
262 pshA.push_back(Geom::Point(uniform()*400, uniform()*400));
263 handles.push_back(&pshB);
264 for (unsigned int i = 0; i <= B_bez_ord; ++i)
265 pshB.push_back(Geom::Point(uniform()*400, uniform()*400));
266
267 }
268
269 private:
270 unsigned int A_bez_ord, B_bez_ord;
271 PointSetHandle pshA, pshB;
272 intersections_info xs;
273 };
274
275
main(int argc,char ** argv)276 int main(int argc, char **argv)
277 {
278 unsigned int A_bez_ord = 4;
279 unsigned int B_bez_ord = 6;
280 if(argc > 1)
281 sscanf(argv[1], "%d", &A_bez_ord);
282 if(argc > 2)
283 sscanf(argv[2], "%d", &B_bez_ord);
284
285
286 init( argc, argv, new IntersectImplicit(A_bez_ord, B_bez_ord));
287 return 0;
288 }
289
290
291 /*
292 Local Variables:
293 mode:c++
294 c-file-style:"stroustrup"
295 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
296 indent-tabs-mode:nil
297 fill-column:99
298 End:
299 */
300 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
301