1 #include <iostream>
2 #include <2geom/path.h>
3 #include <2geom/svg-path-parser.h>
4 #include <2geom/path-intersection.h>
5 #include <2geom/basic-intersection.h>
6 #include <2geom/pathvector.h>
7 #include <2geom/exception.h>
8 #include <2geom/sbasis-geometric.h>
9 #include <2geom/path-intersection.h>
10 #include <2geom/nearest-time.h>
11 #include <2geom/line.h>
12 #include <2geom/bezier-to-sbasis.h>
13 #include <2geom/sbasis-to-bezier.h>
14
15 #include <cstdlib>
16 #include <map>
17 #include <vector>
18 #include <algorithm>
19 #include <optional>
20
21 #include <toys/path-cairo.h>
22 #include <toys/toy-framework-2.h>
23 #include <2geom/bezier-to-sbasis.h>
24 #include <2geom/ord.h>
25
26 #include <2geom/conicsec.h>
27
xAx_to_RatQuads(Geom::xAx const &,Geom::Rect const &)28 std::vector<Geom::RatQuad> xAx_to_RatQuads(Geom::xAx const &/*C*/, Geom::Rect const &/*bnd*/) {
29 // find points on boundary
30 // if there are exactly 0 points return
31 // if there are exactly 2 points fit ratquad and return
32 // if there are an odd number, split bnd on the point with the smallest dot(unit_vector(grad), rect_edge)
33 // sort into clockwise order ABCD
34 // compute corresponding tangents
35 // test boundary points against the line through A
36 // if all on one side
37 //
38 // if A,X and Y,Z
39 // ratquad from A,X and Y,Z
40 return std::vector<Geom::RatQuad>();
41 }
42
43
44
45 using namespace Geom;
46 using namespace std;
47
48
49 // File: convert.h
50 #include <sstream>
51 #include <stdexcept>
52
53 class BadConversion : public std::runtime_error {
54 public:
BadConversion(const std::string & s)55 BadConversion(const std::string& s)
56 : std::runtime_error(s)
57 { }
58 };
59
60 template <typename T>
stringify(T x)61 inline std::string stringify(T x)
62 {
63 std::ostringstream o;
64 if (!(o << x))
65 throw BadConversion("stringify(T)");
66 return o.str();
67 }
68
69 namespace Geom{
70 xAx degen;
71 };
72
draw_hull(cairo_t * cr,RatQuad rq)73 void draw_hull(cairo_t*cr, RatQuad rq) {
74 cairo_move_to(cr, rq.P[0]);
75 cairo_line_to(cr, rq.P[1]);
76 cairo_line_to(cr, rq.P[2]);
77 cairo_stroke(cr);
78 }
79
80
81
draw(cairo_t * cr,xAx C,Rect bnd)82 void draw(cairo_t* cr, xAx C, Rect bnd) {
83 if(bnd[1].extent() < 5) return;
84 vector<double> prev_rts;
85 double py = bnd[Y].min();
86 for(int i = 0; i < 100; i++) {
87 double t = i/100.;
88 double y = bnd[Y].valueAt(t);
89 vector<double> rts = C.roots(Point(1, 0), Point(0, y));
90 int top = 0;
91 for(unsigned j = 0; j < rts.size(); j++) {
92 if(bnd[0].contains(rts[j])) {
93 rts[top] = rts[j];
94 top++;
95 }
96 }
97 rts.erase(rts.begin()+top, rts.end());
98
99 if(rts.size() == prev_rts.size()) {
100 for(unsigned j = 0; j < rts.size(); j++) {
101 cairo_move_to(cr, prev_rts[j], py);
102 cairo_line_to(cr, rts[j], y);
103 cairo_stroke(cr);
104 }
105 /* } else if(prev_rts.size() == 1) {
106 for(unsigned j = 0; j < rts.size(); j++) {
107 cairo_move_to(cr, prev_rts[0], py);
108 cairo_line_to(cr, rts[j], y);
109 cairo_stroke(cr);
110 }
111 } else if(rts.size() == 1) {
112 for(unsigned j = 0; j < prev_rts.size(); j++) {
113 cairo_move_to(cr, prev_rts[j], py);
114 cairo_line_to(cr, rts[0], y);
115 cairo_stroke(cr);
116 }*/
117 } else {
118 draw(cr, C, Rect(bnd[0], Interval(py, y)));
119 /*for(unsigned j = 0; j < rts.size(); j++) {
120 cairo_move_to(cr, rts[j], y);
121 cairo_rel_line_to(cr, 1,1);
122 }*/
123 }
124 prev_rts = rts;
125 py = y;
126 }
127 }
128
129 template <typename T>
det(T a,T b,T c,T d)130 static T det(T a, T b, T c, T d) {
131 return a*d - b*c;
132 }
133
134 template <typename T>
det(T M[2][2])135 static T det(T M[2][2]) {
136 return M[0][0]*M[1][1] - M[1][0]*M[0][1];
137 }
138
139 template <typename T>
det3(T M[3][3])140 static T det3(T M[3][3]) {
141 return ( M[0][0] * det(M[1][1], M[1][2],
142 M[2][1], M[2][2])
143 -M[1][0] * det(M[0][1], M[0][2],
144 M[2][1], M[2][2])
145 +M[2][0] * det(M[0][1], M[0][2],
146 M[1][1], M[1][2]));
147 }
148
149 class Conic6: public Toy {
150 PointSetHandle C1H, C2H;
151 std::vector<Slider> sliders;
152 Point mouse_sampler;
153
mouse_moved(GdkEventMotion * e)154 void mouse_moved(GdkEventMotion* e) override {
155 mouse_sampler = Point(e->x, e->y);
156 Toy::mouse_moved(e);
157 }
158
draw(cairo_t * cr,std::ostringstream * notify,int width,int height,bool save,std::ostringstream * timer_stream)159 void draw(cairo_t *cr, std::ostringstream *notify, int width, int height, bool save, std::ostringstream *timer_stream) override {
160 cairo_set_source_rgba (cr, 0., 0., 0, 1);
161 cairo_set_line_width (cr, 1);
162 Rect screen_rect(Interval(10, width-10), Interval(10, height-10));
163
164 Geom::xAx C1 = xAx::fromPoints(C1H.pts);
165 ::draw(cr, C1, screen_rect);
166 *notify << C1;
167
168 Geom::xAx C2 = xAx::fromPoints(C2H.pts);
169 ::draw(cr, C2, screen_rect);
170 *notify << C2;
171
172
173 SBasis T(Linear(-1,1));
174 SBasis S(Linear(1,1));
175 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},
176 {(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},
177 {(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]}};
178
179 SBasis D = det3(C);
180 std::vector<double> rts = Geom::roots(D);
181 if(rts.empty()) {
182 T = Linear(1,1);
183 S = Linear(-1,1);
184 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},
185 {(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},
186 {(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]}};
187
188 D = det3(C);
189 rts = Geom::roots(D);
190 }
191 // at this point we have a T and S and perhaps some roots that represent our degenerate conic
192 // Let's just pick one randomly (can we do better?)
193 //for(unsigned i = 0; i < rts.size(); i++) {
194 if(!rts.empty()) {
195 cairo_save(cr);
196
197 unsigned i = 0;
198 double t = T.valueAt(rts[i]);
199 double s = S.valueAt(rts[i]);
200 *notify << t << "; " << s << std::endl;
201 /*double C0[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},
202 {(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},
203 {(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]}};*/
204 xAx xC0 = C1*t + C2*s;
205 //::draw(cr, xC0, screen_rect); // degen
206
207 std::optional<Point> oB0 = xC0.bottom();
208
209 Point B0 = *oB0;
210 //*notify << B0 << " = " << C1.gradient(B0);
211 draw_circ(cr, B0);
212
213 Point n0, n1;
214 // Are these just the eigenvectors of A11?
215 if(fabs(xC0.c[0]) > fabs(xC0.c[2])) {
216 double b = 0.5*xC0.c[1]/xC0.c[0];
217 double c = xC0.c[2]/xC0.c[0];
218 double d = std::sqrt(b*b-c);
219 n0 = Point(1, b+d);
220 n1 = Point(1, b-d);
221 } else {
222
223 double b = 0.5*xC0.c[1]/xC0.c[2];
224 double c = xC0.c[0]/xC0.c[2];
225 double d = std::sqrt(b*b-c);
226 n0 = Point(b+d, 1);
227 n1 = Point(b-d, 1);
228 }
229 cairo_set_source_rgb(cr, 0.7, 0.7, 0.7);
230
231 Line L0 = Line::from_origin_and_vector(B0, rot90(n0));
232 draw_line(cr, L0, screen_rect);
233 Line L1 = Line::from_origin_and_vector(B0, rot90(n1));
234 draw_line(cr, L1, screen_rect);
235
236 cairo_set_source_rgb(cr, 1, 0., 0.);
237 rts = C1.roots(L0);
238 for(double rt : rts) {
239 Point P = L0.pointAt(rt);
240 draw_cross(cr, P);
241 *notify << C1.valueAt(P) << "; " << C2.valueAt(P) << "\n";
242 }
243 rts = C1.roots(L1);
244 for(double rt : rts) {
245 Point P = L1.pointAt(rt);
246 draw_cross(cr, P);
247 *notify << C1.valueAt(P) << "; "<< C2.valueAt(P) << "\n";
248 }
249 cairo_stroke(cr);
250 cairo_restore(cr);
251 }
252
253 ::draw(cr, C1*sliders[0].value() + C2*sliders[1].value(), screen_rect);
254
255 std::vector<Point> res = intersect(C1, C2);
256 for(auto & re : res) {
257 draw_circ(cr, re);
258 }
259
260 cairo_stroke(cr);
261
262 //*notify << "w = " << w << "; lambda = " << rq.lambda() << "\n";
263 Toy::draw(cr, notify, width, height, save, timer_stream);
264 }
265
266 public:
Conic6()267 Conic6() {
268 for(int j = 0; j < 5; j++){
269 C1H.push_back(uniform()*400, 100+ uniform()*300);
270 C2H.push_back(uniform()*400, 100+ uniform()*300);
271 }
272 handles.push_back(&C1H);
273 handles.push_back(&C2H);
274 sliders.emplace_back(-1.0, 1.0, 0, 0.0, "a");
275 sliders.emplace_back(-1.0, 1.0, 0, 0.0, "b");
276 sliders.emplace_back(0.0, 5.0, 0, 0.0, "c");
277 handles.push_back(&(sliders[0]));
278 handles.push_back(&(sliders[1]));
279 handles.push_back(&(sliders[2]));
280 sliders[0].geometry(Point(50, 20), 250);
281 sliders[1].geometry(Point(50, 50), 250);
282 sliders[2].geometry(Point(50, 80), 250);
283 }
284
first_time(int,char **)285 void first_time(int /*argc*/, char**/* argv*/) override {
286
287 }
288 };
289
main(int argc,char ** argv)290 int main(int argc, char **argv) {
291 init(argc, argv, new Conic6());
292 return 0;
293 }
294
295 /*
296 Local Variables:
297 mode:c++
298 c-file-style:"stroustrup"
299 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
300 indent-tabs-mode:nil
301 fill-column:99
302 End:
303 */
304 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4:fileencoding=utf-8:textwidth=99 :
305