1 /*
2 * $Id: pl3d.i,v 1.1 2005-09-18 22:06:03 dhmunro Exp $
3 * Viewing transforms and other aids for 3D plotting.
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 /* General overview:
12
13 (1) Viewing transform machinery. Arguably the simplest model
14 is the CAD/CAM notion that the object you see is oriented
15 as you see it in the current picture. You can then move
16 it left, right, up, down, or toward or away from you,
17 or you can rotate it about any of the three axes (horizontal,
18 vertical, or out of the screen). The xyz coordinates of the
19 object remains unchanged throughout all of this, but this
20 object coordinate system changes relative to the fixed
21 xyz of the viewer, in which x is always to the right, y is
22 up, and z is directed out of the screen. Initially, the
23 two coordinate systems coincide.
24
25 rot3, xangle,yangle,zangle
26 Rotate the object about viewer's x-axis by xangle, then
27 about viewer's y-axis by yangle, then about viewer's
28 z-axis by zangle
29 mov3, xchange,ychange,zchange
30 Move the object by the specified amounts.
31
32 setz3, zcamera
33 The "camera" is located at (0,0,zcamera) in the viewer's
34 coordinate system, looking in the minus-z direction.
35 Initially, zcamera is very large, and the magnification
36 factor is correspondingly large, giving an isometric view.
37 Decreasing zcamera makes the perspective more extreme.
38 If parts of the object are behind the camera, strange things
39 may happen.
40
41 undo3
42 undo3, n
43 Undo the last N (default 1) viewpoint commands (rot3, mov3,
44 or setz3). Up to 100 viewpoint changes are remembered.
45 viewpoint= save3()
46 ...
47 restore3, viewpoint
48 The current viewpoint transformation can be saved and later
49 restored.
50
51 gnomon, on_off
52 Toggle the gnomon (a simple display showing the orientation
53 of the xyz axes of the object).
54 */
55
56 /* ------------------------------------------------------------------------ */
57
rot3(xa,ya,za)58 func rot3(xa,ya,za)
59 /* DOCUMENT rot3, xa,ya,za
60 rotate the current 3D plot by XA about viewer's x-axis,
61 YA about viewer's y-axis, and ZA about viewer's z-axis.
62 SEE ALSO: orient3, mov3, aim3, setz3, undo3, save3, restore3, light3
63 */
64 {
65 if (is_void(xa)) xa= 0.;
66 if (is_void(ya)) ya= 0.;
67 if (is_void(za)) za= 0.;
68 x= [1.,0.,0.];
69 y= [0.,1.,0.];
70 z= [0.,0.,1.];
71 _rot3, za, x, y;
72 _rot3, ya, z, x;
73 _rot3, xa, y, z;
74 _setrot3, [x,y,z](,+)*_getrot3()(+,);
75 }
76
77 func _rot3(a, &x, &y)
78 {
79 ca= cos(a);
80 sa= sin(a);
81 a= x;
82 x= ca*a + sa*y;
83 y= -sa*a + ca*y;
84 }
85
mov3(xa,ya,za)86 func mov3(xa,ya,za)
87 /* DOCUMENT mov3, xa,ya,za
88 move the current 3D plot by XA along viewer's x-axis,
89 YA along viewer's y-axis, and ZA along viewer's z-axis.
90 SEE ALSO: rot3, orient3, setz3, undo3, save3, restore3, light3
91 */
92 {
93 if (is_void(xa)) xa= 0.;
94 if (is_void(ya)) ya= 0.;
95 if (is_void(za)) za= 0.;
96 _setorg3, _getorg3() - _getrot3()(+,)*[xa,ya,za](+);
97 }
98
aim3(xa,ya,za)99 func aim3(xa,ya,za)
100 /* DOCUMENT aim3, xa,ya,za
101 move the current 3D plot to put the point (XA,YA,ZA) in object
102 coordinates at the point (0,0,0) -- the aim point -- in the
103 viewer's coordinates. If any of XA, YA, or ZA is nil, it defaults
104 to zero.
105 SEE ALSO: mov3, rot3, orient3, setz3, undo3, save3, restore3, light3
106 */
107 {
108 if (is_void(xa)) xa= 0.;
109 if (is_void(ya)) ya= 0.;
110 if (is_void(za)) za= 0.;
111 _setorg3, [xa,ya,za];
112 }
113
setz3(zc)114 func setz3(zc)
115 /* DOCUMENT setz3, zc
116 Set the camera position to z=ZC (x=y=0) in the viewer's coordinate
117 system. If ZC is nil, set the camera to infinity (default).
118 SEE ALSO: rot3, orient3, undo3, save3, restore3, light3
119 */
120 {
121 if (!is_void(zc)) {
122 if (dimsof(zc)(1)) error, "camera position must be scalar";
123 zc= double(zc);
124 }
125 _setzc3, zc;
126 }
127
orient3(phi,theta)128 func orient3(phi, theta)
129 /* DOCUMENT orient3, phi, theta
130 or orient3, phi
131 or orient3, , theta
132 or orient3
133 Set the "orientation" of the object to (PHI,THETA). "Orientations"
134 are a subset of the possible rotation matrices in which the z-axis
135 of the object appears vertical on the screen (that is, the object
136 z-axis projects onto the viewer y-axis). The THETA angle is the
137 angle from the viewer y-axis to the object z-axis, positive if
138 the object z-axis is tilted toward you (toward viewer +z). PHI is
139 zero when the object x-axis coincides with the viewer x-axis. If
140 neither PHI nor THETA is specified, PHI defaults to -pi/4 and
141 THETA defaults to pi/6. If only one of PHI or THETA is specified,
142 the other remains unchanged, unless the current THETA is near pi/2,
143 in which case THETA returns to pi/6, or unless the current
144 orientation does not have a vertical z-axis, in which case the
145 unspecified value returns to its default.
146
147 Unlike rot3, orient3 is not a cumulative operation.
148
149 SEE ALSO: rot3, mov3, aim3, save3, restore3, light3, limit3
150 */
151 {
152 if (is_void(theta) && is_void(phi)) {
153 theta= _orient3_theta;
154 phi= _orient3_phi;
155 } else if (is_void(theta) || is_void(phi)) {
156 z= _getrot3()(,+)*[0.,0.,1.](+);
157 if (abs(z(1))>1.e-6) {
158 /* object z-axis not aligned with viewer y-axis */
159 if (is_void(theta)) theta= _orient3_theta;
160 else phi= _orient3_phi;
161 } else if (is_void(theta)) {
162 if (abs(z(2))<1.e-6) theta= _orient3_theta;
163 else theta= atan(z(3),z(2));
164 } else /*if (is_void(phi))*/ {
165 y= [0.,z(3),-z(2)]; /* in object xy-plane */
166 x= _getrot3()(,+)*[1.,0.,0.](+);
167 phi= atan(y(+)*x(+), x(1));
168 }
169 }
170 x= [1.,0.,0.];
171 y= [0.,1.,0.];
172 z= [0.,0.,1.];
173 _rot3, theta, y, z;
174 _rot3, phi, z, x;
175 _setrot3, [x,-z,y];
176 }
177
178 /* unless user has supplied alternative defaults, set orient3 defaults */
179 if (is_void(_orient3_theta)) _orient3_theta= pi/6;
180 if (is_void(_orient3_phi)) _orient3_phi= -pi/4;
181
182 func limit3(xn,xx,yn,yx,zn,zx,aspect=)
183 /* DOCUMENT limit3, xmin,xmax, ymin,ymax
184 or limit3, xmin,xmax, ymin,ymax, zmin,zmax
185 Set the 3D axis limits for use with the cage.
186 Use keyword aspect=[ax,ay,az] to set the aspect ratios of the
187 cage to ax:ay:az -- that is, the ratios of the lengths of the
188 cage axes will become ax:ay:az.
189 SEE ALSO: cage3, range3, plwf (include plwf.i), orient3
190 */
191 {
192 local limits;
193 eq_nocopy, limits, _car(_draw3_list, _draw3_nll+1);
194 if (dimsof(xn)(1)>1) {
195 if (dimsof(xn)(1)!=2 || anyof(dimsof(xn)(2:3)!=3))
196 error, "bad limit3 arguments";
197 _setscl3, xn;
198 return limits;
199 }
200 void= [[is_void(xn),is_void(xx)],[is_void(yn),is_void(yx)],
201 [is_void(zn),is_void(zx)]];
202 if (nallof(void) || !is_void(aspect)) {
203 if (is_void(limits)) {
204 if (anyof(void))
205 error, "no xyz limits currently set -- you must set all six";
206 lims= array([0.,1.,1.],3);
207 } else {
208 lims= limits;
209 }
210 if (!void(1,1)) lims(1,1)= xn;
211 if (!void(2,1)) lims(2,1)= xx;
212 if (!void(1,2)) lims(1,2)= yn;
213 if (!void(2,2)) lims(2,2)= yx;
214 if (!void(1,3)) lims(1,3)= zn;
215 if (!void(2,3)) lims(2,3)= zx;
216 if (!is_void(aspect)) lims(3,)= aspect;
217 _setscl3, lims;
218 }
219 return limits;
220 }
221
222 func range3(zn,zx,aspect=)
223 /* DOCUMENT range3, zmin,zmax
224 Set the 3D axis z limits for use with the cage.
225 Use keyword aspect=[ax,ay,az] to set the aspect ratios of the
226 cage to ax:ay:az -- that is, the ratios of the lengths of the
227 cage axes will become ax:ay:az.
228 SEE ALSO: cage3, limit3, plwf (include plwf.i), orient3
229 */
230 {
231 return limit3(,,,,zn,zx,aspect=aspect);
232 }
233
save3(void)234 func save3(void)
235 /* DOCUMENT view= save3()
236 Save the current 3D viewing transformation and lighting.
237 SEE ALSO: restore3, rot3, mov3, aim3, light3
238 */
239 {
240 return _cpy(_draw3_list, _draw3_n);
241 }
242
restore3(view)243 func restore3(view)
244 /* DOCUMENT restore3, view
245 Restore a previously saved 3D viewing transformation and lighting.
246 If VIEW is nil, rotate object to viewer's coordinate system.
247 SEE ALSO: restore3, rot3, mov3, aim3, light3
248 */
249 {
250 if (!is_void(view)) view= _cpy(view);
251 else view= _cat(_cpy(_draw3_view), _cpy(_light3_list));
252 _cdr, view, _draw3_n, _cdr(_draw3_list, _draw3_n, []);
253 _draw3_list= view;
254 _undo3_set, restore3, old;
255 }
256
257 /* set default viewing direction if user hasn't already done so */
258 if (is_void(_draw3_view)) _draw3_view= _lst(unit(3), [0.,0.,0.], []);
259 _draw3_nv= _len(_draw3_view);
260
261 /* ------------------------------------------------------------------------ */
262
263 func light3(ambient=,diffuse=,specular=,spower=,sdir=)
264 /* DOCUMENT light3, ambient=a_level,
265 diffuse=d_level,
266 specular=s_level,
267 spower=n,
268 sdir=xyz
269 Sets lighting properties for 3D shading effects.
270 A surface will be shaded according to its to its orientation
271 relative to the viewing direction.
272
273 The ambient level A_LEVEL is a light level (arbitrary units)
274 that is added to every surface independent of its orientation.
275
276 The diffuse level D_LEVEL is a light level which is proportional
277 to cos(theta), where theta is the angle between the surface
278 normal and the viewing direction, so that surfaces directly
279 facing the viewer are bright, while surfaces viewed edge on are
280 unlit (and surfaces facing away, if drawn, are shaded as if they
281 faced the viewer).
282
283 The specular level S_LEVEL is a light level proportional to a high
284 power spower=N of 1+cos(alpha), where alpha is the angle between
285 the specular reflection angle and the viewing direction. The light
286 source for the calculation of alpha lies in the direction XYZ (a
287 3 element vector) in the viewer's coordinate system at infinite
288 distance. You can have ns light sources by making S_LEVEL, N, and
289 XYZ (or any combination) be vectors of length ns (3-by-ns in the
290 case of XYZ). (See source code for specular_hook function
291 definition if powers of 1+cos(alpha) aren't good enough for you.)
292
293 With no arguments, return to the default lighting.
294
295 EXAMPLES:
296 light3, diffuse=.1, specular=1., sdir=[0,0,-1]
297 (dramatic "tail lighting" effect)
298 light3, diffuse=.5, specular=1., sdir=[1,.5,1]
299 (classic "over your right shoulder" lighting)
300 light3, ambient=.1,diffuse=.1,specular=1.,
301 sdir=[[0,0,-1],[1,.5,1]],spower=[4,2]
302 (two light sources combining previous effects)
303
304 SEE ALSO: rot3, save3, restore3
305 */
306 {
307 old= _cpy(_cdr(_draw3_list,_draw3_nv),5);
308
309 flags= 0;
310 if (!is_void(ambient)) {
311 if (dimsof(ambient)(1)) error, "ambient light level must be scalar";
312 flags|= 1;
313 _car, _draw3_list, _draw3_nv+1, double(ambient);
314 }
315 if (!is_void(diffuse)) {
316 if (dimsof(diffuse)(1)) error, "diffuse light level must be scalar";
317 flags|= 2;
318 _car, _draw3_list, _draw3_nv+2, double(diffuse);
319 }
320
321 if (!is_void(specular)) flags|= 4;
322 else specular= _car(_draw3_list, _draw3_nv+3);
323 if (!is_void(spower)) flags|= 8;
324 else spower= _car(_draw3_list, _draw3_nv+4);
325 if (!is_void(sdir)) {
326 dims= dimsof(sdir);
327 if (dims(1)<1 || dims(2)!=3)
328 error, "lighting direction must be 3 vector or 3-by-ns array"
329 flags|= 16;
330 } else {
331 sdir= _car(_draw3_list, _draw3_nv+5);
332 }
333 if (flags&28) {
334 if (is_void(dimsof(specular,spower,sdir(1,..))))
335 error, "specular, spower, and sdir not conformable";
336 if (flags&4) _car, _draw3_list, _draw3_nv+3, double(specular);
337 if (flags&8) _car, _draw3_list, _draw3_nv+4, spower;
338 if (flags&16) _car, _draw3_list, _draw3_nv+5, double(sdir);
339 }
340
341 if (!flags) {
342 for (i=1 ; i<=5 ; ++i)
343 _car, _draw3_list, _draw3_nv+i, _car(_light3_list,i);
344 }
345
346 _undo3_set, _light3, old;
347 }
348
_light3(arg)349 func _light3(arg)
350 {
351 for (i=1 ; i<=5 ; ++i)
352 _car, _draw3_list, _draw3_nv+i, _car(arg,i);
353 }
354
355 /* set default values if user hasn't already done so */
356 if (is_void(_light3_ambient)) _light3_ambient= 0.2;
357 if (is_void(_light3_diffuse)) _light3_diffuse= 1.0;
358 if (is_void(_light3_specular)) _light3_specular= 0.0;
359 if (is_void(_light3_spower)) _light3_spower= 2;
360 if (is_void(_light3_sdir)) _light3_sdir= [1.0, 0.5, 1.0]/sqrt(2.25);
361 _light3_list= _lst(_light3_ambient,_light3_diffuse,_light3_specular,
362 _light3_spower,_light3_sdir);
363
get3_light(xyz,nxyz)364 func get3_light(xyz, nxyz)
365 /* DOCUMENT get3_light(xyz, nxyz)
366 or get3_light(xyz)
367
368 return 3D lighting for polygons with vertices XYZ. If NXYZ is
369 specified, XYZ should be 3-by-sum(nxyz), with NXYZ being the
370 list of numbers of vertices for each polygon (as for the plfp
371 function). If NXYZ is not specified, XYZ should be a quadrilateral
372 mesh, 3-by-ni-by-nj (as for the plf function). In the first case,
373 the return value is numberof(NXYZ); in the second case, the
374 return value is (ni-1)-by-(nj-1).
375
376 The parameters of the lighting calculation are set by the
377 light3 function.
378
379 SEE ALSO: light3, set3_object, get3_normal, get3_centroid
380 */
381 {
382 list= _cdr(_draw3_list, _draw3_nv);
383 ambient= _nxt(list);
384 diffuse= _nxt(list);
385 specular= _nxt(list);
386 spower= _nxt(list);
387 sdir= _nxt(list);
388
389 /* get normal */
390 normal= get3_normal(xyz, nxyz);
391
392 /* get direction to viewer's eye (camera) */
393 zc= _getzc3();
394 if (is_void(zc)) {
395 view= [0.,0.,1.];
396 } else {
397 view= [0.,0.,zc]-get3_centroid(xyz, nxyz);
398 m1= abs(view(1,..),view(2,..),view(3,..))(-,..);
399 m1= m1 + (m1==0.0);
400 view/= m1;
401 }
402
403 /* do lighting calculation */
404 nv= (normal*view)(sum,..);
405 light= ambient + diffuse*abs(nv);
406 if (anyof(specular)) {
407 sdir= sdir(,*);
408 sdir/= abs(sdir(1,),sdir(2,),sdir(3,))(-,);
409 sv= sdir(+,*)*view(+,..);
410 sn= sdir(+,*)*normal(+,..);
411 m1= max(sn*nv(-,) - 0.5*sv + 0.5, 1.e-30); /* max(1+cos(alpha),0) */
412 if (is_func(specular_hook))
413 m1= specular_hook(m1, abs(nv)(-,..), spower);
414 else
415 m1= m1^spower;
416 light+= (specular(*)*m1)(sum,..);
417 }
418
419 return light;
420 }
421
get3_normal(xyz,nxyz)422 func get3_normal(xyz, nxyz)
423 /* DOCUMENT get3_normal(xyz, nxyz)
424 or get3_normal(xyz)
425
426 return 3D normals for polygons with vertices XYZ. If NXYZ is
427 specified, XYZ should be 3-by-sum(nxyz), with NXYZ being the
428 list of numbers of vertices for each polygon (as for the plfp
429 function). If NXYZ is not specified, XYZ should be a quadrilateral
430 mesh, 3-by-ni-by-nj (as for the plf function). In the first case,
431 the return value is 3-by-numberof(NXYZ); in the second case, the
432 return value is 3-by-(ni-1)-by-(nj-1).
433
434 The normals are constructed from the cross product of the lines
435 joining the midpoints of two edges which as nearly quarter the
436 polygon as possible (the medians for a quadrilateral). No check
437 is made that these not be parallel; the returned "normal" is
438 [0,0,0] in that case. Also, if the polygon vertices are not
439 coplanar, the "normal" has no precisely definable meaning.
440
441 SEE ALSO: get3_centroid, get3_light
442 */
443 {
444 if (is_void(nxyz)) {
445 /* if no polygon list is given, assume xyz is 2D mesh */
446 /* form normal as cross product of medians */
447 m1= xyz(,zcen,dif);
448 m2= xyz(,dif,zcen);
449
450 } else {
451 /* with polygon list, more elaborate calculation required */
452 frst= nxyz(psum)-nxyz+1;
453
454 /* form normal by getting two approximate diameters
455 * (reduces to above medians for quads) */
456 n2= (nxyz+1)/2;
457 zero= frst-1;
458 c0= 0.5*(xyz(,zero+1)+xyz(,zero+2));
459 i= zero+n2; /* n2>=2, nxyz>=3 */
460 c1= 0.5*(xyz(,i)+xyz(,i+1));
461 i= 1+n2/2;
462 c2= 0.5*(xyz(,zero+i)+xyz(,zero+i+1));
463 i= min(i+n2, nxyz);
464 c3= 0.5*(xyz(,zero+i)+xyz(,zero+i%nxyz+1));
465 m1= c1 - c0;
466 m2= c3 - c2;
467 }
468
469 /* poly normal is cross product of two medians (or diameters) */
470 normal= m1;
471 normal(1,..)= n1= m1(2,..)*m2(3,..) - m1(3,..)*m2(2,..);
472 normal(2,..)= n2= m1(3,..)*m2(1,..) - m1(1,..)*m2(3,..);
473 normal(3,..)= n3= m1(1,..)*m2(2,..) - m1(2,..)*m2(1,..);
474 m1= abs(n1,n2,n3)(-,..);
475 m1= m1 + (m1==0.0);
476 normal/= m1;
477
478 return normal;
479 }
480
get3_centroid(xyz,nxyz)481 func get3_centroid(xyz, nxyz)
482 /* DOCUMENT get3_centroid(xyz, nxyz)
483 or get3_centroid(xyz)
484
485 return 3D centroids for polygons with vertices XYZ. If NXYZ is
486 specified, XYZ should be 3-by-sum(nxyz), with NXYZ being the
487 list of numbers of vertices for each polygon (as for the plfp
488 function). If NXYZ is not specified, XYZ should be a quadrilateral
489 mesh, 3-by-ni-by-nj (as for the plf function). In the first case,
490 the return value is 3-by-numberof(NXYZ); in the second case, the
491 return value is 3-by-(ni-1)-by-(nj-1).
492
493 The centroids are constructed as the mean value of all vertices
494 of each polygon.
495
496 SEE ALSO: get3_normal, get3_light
497 */
498 {
499 if (is_void(nxyz)) {
500 /* if no polygon list is given, assume xyz is 2D mesh */
501 centroid= xyz(,zcen,zcen);
502
503 } else {
504 /* with polygon list, more elaborate calculation required */
505 last= nxyz(psum);
506 list= histogram(1+last)(1:-1);
507 list(1)+= 1;
508 list= list(psum);
509 centroid= array(0.0, 3, numberof(nxyz));
510 centroid(1,)= histogram(list, xyz(1,));
511 centroid(2,)= histogram(list, xyz(2,));
512 centroid(3,)= histogram(list, xyz(3,));
513 centroid/= double(nxyz);
514 }
515
516 return centroid;
517 }
518
519 /* ------------------------------------------------------------------------ */
520
521 func get3_xy(xyz, &x, &y, &z, getz)
522 /* DOCUMENT get3_xy, xyz, x, y
523 or get3_xy, xyz, x, y, z, 1
524
525 Given 3-by-anything coordinates XYZ, return X and Y in viewer's
526 coordinate system (set by rot3, mov3, orient3, etc.). If the
527 fifth argument is present and non-zero, also return Z (for use
528 in sort3d or get3_light, for example). If the camera position
529 has been set to a finite distance with setz3, the returned
530 coordinates will be tangents of angles for a perspective
531 drawing (and Z will be scaled by 1/zc).
532
533 SEE ALSO: sort3d, get3_light, rot3, setz3, set3_object
534 */
535 {
536 /* rotate and translate to viewer's coordinate system */
537 xyz= _getrot3()(,+)*(_getscl3()*xyz - _getorg3())(+,..);
538 x= xyz(1,..);
539 y= xyz(2,..);
540
541 /* do optional perspective projection */
542 zc= _getzc3();
543 if (!is_void(zc)) {
544 z= xyz(3,..);
545 zc= max(zc-z, 0.0); /* protect behind camera */
546 zc+= (zc==0.0)*1.e-35; /* avoid zero divide */
547 x/= zc;
548 y/= zc;
549 if (getz) z/= zc;
550 } else if (getz) {
551 z= xyz(3,..);
552 }
553 }
554
_getrot3(void)555 func _getrot3(void)
556 {
557 return _car(_draw3_list, 1);
558 }
_getorg3(void)559 func _getorg3(void)
560 {
561 local limits, org3;
562 eq_nocopy, limits, _car(_draw3_list, _draw3_nll+1);
563 eq_nocopy, org3, _car(_draw3_list, 2);
564 if (is_void(limits)) return org3;
565 return org3 + limits(avg:1:2,)*limits(3,)/limits(ptp:1:2,);
566 }
_getzc3(void)567 func _getzc3(void)
568 {
569 return _car(_draw3_list, 3);
570 }
_getscl3(void)571 func _getscl3(void)
572 {
573 local limits;
574 eq_nocopy, limits, _car(_draw3_list, _draw3_nll+1);
575 if (is_void(limits)) return [1.,1.,1.];
576 return limits(3,)/limits(ptp:1:2,);
577 }
578
_setrot3(x)579 func _setrot3(x)
580 {
581 _undo3_set, _setrot3, _car(_draw3_list, 1, x);
582 }
_setorg3(x)583 func _setorg3(x)
584 {
585 _undo3_set, _setorg3, _car(_draw3_list, 2, x);
586 }
_setzc3(x)587 func _setzc3(x)
588 {
589 _undo3_set, _setzc3, _car(_draw3_list, 3, x);
590 }
_setscl3(x)591 func _setscl3(x)
592 { /* not undoable, cleared by clear3 */
593 _car, _draw3_list, _draw3_nll+1, x;
594 draw3_trigger;
595 }
596
_undo3_set(fnc,arg)597 func _undo3_set(fnc, arg)
598 {
599 if (!_in_undo3) {
600 if (_len(_undo3_list)>=2*_undo3_limit)
601 _cdr, _undo3_list, 2*_undo3_limit-2, [];
602 _undo3_list= _cat(_lst(fnc,arg), _undo3_list);
603 }
604 draw3_trigger;
605 }
606
607 _undo3_limit= 100;
608
undo3(n)609 func undo3(n)
610 /* DOCUMENT undo3
611 or undo3, n
612 Undo the effects of the last N (default 1) rot3, orient3, mov3, aim3,
613 setz3, or light3 commands.
614 */
615 {
616 if (is_void(n)) n= 1;
617 n= 2*(n-1);
618 if (n<0 || n>_len(_undo3_list))
619 error, "not that many items in undo list";
620 _in_undo3= 1; /* flag to skip _undo3_set */
621 /* perhaps should save discarded items in a redo list? */
622 if (n) _undo3_list= _cdr(_undo3_list, n);
623 for (; n>=0 ; n-=2) {
624 fnc= _nxt(_undo3_list);
625 arg= _nxt(_undo3_list);
626 fnc, arg;
627 }
628 draw3_trigger;
629 }
630
set3_object(fnc,arg)631 func set3_object(fnc, arg)
632 /* DOCUMENT set3_object, drawing_function, _lst(arg1,arg2,...)
633
634 set up to trigger a call to draw3, adding a call to the
635 3D display list of the form:
636
637 DRAWING_FUNCTION, _lst(ARG1, ARG2, ...)
638
639 When draw3 calls DRAWING_FUNCTION, the external variable _draw3
640 will be non-zero, so DRAWING_FUNCTION can be written like this:
641
642 func drawing_function(arg1,arg2,...)
643 {
644 require, "pl3d.i";
645 if (_draw3) {
646 list= arg1;
647 arg1= _nxt(list);
648 arg2= _nxt(list);
649 ...
650 ...<calls to get3_xy, sort3d, get3_light, etc.>...
651 ...<calls to graphics functions plfp, plf, etc.>...
652 return;
653 }
654 ...<verify args>...
655 ...<do orientation and lighting independent calcs>...
656 set3_object, drawing_function, _lst(arg1,arg2,...);
657 }
658
659 SEE ALSO: get3_xy, get3_light, sort3d
660 */
661 {
662 _draw3_list= _cat(_draw3_list, _lst(fnc,arg));
663 draw3_trigger;
664 }
665
666 /* ------------------------------------------------------------------------ */
667
draw3(called_as_idler)668 func draw3(called_as_idler)
669 /* DOCUMENT draw3
670 Draw the current 3D display list.
671 (Ordinarily triggered automatically when the drawing changes.)
672 */
673 {
674 if (_draw3_changes) {
675 if (called_as_idler) fma;
676 /* the first _draw3_n elements of _draw3_list are the viewing
677 * transforms, lighting, etc.
678 * thereafter, elements are (function,argument-list) pairs
679 * the _draw3 flag alerts the functions that these are the draw
680 * calls rather than the interactive setup calls */
681 limits, square=1;
682 if (_cage3) {
683 local lims;
684 eq_nocopy, lims, _car(_draw3_list, _draw3_nll+1);
685 if (!is_void(lims)) _3cage, transpose(lims(1:2,));
686 }
687 _draw3= 1;
688 for (list=_cdr(_draw3_list, _draw3_n) ; list ; list=_cdr(list)) {
689 fnc= _car(list);
690 list= _cdr(list);
691 fnc, _car(list);
692 }
693 if (_gnomon) _gnomon_draw;
694 _draw3_changes= [];
695 }
696 }
697
698 _draw3_nll= _draw3_nv+_len(_light3_list);
699 _draw3_list= _cat(_cpy(_draw3_view), _cpy(_light3_list), _lst([]));
700 _draw3_n= _len(_draw3_list);
701
702 func draw3_trigger
703 {
704 /* arrange to call draw3 when everything else is finished */
705 set_idler, _draw3_idler;
706 extern _draw3_changes;
707 _draw3_changes= 1;
708 }
709
710 func _draw3_idler
711 {
712 draw3, 1;
713 }
714
clear3(void)715 func clear3(void)
716 /* DOCUMENT clear3
717 Clear the current 3D display list.
718 */
719 {
720 _cdr, _draw3_list, _draw3_n, [];
721 _car, _draw3_list, _draw3_nll+1, []; /* clear limits */
722 }
723
window3(n)724 func window3(n)
725 /* DOCUMENT window3
726 or window3, n
727 initialize style="nobox.gs" window for 3D graphics
728 */
729 {
730 window, n, wait=1, style="nobox.gs", legends=0;
731 }
732
733 /* ------------------------------------------------------------------------ */
734
gnomon(on)735 func gnomon(on)
736 /* DOCUMENT gnomon
737 or gnomon, onoff
738 Toggle the gnomon display. If ONOFF is non-nil and non-zero,
739 turn on the gnomon. If ONOFF is zero, turn off the gnomon.
740
741 The gnomon shows the X, Y, and Z axis directions in the
742 object coordinate system. The directions are labeled.
743 The gnomon is always infinitely far behind the object
744 (away from the camera).
745
746 There is a mirror-through-the-screen-plane ambiguity in the
747 display which is resolved in two ways: (1) The (X,Y,Z)
748 coordinate system is right-handed, and (2) If the tip of an
749 axis projects into the screen, it's label is drawn in opposite
750 polarity to the other text on the screen.
751 */
752 {
753 old= _gnomon;
754 if (is_void(on)) _gnomon~= 1;
755 else if (on) _gnomon= 1;
756 else _gnomon= 0;
757 if (old!=_gnomon) draw3_trigger;
758 }
759
760 if (is_void(_gnomon)) _gnomon= 0;
761
_gnomon_draw(void)762 func _gnomon_draw(void)
763 {
764 o= [0.,0.,0.];
765 x1= [1.,0.,0.];
766 y1= [0.,1.,0.];
767 z1= [0.,0.,1.];
768 xyz= _getrot3()(,+)*[[o,x1],[o,y1],[o,z1]](+,,,);
769 xyz*= 0.0013*_gnomon_scale;
770 x1= xyz(1,,);
771 y1= xyz(2,,);
772 z1= xyz(3,2,);
773 x0= x1(1,);
774 x1= x1(2,);
775 y0= y1(1,);
776 y1= y1(2,);
777 wid= min(_gnomon_scale/18.,6.);
778 if (wid<0.5) wid= 0.0;
779 plsys, 0;
780 pldj, x0+_gnomon_x, y0+_gnomon_y, x1+_gnomon_x, y1+_gnomon_y,
781 width=wid, type=1, legend=string(0);
782 plsys, 1;
783
784 /* compute point size of labels (1/3 of axis length) */
785 pts= [8,10,12,14,18,24](digitize(_gnomon_scale/3.0,
786 [9,11,13,16,21]));
787 if (_gnomon_scale < 21.) {
788 x1*= 21./_gnomon_scale;
789 y1*= 21./_gnomon_scale;
790 }
791 /* label positions: first find shortest axis */
792 xy= abs(x1,y1);
793 i= xy(mnx);
794 jk= [[2,3],[3,1],[1,2]](,i);
795 if (xy(i)<1.e-7*xy(sum)) { /* guarantee not exactly zero */
796 x1(i)= -1.e-6*x1(jk)(sum);
797 y1(i)= -1.e-6*y1(jk)(sum);
798 xy(i)= abs(x1(i),y1(i));
799 }
800 xyi= xy(i);
801 /* next find axis nearest to shortest */
802 j= jk(1);
803 k= jk(2);
804 if (abs(x1(j)*y1(i)-y1(j)*x1(i))*xy(k) >
805 abs(x1(k)*y1(i)-y1(k)*x1(i))*xy(j)) {
806 jk= j; j= k; k= jk;
807 }
808 /* furthest axis first - move perpendicular to nearer axis */
809 xk= -y1(j);
810 yk= x1(j);
811 xy= abs(xk,yk);
812 xk/= xy;
813 yk/= xy;
814 if (xk*x1(k)+yk*y1(k) < 0.0) { xk= -xk; yk= -yk; }
815 /* nearer axis next - move perpendicular to furthest axis */
816 xj= -y1(k);
817 yj= x1(k);
818 xy= abs(xj,yj);
819 xj/= xy;
820 yj/= xy;
821 if (xj*x1(j)+yj*y1(j) < 0.0) { xj= -xj; yj= -yj; }
822 /* shortest axis last - move perpendicular to nearer */
823 xi= -y1(j);
824 yi= x1(j);
825 xy= abs(xi,yi);
826 xi/= xy;
827 yi/= xy;
828 if (xi*x1(i)+yi*y1(i) < 0.0) { xi= -xi; yi= -yi; }
829
830 /* shortest axis label may need adjustment */
831 d= 0.0013*pts;
832 if (xyi < d) {
833 /* just center it in correct quadrant */
834 jk= sign(xi*xj+yi*yj);
835 yi= sign(xi*xk+yi*yk);
836 xi= jk*xj + yi*xk;
837 yi= jk*yj + yi*yk;
838 jk= abs(xi, yi);
839 xi/= jk;
840 yi/= jk;
841 }
842
843 x= y= [0.,0.,0.];
844 x([i,j,k])= [xi,xj,xk];
845 y([i,j,k])= [yi,yj,yk];
846 x*= d;
847 y*= d;
848 x+= x1 + _gnomon_x;
849 y+= y1 + _gnomon_y;
850 chr= ["X","Y","Z"];
851 _gnomon_text, chr(i), x(i),y(i), pts, z1(i)<-1.e-6;
852 _gnomon_text, chr(j), x(j),y(j), pts, z1(j)<-1.e-6;
853 _gnomon_text, chr(k), x(k),y(k), pts, z1(k)<-1.e-6;
854 }
855
856 /* recommended _gnomon_scale: 24, 30, 36, 42, 54, or 72 */
857 if (is_void(_gnomon_scale))
858 _gnomon_scale= 30. /* X,Y,Z axis lengths in points */
859 if (is_void(_gnomon_x))
860 _gnomon_x= 0.18; /* gnomon origin in system 0 coordinates */
861 if (is_void(_gnomon_y))
862 _gnomon_y= 0.42;
863
_gnomon_text(chr,x,y,pts,invert)864 func _gnomon_text(chr, x, y, pts, invert)
865 {
866 /* pts= 8, 10, 12, 14, 18, or 24 */
867 col= "fg";
868 if (invert) {
869 plsys, 0;
870 plg,[y,y],[x,x],
871 type=1,width=2.2*pts,color=col,marks=0,legend=string(0);
872 plsys, 1;
873 col= "bg";
874 }
875 plt, chr, x,y,
876 justify="CH",color=col,height=pts,font="helvetica",opaque=0;
877 }
878
879 /* ------------------------------------------------------------------------ */
880
881 func sort3d(z, npolys, &list, &vlist)
882 /* DOCUMENT sort3d(z, npolys, &list, &vlist)
883
884 given Z and NPOLYS, with numberof(Z)==sum(npolys), return
885 LIST and VLIST such that Z(VLIST) and NPOLYS(LIST) are
886 sorted from smallest average Z to largest average Z, where
887 the averages are taken over the clusters of length NPOLYS.
888 Within each cluster (polygon), the cyclic order of Z(VLIST)
889 remains unchanged, but the absolute order may change.
890
891 This sorting order produces correct or nearly correct order
892 for a plfp command to make a plot involving hidden or partially
893 hidden surfaces in three dimensions. It works best when the
894 polys form a set of disjoint closed, convex surfaces, and when
895 the surface normal changes only very little between neighboring
896 polys. (If the latter condition holds, then even if sort3d
897 mis-orders two neighboring polys, their colors will be very
898 nearly the same, and the mistake won't be noticeable.) A truly
899 correct 3D sorting routine is impossible, since there may be no
900 rendering order which produces correct surface hiding (some polys
901 may need to be split into pieces in order to do that). There
902 are more nearly correct algorithms than this, but they are much
903 slower.
904
905 SEE ALSO: get3_xy
906 */
907 {
908 /* first compute z, the z-centroid of every poly
909 * get a list the same length as x, y, or z which is 1 for each
910 * vertex of poly 1, 2 for each vertex of poly2, etc.
911 * the goal is to make nlist with histogram(nlist)==npolys */
912 nlist= histogram(1+npolys(psum))(1:-1);
913 nlist(1)+= 1; /* another problem with 1-origin indexing */
914 nlist= nlist(psum);
915 /* now sum the vertex values and divide by the number of vertices */
916 z= histogram(nlist, double(z))/npolys;
917
918 /* sort the polygons from smallest z to largest z
919 * npolys(list) is the sorted list of lengths */
920 list= sort(z);
921
922 /* next, find the list which sorts the polygon vertices
923 * first, find a list vlist such that sort(vlist) is above list */
924 vlist= list;
925 vlist(list)= indgen(numberof(list));
926 /* then reset the nlist values to that pre-sorted order, so that
927 * sort(nlist) will be the required vertex sorting list */
928 nlist= vlist(nlist);
929 /* the final hitch is to ensure that the vertices within each polygon
930 * remain in their initial order (sort scrambles equal values)
931 * since the vertices of a polygon can be cyclically permuted,
932 * it suffices to add a sawtooth function to a scaled nlist to
933 * produce a list in which each cluster of equal values will retain
934 * the same cyclic order after the sort
935 * (note that the more complicated msort routine would leave the
936 * clusters without even a cyclic permutation, if that were
937 * necessary) */
938 nmax= max(npolys); /* this must never be so large that
939 * numberof(npolys)*nmax > 2e9 */
940 vlist= sort(nmax*nlist + indgen(numberof(nlist))%nmax);
941 /* primary sort key ^ secondary key ^ */
942 }
943
944 /* ------------------------------------------------------------------------ */
945
946 func spin3(nframes, axis, tlimit=, dtmin=, bracket_time=)
947 /* DOCUMENT spin3
948 or spin3, nframes
949 or spin3, nframes, axis
950 Spin the current 3D display list about AXIS over NFRAMES. Keywords
951 tlimit= the total time allowed for the movie in seconds (default 60),
952 dtmin= the minimum allowed interframe time in seconds (default 0.0),
953 bracket_time= (as for movie function in movie.i)
954
955 The default AXIS is [-1,1,0] and the default NFRAMES is 30.
956
957 SEE ALSO: rot3
958 */
959 {
960 require, "movie.i";
961 if (is_void(nframes)) nframes= 30;
962 dtheta= 2*pi/(nframes-1);
963 if (is_void(axis)) axis= [-1.,1.,0.];
964 theta= acos(axis(3)/abs(axis(1),axis(2),axis(3)));
965 phi= atan(axis(2),axis(1)+(!axis(1)&&!axis(2)));
966 orig= save3();
967 movie, _spin3, tlimit, dtmin, bracket_time;
968 restore3, orig;
969 }
970
_spin3(i)971 func _spin3(i)
972 {
973 if (i>=nframes) return 0;
974 rot3,,,-phi
975 rot3,,-theta,dtheta;
976 rot3,,theta,phi;
977 draw3;
978 return 1;
979 }
980
981 /* ------------------------------------------------------------------------ */
982
cage3(on)983 func cage3(on)
984 /* DOCUMENT cage3
985 or cage3, onoff
986 Toggle the cage display. If ONOFF is non-nil and non-zero,
987 turn on the cage. If ONOFF is zero, turn off the cage.
988
989 The cage draws a rectangular box "behind" the 3D object and
990 attempts to put ticks and labels around the edge of the box.
991
992 SEE ALSO: limit3, plwf (include plwf.i)
993 */
994 {
995 old= _cage3;
996 if (is_void(on)) _cage3~= 1;
997 else if (on) _cage3= 1;
998 else _cage3= 0;
999 if (old!=_cage3) draw3_trigger;
1000 }
1001
1002 if (is_void(_cage3)) _cage3= 0;
1003
1004 /*
1005 Idea:
1006
1007 Draw hexagon, representing the free edges of the three "backplanes"
1008 of a cube surrounding the 3D object being plotted. Ticks should
1009 be perpendicular to the backplane associated with the edge. Try
1010 to put numeric labels on the "leftmost" and "bottommost" edges,
1011 while "X", "Y", and "Z" markers can go on the"topmost" and
1012 "rightmost" edges. In the pl3d.i interface, this "ticked box"
1013 should optionally replace the gnomon, since there is no reason
1014 to have both. Should also draw (optionally?) the common edges of
1015 the three backplanes.
1016 */
1017
_3cage(xyzlim)1018 func _3cage(xyzlim)
1019 {
1020 /* make 3x2x2x2 = (vector,dx,dy,dz) coordinates of the box,
1021 starting from 3x2 (vector,min_max) limits values */
1022 dxyz= xyzlim(,2)-xyzlim(,1);
1023 xyz= xyzlim(,1,-:1:2,-:1:2,-:1:2);
1024 xyz(1,2,,)+= dxyz(1);
1025 xyz(2,,2,)+= dxyz(2);
1026 xyz(3,,,2)+= dxyz(3);
1027
1028 /* transform to three 2x2x2 arrays (x,y,z) in camera coord system */
1029 local x, y, z;
1030 get3_xy, xyz, x, y, z, 1;
1031
1032 /* find the three orthogonal backplanes
1033 owing to perspective, there may be up to five back facing planes
1034 of the box -- choose the most rearward, or the smaller object
1035 coordinate value (xn,yn,zn) in case of ties */
1036 i= z(,avg,avg)(mnx);
1037 j= z(avg,,avg)(mnx);
1038 k= z(avg,avg,)(mnx);
1039 /* now the common point of the backplanes is (i,j,k), the three common
1040 edges connect to the points (3-i,j,k), (i,3-j,k), and (i,j,3-k),
1041 and the three non-common points of the backplanes are
1042 (3-i,j,3-k), (3-i,3-j,k), and (i,3-j,3-k) */
1043
1044 /* it is best to draw these as polylines, since polylines have
1045 round endcaps, while disjoint lines have square caps
1046 -- if the common edges are omitted, there is just a single
1047 closed hexagon. otherwise, we draw two open polylines,
1048 one starting at the common point then traversing the
1049 hexagon (8 points), the other consisting of the two
1050 remaining common edges (3 points) */
1051 ip= 3-i;
1052 j= 2*(j-1);
1053 jp= 2-j;
1054 k= 4*(k-1);
1055 kp= 4-k;
1056 if (omit3_common)
1057 list= [ip+j+k, ip+jp+k, i+jp+k, i+jp+kp, i+j+kp, ip+j+kp];
1058 else
1059 list= [i+j+k, ip+j+k, ip+jp+k, i+jp+k, i+jp+kp, i+j+kp, ip+j+kp, ip+j+k];
1060 x1= x(list);
1061 y1= y(list);
1062 plg, y1, x1, closed=omit3_common, marks=0, legend="",
1063 color=_3kcolor, type=_3ktype, width=_3kwidth;
1064 xmin= min(x1);
1065 xmax= max(x1);
1066 ymin= min(y1);
1067 ymax= max(y1);
1068
1069 if (!omit3_common) {
1070 list= [i+jp+k, i+j+k, i+j+kp];
1071 plg, y(list), x(list), closed=0, marks=0, legend="",
1072 color=_3kcolor, type=_3ktype, width=_3kwidth;
1073 }
1074
1075 nlabel= 0;
1076 if (!omit3_ticks) {
1077 /* innermost index is [back,front] point,
1078 next index is [first, second] edge,
1079 third index is [x,y,z] */
1080 list= [[[i+j+kp,ip+j+kp],[i+jp+k,ip+jp+k]],
1081 [[i+j+kp,i+jp+kp],[ip+j+k,ip+jp+k]],
1082 [[ip+j+k,ip+j+kp],[i+jp+k,i+jp+kp]]];
1083 /* reorder first index so that first point is minimum */
1084 n= list(2,1,)<list(1,1,);
1085 for (m=1 ; m<=3 ; ++m) if (n(m)) list(,,m)= list(2:1:-1,,m);
1086 /* reorder second index so that first edge is the one which
1087 will get any labels -- leftmost or lower */
1088 x1= x(list);
1089 y1= y(list);
1090 horiz= (abs(y1(ptp,,))>abs(x1(ptp,,)));
1091 x2= x1(avg,,);
1092 y2= y1(avg,,);
1093 dir= sign(x2-x2(2:1:-1,))*horiz + sign(y2-y2(2:1:-1,))*(!horiz);
1094 for (m=1 ; m<=3 ; ++m) {
1095 if (dir(1,m)<0.) {
1096 if (dir(2,m)>0. || horiz(2,m)) continue;
1097 } else {
1098 if (dir(2,m)>0. && horiz(1,m)) continue;
1099 }
1100 list(,,m)= list(,2:1:-1,m);
1101 horiz(,m)= horiz(2:1:-1,m);
1102 dir(,m)= dir(2:1:-1,m);
1103 }
1104 xyz0= xyz(,i+j+k); /* back corner */
1105 xyz= xyz(,list); /* 3x2x2x3 edge coordinates, with last
1106 three indices as above */
1107
1108 /* compute projected lengths of the labeled edges */
1109 n= list(,1,);
1110 len= abs(x(n)(ptp,),y(n)(ptp,),z(n)(ptp,));
1111 mxlen= max(len);
1112 len/= mxlen; /* ...relative to longest edge */
1113
1114 /* identify the two shared corners of the labeled edges --
1115 list(,1,) is 6 numbers, 2 are unique and 2 others shared;
1116 both endpoints of one of the three edges are shared */
1117 ms1= (n(*)==n(-,,))(sum,,)-1; /* 0 for unique, 1 for shared pt */
1118 ms1= where(ms1); /* four indices, three possible pairings */
1119 if (n(ms1(1))==n(ms1(2))) {
1120 ms2= [ms1(3),ms1(4)];
1121 ms1= [ms1(1),ms1(2)];
1122 } else if (n(ms1(1))==n(ms1(3))) {
1123 ms2= [ms1(2),ms1(4)];
1124 ms1= [ms1(1),ms1(3)];
1125 } else {
1126 ms2= [ms1(2),ms1(3)];
1127 ms1= [ms1(1),ms1(4)];
1128 }
1129
1130 /* First idea was to make the ticks lie along the perpendicular
1131 (in 3D) to the backplane of the associated edge. This looks
1132 slick for orientations in which the hexagon is near regular,
1133 but very ugly when some angles project to near 0 or 180
1134 degrees.
1135
1136 The second idea is more robust, but not quite as nice at its
1137 best: ticks are always either horizontal or vertical (after
1138 projection), with the choice made according to whether the
1139 edge is more nearly vertical or horizontal (respectively).
1140 */
1141
1142 local xmajor, xminor, xlabel;
1143 unit= array(0.,3,3);
1144 unit(1:9:4,1)= 1.;
1145 llen= [0,0,0];
1146 xptr= yptr= lptr= array(pointer, 3);
1147 xylabs= array(0.,2,2,3);
1148 for (m=1 ; m<=3 ; ++m) {
1149 /* longest labeled edge gets specified maximum number of
1150 major ticks -- other two edges get scaled number */
1151 if (_3ticks(xyzlim(m,1), xyzlim(m,2), _3nmajor*len(m)))
1152 continue; /* (no ticks on very short axes) */
1153
1154 /* expand coordinate lists to vector lists */
1155 mask= unit(,m);
1156 xmajor= mask*xmajor(-,);
1157 xminor= mask*xminor(-,);
1158 mask= 1.-mask;
1159
1160 /* do upper or right edge, then lower or left edge */
1161 for (n=2 ; n>=1 ; --n) {
1162 xyz0= xyz(,1,n,m)*mask;
1163
1164 /* first do minor ticks */
1165 tend= xyz0 + xminor;
1166 get3_xy, tend, x1, y1;
1167 if (horiz(n,m)) tdir= [mxlen*dir(n,m),0.];
1168 else tdir= [0.,mxlen*dir(n,m)];
1169 x2= x1 + tdir(1)*_3lminor;
1170 y2= y1 + tdir(2)*_3lminor;
1171 pldj, x1,y1,x2,y2, legend="",
1172 type=_3ktype,color=_3kcolor,width=_3kwidth;
1173 xmin= min(min(x2),xmin);
1174 xmax= max(max(x2),xmax);
1175 ymin= min(min(y2),ymin);
1176 ymax= max(max(y2),ymax);
1177
1178 /* then do major ticks */
1179 tend= xyz0 + xmajor;
1180 get3_xy, tend, x1, y1;
1181 x2= x1 + tdir(1)*_3lmajor;
1182 y2= y1 + tdir(2)*_3lmajor;
1183 pldj, x1,y1,x2,y2, legend="",
1184 type=_3ktype,color=_3kcolor,width=_3kwidth;
1185 xmin= min(min(x2),xmin);
1186 xmax= max(max(x2),xmax);
1187 ymin= min(min(y2),ymin);
1188 ymax= max(max(y2),ymax);
1189 }
1190
1191 /* save the tips of the major ticks for the labeled edge */
1192 xylabs(,,m)= [[x2(1),y2(1)],[x2(0),y2(0)]];
1193 xptr(m)= &x2;
1194 yptr(m)= &y2;
1195 lptr(m)= &xlabel;
1196 llen(m)= max(strlen(xlabel));
1197 nlabel+= numberof(xlabel);
1198 }
1199
1200 if (!omit_labels && nlabel) {
1201 /* estimate horizontal and vertical size of labels in
1202 world coordinates -- obviously this makes assumptions about
1203 the size of the viewport and won't work if the cage is zoomed
1204 to a very different size
1205 the default numbers assume nobox.gs style (6 inch viewport) */
1206 scale= 6. /*inches/viewport*/ * 72.27 /*points/inch*/ /
1207 (2.2*mxlen) /*units/viewport*/; /* net: points/unit */
1208 scale*= _3xfudge; /* get serious */
1209 tvert= _3xheight/scale; /* text size in world coordinates */
1210 thoriz= 0.6*tvert*llen; /* assume text shape about 9x15 */
1211
1212 /* locate and discard any labels
1213 which interfere near the shared corners */
1214 hflags= horiz(1,);
1215 _3interference, ms1;
1216 _3interference, ms2;
1217
1218 dir= dir(1,);
1219 for (m=1 ; m<=3 ; ++m) {
1220 horiz= hflags(m);
1221 tdir= dir(m);
1222 offset= _3xoffset*mxlen*tdir;
1223 eq_nocopy, x2, *xptr(m);
1224 eq_nocopy, y2, *yptr(m);
1225 eq_nocopy, xlabel, *lptr(m);
1226 if (numberof(xlabel)) {
1227 if (horiz) x2+= offset;
1228 else y2+= offset; /* add term in thoriz(m)? */
1229 if (horiz) justify= (tdir<0.)? 15 : 13;
1230 else justify= (tdir<0.)? 10 : 18;
1231 for (n=1 ; n<=numberof(xlabel) ; ++n) {
1232 plt, xlabel(n), x2(n),y2(n), tosys=1, justify=justify,
1233 color=_3xcolor,font=_3xfont,height=_3xheight;
1234 }
1235 if (horiz) {
1236 xmin= min(min(x2-thoriz(m)),xmin);
1237 } else {
1238 ymin= min(min(y2-tvert),ymin);
1239 xmin= min(min(x2-0.5*thoriz(m)),xmin);
1240 xmax= max(max(x2+0.5*thoriz(m)),xmax);
1241 }
1242 }
1243 }
1244
1245 /* autoscaling will clip text */
1246 dx= xmax-xmin;
1247 dy= ymax-ymin;
1248 dd= 0.5*(dx-dy);
1249 if (dd>0.) {
1250 ymax+= dd;
1251 ymin-= dd;
1252 } else {
1253 xmax-= dd;
1254 xmin+= dd;
1255 }
1256 limits, xmin,xmax,ymin,ymax;
1257 }
1258 }
1259 }
1260
1261 /* tick parameters */
1262 _3nmajor= 4.5; /* 1/(max allowed tick density on longest edge) */
1263 _3lmajor= 0.06; /* length of major ticks as fraction of longest edge */
1264 _3lminor= 0.03; /* length of minor ticks as fraction of longest edge */
1265 _3ktype= 1;
1266 _3kwidth= 0.0;
1267 _3kcolor= -2;
1268 /* label parameters */
1269 _3xoffset= 0.02; /* label offset from tick as fraction of longest edge */
1270 _3xfudge= 1.0; /* fudge factor for points/(longest edge) calc */
1271 _3xfont= 8; /* helvetica */
1272 _3xheight= 14.0 /* point size of label text */
1273 _3xcolor= -2;
1274
_3interference(ms)1275 func _3interference(ms)
1276 {
1277 ms-= 1;
1278 m= ms/2 + 1;
1279 m1= m(1);
1280 m2= m(2);
1281 horiz= hflags(m1);
1282
1283 /* interference only possible if ticks have like orientation */
1284 if (hflags(m2)!=horiz) return;
1285
1286 n= 1-(ms%2); /* 1 for 1st index, 0 for last index */
1287 n1= n(1);
1288 n2= n(2);
1289
1290 local x1, y1, x2, y2;
1291 if (horiz) {
1292 eq_nocopy, y1, *yptr(m1);
1293 eq_nocopy, y2, *yptr(m2);
1294 if (!is_void(y1) && !is_void(y2))
1295 interfere= (abs(y1(n1)-y2(n2))<tvert);
1296 } else {
1297 eq_nocopy, x1, *xptr(m1);
1298 eq_nocopy, x2, *xptr(m2);
1299 if (!is_void(x1) && !is_void(x2))
1300 interfere= (abs(x1(n1)-x2(n2))<0.5*(thoriz(m1)+thoriz(m2)));
1301 }
1302
1303 if (interfere) {
1304 _3remove, m1, n1;
1305 _3remove, m2, n2;
1306 }
1307 }
1308
_3remove(m,n,a)1309 func _3remove(m, n, a)
1310 {
1311 if (!is_void(m)) {
1312 xptr(m)= &_3remove(,n,*xptr(m));
1313 yptr(m)= &_3remove(,n,*yptr(m));
1314 lptr(m)= &_3remove(,n,*lptr(m));
1315 } else if (numberof(a)>1) {
1316 if (n) return a(2:0);
1317 else return a(1:-1);
1318 }
1319 }
1320
_3ticks(xmin,xmax,n)1321 func _3ticks(xmin, xmax, n)
1322 {
1323 extern xmajor, xminor, xlabel; /* results */
1324 xmajor= xminor= xlabel= [];
1325
1326 /* dx is the minimum allowed spacing between ticks */
1327 if (n<1.0) return 1;
1328 dx= abs(xmax-xmin)/double(n);
1329 sdx= sign(xmax-xmin);
1330 xmin*= sdx;
1331 xmax*= sdx;
1332
1333 /* round dx up to the nearest "nice" value */
1334 pwr= 10.^floor(log10(dx)+0.0001);
1335 base= digitize(dx/pwr,[2.0001,5.0001]);
1336 if (base==3) pwr*= 10.;
1337 dx= [2.,5.,1.](base)*pwr;
1338
1339 /* find the major tick values */
1340 xn= ceil(xmin/dx - 0.0001);
1341 xx= floor(xmax/dx + 0.0001);
1342 if (xx<=xn) return 1;
1343 xmajor= (xn+indgen(0:long(xx-xn)))*dx*sdx;
1344 /* sigh -- don't want to print "-0" */
1345 xn= where(xmajor==0.);
1346 if (numberof(xn)) xmajor(xn(1))= 0.0;
1347
1348 /* find the minor tick values */
1349 dxn= dx/(base==1? 4. : 5.);
1350 xn= ceil(xmin/dxn - 0.0001);
1351 xx= floor(xmax/dxn + 0.0001);
1352 xminor= (xn+indgen(0:long(xx-xn)))*dxn*sdx;
1353 /* remove minor ticks which are also major */
1354 xminor= xminor(where(abs(floor(xminor/dxn+0.01)*dxn -
1355 floor(xminor/dx+0.01)*dx) > 0.01*dx));
1356
1357 /* compute major tick labels */
1358 xn= abs(xmajor);
1359 xx= max(xn);
1360 ipwr= floor(log10(xx)+0.0001);
1361 xpwr= 10.^ipwr;
1362 ipwr= long(ipwr);
1363 if (ipwr>3 || min(xn+(!xn))<0.00099999999) {
1364 /* use e format -- normalize to 0<= xmajor <10 */
1365 npwr= xpwr;
1366 mpwr= ipwr;
1367 pwr/= npwr;
1368 ipwr= 0;
1369 xpwr= 1.0;
1370 } else {
1371 /* use f format */
1372 npwr= 1.0;
1373 mpwr= 0;
1374 }
1375 ndecimals= max(long(floor(log10(xpwr/pwr)+0.0001)) - ipwr, 0);
1376 format= swrite(format="%%.%ldf",ndecimals);
1377 xlabel= swrite(format=format, xmajor/npwr);
1378 /* sigh -- "0" often padded with blanks */
1379 xlabel= strtok(xlabel)(1,);
1380 if (mpwr)
1381 xlabel+= swrite(format="%se%+02ld",(ndecimals?"":"."),mpwr);
1382
1383 return 0;
1384 }
1385
1386 /* ------------------------------------------------------------------------ */
1387