1
2 /******************************************************************************
3 * MODULE : grid.cpp
4 * DESCRIPTION: grids for the graphics
5 * COPYRIGHT : (C) 2003 Henri Lesourd
6 *******************************************************************************
7 * This software falls under the GNU general public license version 3 or later.
8 * It comes WITHOUT ANY WARRANTY WHATSOEVER. For details, see the file LICENSE
9 * in the root directory or <http://www.gnu.org/licenses/gpl-3.0.html>.
10 ******************************************************************************/
11
12 #include "grid.hpp"
13 #include "math_util.hpp"
14 #include "curve.hpp"
15
16 /******************************************************************************
17 * General routines
18 ******************************************************************************/
19
20 void
set_aspect(tree aspect)21 grid_rep::set_aspect (tree aspect) {
22 subd= array<SI> (2);
23 subd[0]= 0;
24 subd[1]= 1;
25 col= array<string> (2);
26 col[0]= "#808080";
27 col[1]= "#c0c0c0";
28 if (is_tuple (aspect)) {
29 int i;
30 bool b= false;
31 subd= array<SI> (N(aspect));
32 col= array<string> (N(aspect));
33 for (i=0; i<N(aspect); i++) {
34 if (is_tuple (aspect[i], "axes", 1)) {
35 subd[i]= 0;
36 b= true;
37 }
38 else {
39 subd[i]= as_int (aspect[i][0]);
40 }
41 col[i]= as_string (aspect[i][1]);
42 }
43 if (!b) {
44 array<SI> subd0 (1);
45 array<string> col0 (1);
46 subd0[0]= 0;
47 col0[0]= "#e0e0ff";
48 subd= subd0 << subd;
49 col= col0 << col;
50 }
51 do {
52 b= true;
53 for (i=1; i<N(subd); i++)
54 if (subd[i-1]>subd[i]) {
55 SI j;
56 string c;
57 j= subd[i-1];subd[i-1]= subd[i];subd[i]= j;
58 c= col[i-1];col[i-1]= col[i];col[i]= c;
59 b= false;
60 }
61 }
62 while (!b);
63 }
64 }
65
66 array<grid_curve>
get_curves_around(point p,double delta,frame f)67 grid_rep::get_curves_around (point p, double delta, frame f) {
68 point p2= f (p);
69 point pmin, pmax;
70 pmin= f[point (p2[0]-delta, p2[1]-delta)];
71 pmax= f[point (p2[0]+delta, p2[1]+delta)];
72 return get_curves (pmin, pmax, true);
73 }
74
75 point
find_point_around(point p,double delta,frame f)76 grid_rep::find_point_around (point p, double delta, frame f) {
77 point p2= f (p);
78 point pmin, pmax;
79 pmin= f[point (p2[0]-delta, p2[1]-delta)];
80 pmax= f[point (p2[0]+delta, p2[1]+delta)];
81 return find_closest_point (p, pmin, pmax);
82 }
83
84 /******************************************************************************
85 * The empty grid
86 ******************************************************************************/
87
88 struct empty_grid_rep: public grid_rep {
empty_grid_repempty_grid_rep89 empty_grid_rep ():
90 grid_rep () {}
get_curvesempty_grid_rep91 array<grid_curve> get_curves (point lim1, point lim2, bool b) {
92 (void) lim1; (void) lim2; (void) b;
93 array<grid_curve> res;
94 return res; }
operator treeempty_grid_rep95 operator tree () { return "empty_grid"; }
96 point find_closest_point (point p, point pmin, point pmax);
97 };
98
99 point
find_closest_point(point p,point pmin,point pmax)100 empty_grid_rep::find_closest_point (point p, point pmin, point pmax) {
101 (void) pmin; (void) pmax;
102 /*double x= floor (10.0*p[0] + 0.5);
103 double y= floor (10.0*p[1] + 0.5);
104 return point (x / 10.0, y / 10.0);*/
105 return p;
106 }
107
108 grid
empty_grid()109 empty_grid () {
110 return tm_new<empty_grid_rep> ();
111 }
112
113 /******************************************************************************
114 * Cartesian grids
115 ******************************************************************************/
116
117 struct cartesian_rep: public grid_rep {
118 double step; // Unit length
cartesian_repcartesian_rep119 cartesian_rep (array<SI> subd, array<string> col, point o, double st):
120 grid_rep (subd, col, o), step (st) {}
operator treecartesian_rep121 operator tree () { return "cartesian"; }
122 array<grid_curve> get_curves (point lim1, point lim2, bool b);
123 array<grid_curve> get_curves_around (point p, double delta, frame f);
124 point find_closest_point (point p, point pmin, point pmax);
125 };
126
127 static grid_curve
create_line(double x1,double y1,double x2,double y2,string col)128 create_line (double x1, double y1, double x2, double y2, string col) {
129 array<point> a(2);
130 a[0]= point (x1, y1);
131 a[1]= point (x2, y2);
132 array<path> cip(2);
133 grid_curve res= grid_curve (col, poly_segment (a, cip));
134 return res;
135 }
136
137 array<grid_curve>
get_curves(point lim1,point lim2,bool b)138 cartesian_rep::get_curves (point lim1, point lim2, bool b) {
139 array<grid_curve> res;
140 if (N(subd)<1) return res;
141 double x1= min (lim1[0], lim2[0]);
142 double y1= min (lim1[1], lim2[1]);
143 double x2= max (lim1[0], lim2[0]);
144 double y2= max (lim1[1], lim2[1]);
145 double xo= center[0];
146 double yo= center[1];
147 int i;
148 for (i= N(subd)-1; i>=1; i--) {
149 SI nsub;
150 nsub= subd[i];
151 if (nsub!=0) {
152 double x, y, s;
153 s= max (step/nsub, 1.0e-6);
154 for (x= xo; x<=x2; x+=s)
155 if (x>=x1) {
156 res << create_line (x, y1, x, y2, col[i]);
157 if (b) break;
158 }
159 for (x= xo-s; x>=x1; x-=s)
160 if (x<=x2) {
161 res << create_line (x, y1, x, y2, col[i]);
162 if (b) break;
163 }
164 for (y= yo; y<=y2; y+=s)
165 if (y>=y1) {
166 res << create_line (x1, y, x2, y, col[i]);
167 if (b) break;
168 }
169 for (y= yo-s; y>=y1; y-=s)
170 if (y<=y2) {
171 res << create_line (x1, y, x2, y, col[i]);
172 if (b) break;
173 }
174 }
175 }
176 if (!b && yo>=y1 && yo<=y2)
177 res << create_line (x1, yo, x2, yo, col[0]);
178 if (!b && xo>=x1 && xo<=x2)
179 res << create_line (xo, y1, xo, y2, col[0]);
180 return res;
181 }
182
183 array<grid_curve>
get_curves_around(point p,double delta,frame f)184 cartesian_rep::get_curves_around (point p, double delta, frame f) {
185 double nsub= 0;
186 string c;
187 for (int i=0; i<N(subd); i++)
188 if (subd[i] != 0) { nsub= subd[i]; c= col[i]; }
189 if (nsub == 0) return array<grid_curve> (0);
190
191 point p2 = f (p);
192 point lim1= f[point (p2[0]-delta, p2[1]-delta)];
193 point lim2= f[point (p2[0]+delta, p2[1]+delta)];
194
195 double x1= min (lim1[0], lim2[0]);
196 double y1= min (lim1[1], lim2[1]);
197 double x2= max (lim1[0], lim2[0]);
198 double y2= max (lim1[1], lim2[1]);
199 double xo= center[0];
200 double yo= center[1];
201 double s = step / nsub;
202
203 double X1= xo + floor ((p[0] - xo) / s) * s;
204 double Y1= yo + floor ((p[1] - yo) / s) * s;
205 double X2= xo + ceil ((p[0] - xo) / s) * s;
206 double Y2= yo + ceil ((p[1] - yo) / s) * s;
207 x1= min (x1, X1 - s / 10);
208 y1= min (y1, Y1 - s / 10);
209 x2= max (x2, X2 + s / 10);
210 y2= max (y2, Y2 + s / 10);
211
212 array<grid_curve> res;
213 res << create_line (X1, y1, X1, y2, c);
214 res << create_line (x1, Y1, x2, Y1, c);
215 if (X2 > X1) res << create_line (X2, y1, X2, y2, c);
216 if (Y2 > Y1) res << create_line (x1, Y2, x2, Y2, c);
217 return res;
218 }
219
220 point
find_closest_point(point p,point pmin,point pmax)221 cartesian_rep::find_closest_point (point p, point pmin, point pmax) {
222 double x, y, oldx=0, oldy= 0, ssubd;
223 point res= p;
224 p= p-center;
225 int i;
226 for (i=1; i<N(subd); i++) {
227 ssubd= ((double) subd[i]) / step;
228 if (ssubd==0) continue;
229 x= nearest (p[0]*ssubd);
230 y= nearest (p[1]*ssubd);
231 res= center + point (x/ssubd, y/ssubd);
232 if (i!=1) {
233 if (inside_rectangle (point (oldx, res[1]), pmin, pmax))
234 return point (oldx, res[1]);
235 if (inside_rectangle (point (res[0], oldy), pmin, pmax))
236 return point (res[0], oldy);
237 }
238 oldx= res[0];
239 oldy= res[1];
240 if (inside_rectangle (res, pmin, pmax))
241 return res;
242 }
243 return res;
244 }
245
246 grid
cartesian(array<SI> subd,array<string> col,point o,double step)247 cartesian (array<SI> subd, array<string> col, point o, double step) {
248 return tm_new<cartesian_rep> (subd, col, o, step);
249 }
250
251 /******************************************************************************
252 * Polar grids
253 ******************************************************************************/
254
255 struct polar_rep: public grid_rep {
256 double step; // Radial unit length
257 SI astep; // # angles
polar_reppolar_rep258 polar_rep (array<SI> subd, array<string> col, point o, double st, SI ast):
259 grid_rep (subd, col, o), step (st), astep (ast) {}
operator treepolar_rep260 operator tree () { return "polar"; }
261 array<grid_curve> get_curves (point lim1, point lim2, bool b);
262 point find_closest_point (point p, point pmin, point pmax);
263 };
264
265 static grid_curve
create_arc(double x1,double y1,double x2,double y2,double x3,double y3,string col)266 create_arc (
267 double x1, double y1, double x2, double y2, double x3, double y3, string col)
268 {
269 array<point> a(3);
270 a[0]= point (x1, y1);
271 a[1]= point (x2, y2);
272 a[2]= point (x3, y3);
273 array<path> cip(3);
274 grid_curve res= grid_curve (col, arc (a, cip, true));
275 return res;
276 }
277
278 array<grid_curve>
get_curves(point lim1,point lim2,bool b)279 polar_rep::get_curves (point lim1, point lim2, bool b) {
280 (void) b;
281 array<grid_curve> res;
282 if (N(subd)<1) return res;
283 double x1= min (lim1[0], lim2[0]);
284 double y1= min (lim1[1], lim2[1]);
285 double x2= max (lim1[0], lim2[0]);
286 double y2= max (lim1[1], lim2[1]);
287 double xo= center[0];
288 double yo= center[1];
289 point P1, P2;
290 if (x1<=0 && y1<=0 && x2>=0 && y2>=0) {
291 P1= point (0, 0);
292 P2= point (x2, y2) - point (x1, y1);
293 }
294 else {
295 double ox= (x1 + x2) / 2;
296 double oy= (y1 + y2) / 2;
297 if (oy>=0)
298 P1= point (0, y1>=0 ? y1 : 0);
299 else
300 P1= point (0, y2<=0 ? y2 : 0);
301
302 if (ox>=0 && oy>=0)
303 P2= point (x2, y2);
304 else
305 if (ox<=0 && oy>=0)
306 P2= point (x1, y2);
307 else
308 if (ox<=0 && oy<=0)
309 P2= point (x1, y1);
310 else
311 if (ox>=0 && oy<=0)
312 P2= point (x2, y1);
313 }
314 double r, R1= norm (P1), R2= norm (P2);
315 int i;
316 for (i= N(subd)-1; i>=1; i--) {
317 SI nsub;
318 nsub= subd[i];
319 if (nsub!=0) {
320 SI j;
321 double s= max (step/nsub, 1.0e-6);
322 for (r=0; r<=R2; r+=s)
323 if (r>=R1)
324 res << create_arc (xo+r, yo, xo, yo+r, xo-r, yo, col[i]);
325 for (j=0; j<astep*nsub; j++)
326 res << create_line (xo, yo, xo+R2*cos((2*tm_PI*j)/(astep*nsub)),
327 yo+R2*sin((2*tm_PI*j)/(astep*nsub)),
328 col[i]);
329 }
330 }
331 res << create_line (x1, yo, x2, yo, col[0]);
332 res << create_line (xo, y1, xo, y2, col[0]);
333 return res;
334 }
335
336 point
find_closest_point(point p,point pmin,point pmax)337 polar_rep::find_closest_point (point p, point pmin, point pmax) {
338 double n, a, oldn= 0, olda= 0, ssubd;
339 point res= p;
340 p= p-center;
341 int i;
342 for (i=1; i<N(subd); i++) {
343 ssubd= (double)subd[i];
344 if (ssubd==0) continue;
345 n= nearest (norm(p)*(ssubd/step));
346 if (fnull (norm (p), 1.0e-6))
347 a= 0.0;
348 else
349 a= nearest ((arg(p)/(2*tm_PI))*astep*ssubd);
350 n= n*(step/ssubd);
351 a= a/(astep*ssubd);
352 if (i!=1) {
353 res= center + oldn * point (cos(2*tm_PI*a), sin(2*tm_PI*a));
354 if (inside_rectangle (res, pmin, pmax))
355 return res;
356 res= center + n * point (cos(2*tm_PI*olda), sin(2*tm_PI*olda));
357 if (inside_rectangle (res, pmin, pmax))
358 return res;
359 }
360 oldn= n;
361 olda= a;
362 res= center + n * point (cos(2*tm_PI*a), sin(2*tm_PI*a));
363 if (inside_rectangle (res, pmin, pmax))
364 return res;
365 }
366 return res;
367 }
368
369 grid
polar(array<SI> subd,array<string> col,point o,double step,SI astep)370 polar (array<SI> subd, array<string> col, point o, double step, SI astep) {
371 return tm_new<polar_rep> (subd, col, o, step, astep);
372 }
373
374 /******************************************************************************
375 * Logarithmic grids
376 ******************************************************************************/
377
378 struct logarithmic_rep: public grid_rep {
379 double step; // Unit length
380 SI base;
logarithmic_replogarithmic_rep381 logarithmic_rep (array<SI> subd, array<string> col, point o, double st, SI b):
382 grid_rep (subd, col, o), step (max (st, 1.0e-6)), base (b) {}
operator treelogarithmic_rep383 operator tree () { return "logarithmic"; }
384 array<grid_curve> get_curves (point lim1, point lim2, bool b);
385 point find_closest_point (point p, point pmin, point pmax);
386 };
387
388 array<grid_curve>
get_curves(point lim1,point lim2,bool b)389 logarithmic_rep::get_curves (point lim1, point lim2, bool b) {
390 (void) b;
391 array<grid_curve> res;
392 if (N(subd)<1) return res;
393 double x1= min (lim1[0], lim2[0]);
394 double y1= min (lim1[1], lim2[1]);
395 double x2= max (lim1[0], lim2[0]);
396 double y2= max (lim1[1], lim2[1]);
397 double xo= center[0];
398 double yo= center[1];
399 int i;
400 double x, y;
401 if (N(col)>=3) {
402 for (i=2; i<base; i++) {
403 double dx, dy;
404 dx= dy= step*log((double)i)/log((double)base);
405 for (x=xo; x<=x2; x+=step)
406 res << create_line (x+dx, y1, x+dx, y2, col[2]);
407 for (x=xo-step; x>=x1-step; x-=step)
408 res << create_line (x+dx, y1, x+dx, y2, col[2]);
409 for (y=yo; y<=y2; y+=step)
410 res << create_line (x1, y+dy, x2, y+dy, col[2]);
411 for (y=yo-step; y>=y1-step; y-=step)
412 res << create_line (x1, y+dy, x2, y+dy, col[2]);
413 }
414 }
415 if (N(col)>=2) {
416 for (x=xo; x<=x2; x+=step)
417 res << create_line (x, y1, x, y2, col[1]);
418 for (x=xo; x>=x1; x-=step)
419 res << create_line (x, y1, x, y2, col[1]);
420 for (y=yo; y<=y2; y+=step)
421 res << create_line (x1, y, x2, y, col[1]);
422 for (y=yo; y>=y1; y-=step)
423 res << create_line (x1, y, x2, y, col[1]);
424 }
425 res << create_line (x1, yo, x2, yo, col[0]);
426 res << create_line (xo, y1, xo, y2, col[0]);
427 return res;
428 }
429
430 point
find_closest_point(point p,point pmin,point pmax)431 logarithmic_rep::find_closest_point (point p, point pmin, point pmax) {
432 double x, y, ssubd= ((double) subd[1]) / step;
433 point res= p;
434 if (ssubd!=0) {
435 p= p-center;
436 x= nearest (p[0]*ssubd);
437 y= nearest (p[1]*ssubd);
438 res= center + point (x/ssubd, y/ssubd);
439 if (inside_rectangle (res, pmin, pmax))
440 return res;
441 p= center+p;
442 }
443 double xo, yo;
444 xo= center[0];
445 yo= center[1];
446 double x1, y1, x2, y2;
447 x1= xo-step;
448 y1= yo-step;
449 x2= xo+step;
450 y2= yo+step;
451 p= p-center;
452 double x0, y0;
453 x0= (SI)(p[0]/step);
454 y0= (SI)(p[1]/step);
455 x0*= step;
456 y0*= step;
457 p= p-point(x0,y0);
458 p= center+p;
459 double xm, ym;
460 xm= ym= tm_infinity/2;
461 int i;
462 for (i=1; i<base; i++) {
463 double dx, dy;
464 dx= dy= step*log((double)i)/log((double)base);
465 for (x=xo; x<=x2; x+=step)
466 if (norm(x+dx-p[0])<norm(xm-p[0])) xm= x+dx;
467 for (x=xo-step; x>=x1-step; x-=step)
468 if (norm(x+dx-p[0])<norm(xm-p[0])) xm= x+dx;
469 for (y=yo; y<=y2; y+=step)
470 if (norm(y+dy-p[1])<norm(ym-p[1])) ym= y+dy;
471 for (y=yo-step; y>=y1-step; y-=step)
472 if (norm(y+dy-p[1])<norm(ym-p[1])) ym= y+dy;
473 }
474 p= point (x0+xm, y0+ym);
475 if (ssubd!=0) {
476 if (inside_rectangle (point (res[0], p[1]), pmin, pmax))
477 return point (res[0], p[1]);
478 if (inside_rectangle (point (p[0], res[1]), pmin, pmax))
479 return point (p[0], res[1]);
480 }
481 return p;
482 }
483
484 grid
logarithmic(array<SI> subd,array<string> col,point o,double step,SI base)485 logarithmic (array<SI> subd, array<string> col, point o, double step, SI base) {
486 return tm_new<logarithmic_rep> (subd, col, o, step, base);
487 }
488
489 /******************************************************************************
490 * User interface
491 ******************************************************************************/
492
493 grid
as_grid(tree t)494 as_grid (tree t) {
495 array<SI> subd (0, 1);
496 array<string> col ("black", "black");
497 grid gr= empty_grid ();
498 double step= 1.0;
499 point center= point (0.0, 0.0);
500 if (is_tuple (t, "empty"))
501 gr= empty_grid ();
502 else
503 if (is_tuple (t, "cartesian")) {
504 if (is_tuple (t, "cartesian", 0)) ;
505 else
506 if (is_tuple (t, "cartesian", 1))
507 step= as_double (t[1]);
508 else
509 if (is_tuple (t, "cartesian", 2)) {
510 center= as_point (t[1]);
511 step= as_double (t[2]);
512 }
513 gr= cartesian (subd, col, center, step);
514 }
515 else
516 if (is_tuple (t, "polar")) {
517 SI astep= 8;
518 if (is_tuple (t, "polar", 0)) ;
519 else
520 if (is_tuple (t, "polar", 1))
521 step= as_double (t[1]);
522 else
523 if (is_tuple (t, "polar", 2)) {
524 step= as_double (t[1]);
525 astep= as_int (t[2]);
526 }
527 else
528 if (is_tuple (t, "polar", 3)) {
529 center= as_point (t[1]);
530 step= as_double (t[2]);
531 astep= as_int (t[3]);
532 }
533 gr=polar (subd, col, center, step, astep);
534 }
535 else
536 if (is_tuple (t, "logarithmic")) {
537 SI base= 10;
538 if (is_tuple (t, "logarithmic", 0)) ;
539 else
540 if (is_tuple (t, "logarithmic", 1))
541 step= as_double (t[1]);
542 else
543 if (is_tuple (t, "logarithmic", 2)) {
544 step= as_double (t[1]);
545 base= as_int (t[2]);
546 }
547 else
548 if (is_tuple (t, "logarithmic", 3)) {
549 center= as_point (t[1]);
550 step= as_double (t[2]);
551 base= as_int (t[3]);
552 }
553 gr= logarithmic (subd, col, center, step, base);
554 }
555 return gr;
556 }
557
558 tree
as_tree(grid g)559 as_tree (grid g) {
560 return (tree) g;
561 }
562