chiark / gitweb /
allow turn off hidden line removal
[moebius2.git] / view.c
1 /*
2  * Displays a conformation
3  */
4
5 #include <X11/Xlib.h>
6 #include <X11/Xutil.h>
7
8 #include <sys/types.h>
9 #include <sys/stat.h>
10 #include <sys/poll.h>
11
12 #include "mgraph.h"
13
14 #define MAXTRIS (N*2)
15
16 typedef struct { double vertex[3][D3]; } Triangle;
17
18 static Triangle trisbuffer[MAXTRIS], *displaylist[MAXTRIS];
19 static int ntris;
20 static Vertices conformation;
21
22 static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
23 static GSL_MATRIX(transform);
24
25 static FILE *input_f;
26 static struct stat input_stab;
27 static const char *input_filename;
28
29 static void read_input(void) {
30   int r;
31   
32   if (input_f) fclose(input_f);
33   input_f= fopen(input_filename, "rb");  if (!input_f) diee("input file");
34
35   if (fstat(fileno(input_f), &input_stab)) diee("fstat input file");
36
37   errno= 0;
38   r= fread(&conformation,sizeof(conformation),1,input_f);
39   if (r!=1) diee("fread");
40 }
41
42 static void transform_coordinates(void) {
43   double result[D3];
44   GSL_VECTOR(result);
45   gsl_vector input_gsl= { D3,1 };
46
47   int v, k;
48   
49   FOR_VERTEX(v) {
50     input_gsl.data= &conformation[v][0];
51     GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
52                        0.0, &result_gsl) );
53     K conformation[v][k]= result[k];
54   }
55 }
56
57 static void addtriangle(int va, int vb, int vc) {
58   Triangle *t= &trisbuffer[ntris];
59   int k;
60   
61   assert(ntris < MAXTRIS);
62   K {
63     t->vertex[0][k]= conformation[va][k];
64     t->vertex[1][k]= conformation[vb][k];
65     t->vertex[2][k]= conformation[vc][k];
66   }
67   displaylist[ntris++]= t;
68 }
69
70 static void generate_display_list(void) {
71   int vb, ve[3], e;
72
73   ntris= 0;
74   FOR_VERTEX(vb) {
75     /* We use the two triangles in the parallelogram vb, vb+e0, vb+e1, vb+e2.
76      * We go round each triangle clockwise (although our surface is non-
77      * orientable so it shouldn't matter).
78      */
79     for (e=0; e<3; e++) ve[e]= EDGE_END2(vb,e);
80     if (ve[1]>=0) {
81       if (ve[0]>=0) addtriangle(vb,ve[0],ve[1]);
82       if (ve[2]>=0) addtriangle(vb,ve[1],ve[2]);
83     }
84   }
85 }    
86
87 static int dl_compare(const void *tav, const void *tbv) {
88   int i;
89   const Triangle *const *tap= tav, *ta= *tap;
90   const Triangle *const *tbp= tbp, *tb= *tbp;
91   double za=0, zb=0;
92   for (i=0; i<3; i++) {
93     za += ta->vertex[i][2];
94     zb += tb->vertex[i][2];
95   }
96   return za > zb ? -1 :
97          za < zb ? +1 : 0;
98 }
99
100 static void sort_display_list(void) {
101   qsort(displaylist, ntris, sizeof(*displaylist), dl_compare);
102 }
103
104 /*---------- X stuff ----------*/
105
106 #define WSZ 400
107
108 typedef struct { GC fillgc, linegc; } DrawingMode;
109
110 static Display *display;
111 static Pixmap pixmap, doublebuffers[2];
112 static Window window;
113
114 static DrawingMode dmred, dmblue, dmwhite;
115 static const DrawingMode *dmcurrent;
116 static int wwidth=WSZ, wheight=WSZ, wmindim=WSZ, wmaxdim=WSZ;
117 static int ncut, currentbuffer, x11depth, x11screen, wireframe;
118 XVisualInfo visinfo;
119
120 static double sizeadj_scale= 0.3, eyes_apart, scale_wmindim;
121 static double eye_z= -10, eye_x=0;
122 static double cut_z= -9;
123 static const double eyes_apart_preferred=0.05, eyes_apart_min= -0.02;
124
125
126 static void drawtriangle(const Triangle *t) {
127   XPoint points[4];
128   int i;
129
130   for (i=0; i<3; i++) {
131     double *v= t->vertex[i];
132     double x= v[0];
133     double y= v[1];
134     double z= v[2];
135
136     if (z < cut_z) { ncut++; return; }
137     
138     double zezezp= eye_z / (eye_z - z);
139     points[i].x= scale_wmindim * (zezezp * (x - eye_x) + eye_x) + wwidth/2;
140     points[i].y= scale_wmindim * (zezezp *  y                 ) + wheight/2;
141   }
142   points[3]= points[0];
143
144   if (!wireframe)
145     XA( XFillPolygon(display,pixmap, dmcurrent->fillgc,
146                      points,3,Convex,CoordModeOrigin) );
147   XA( XDrawLines(display,pixmap, dmcurrent->linegc,
148                  points, 4,CoordModeOrigin) );
149 }
150
151 static const unsigned long core_event_mask=
152   ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask|
153   KeyPressMask;
154
155 static void mkpixmaps(void) {
156   for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
157     XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,x11depth) );
158     doublebuffers[currentbuffer]= pixmap;
159   }
160   currentbuffer= 0;
161 }
162
163 static void mkgcs(DrawingMode *dm, unsigned long planes) {
164   XGCValues gcv;
165
166   gcv.function= GXcopy;
167   gcv.foreground= WhitePixel(display,x11screen);
168   gcv.plane_mask= planes;
169   dm->linegc= XCreateGC(display,pixmap,
170                         GCFunction|GCForeground|GCPlaneMask,
171                         &gcv);
172
173   gcv.function= GXclear;
174   dm->fillgc= XCreateGC(display,pixmap,
175                         GCFunction|GCPlaneMask,
176                         &gcv);
177 }
178
179 static void display_prepare(void) {
180   XSetWindowAttributes wa;
181   XSizeHints hints;
182   
183   XA( display= XOpenDisplay(0) );
184   x11screen= DefaultScreen(display);
185   x11depth= DefaultDepth(display,x11screen);
186   XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&visinfo) );
187   
188   wa.event_mask= core_event_mask;
189   XA( window= XCreateWindow(display, DefaultRootWindow(display),
190                             0,0, wwidth,wheight, 0,x11depth,
191                             InputOutput, visinfo.visual,
192                             CWEventMask, &wa) );
193
194   hints.flags= USPosition;
195   hints.x= 10;
196   hints.y= 10;
197   XSetWMNormalHints(display,window,&hints);
198
199   mkpixmaps();
200
201   mkgcs(&dmwhite, AllPlanes);
202   mkgcs(&dmblue, visinfo.blue_mask);
203   mkgcs(&dmred, visinfo.red_mask);
204 }
205
206 static void drawtriangles(const DrawingMode *dm) {
207   Triangle *const *t;
208   int i;
209
210   dmcurrent= dm;
211   for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
212     drawtriangle(*t);
213 }
214
215 static void display_conformation(void) {
216   pixmap= doublebuffers[currentbuffer];
217
218   XA( XFillRectangle(display,pixmap,dmwhite.fillgc,0,0,wwidth,wheight) );
219
220   if (eyes_apart > 0) {
221     const double stationary= 0.07;
222
223     eye_x= eyes_apart < eyes_apart_preferred
224               ? eyes_apart :
225            eyes_apart < (eyes_apart_preferred + stationary)
226               ? eyes_apart_preferred
227               : eyes_apart - stationary;
228     eye_x /= sizeadj_scale;
229     drawtriangles(&dmblue);
230     eye_x= -eye_x;
231     drawtriangles(&dmred);
232   } else {
233     drawtriangles(&dmwhite);
234     printf("shown, %d/%d triangles cut\n", ncut, ntris);
235   }
236   
237   XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
238   XA( XClearWindow(display,window) );
239   currentbuffer= !currentbuffer;
240 }
241
242 static void show(void) {
243   scale_wmindim= sizeadj_scale * wmindim;
244   read_input();
245   transform_coordinates();
246   generate_display_list();
247   sort_display_list();
248   display_conformation();
249 }
250
251 typedef struct {
252   const char *name;
253   void (*start)(const XButtonEvent *e);
254   void (*delta)(double dx, double dy);
255   void (*conclude)(void);
256   void (*abandon)(void);
257 } Drag;
258
259 #define DRAG(x)                                 \
260   static const Drag drag_##x= {                 \
261     #x, drag_##x##_start, drag_##x##_delta,     \
262     drag_##x##_conclude, drag_##x##_abandon     \
263   }
264
265 #define DRAG_SAVING(x, thing, hook)                     \
266   static typeof(thing) original_##thing;                \
267   static void drag_##x##_start(const XButtonEvent *e) { \
268     memcpy(&original_##thing, &thing, sizeof(thing));   \
269     hook;                                               \
270   }                                                     \
271   static void drag_##x##_conclude(void) { }             \
272   static void drag_##x##_abandon(void) {                \
273     memcpy(&thing, &original_##thing, sizeof(thing));   \
274     show();                                             \
275   }                                                     \
276   DRAG(x)
277
278 static void drag_none_start(const XButtonEvent *e) { }
279 static void drag_none_delta(double dx, double dy) { }
280 static void drag_none_conclude(void) { }
281 static void drag_none_abandon(void) { }
282 DRAG(none);
283
284 static void pvectorcore(const char *n, double v[D3]) {
285   int k;
286   printf("%10s [ ",n);
287   K printf("%# 10.10f ",v[k]);
288   printf("]\n");
289 }
290 static void pvector(const char *n, double v[D3]) {
291   pvectorcore(n,v);
292   putchar('\n');
293 }
294 static void pmatrix(const char *n, double m[D3][D3]) {
295   int j;
296   for (j=0; j<D3; j++) { pvectorcore(n,m[j]); n=""; }
297   putchar('\n');
298 }
299 #define PMATRIX(x) pmatrix(#x,x);
300
301 static double drag_transform_conv_x_z= 0;
302 static double drag_transform_conv_y_z= 0;
303
304 static void drag_transform_prep(const XButtonEvent *e) {
305   static const double factor= 2.5;
306   drag_transform_conv_x_z= MAX( MIN(e->y * factor / wheight - (factor/2),
307                                     1.0), -1.0);
308   drag_transform_conv_y_z= MAX( MIN(e->x * factor / wwidth  - (factor/2),
309                                     1.0), -1.0);
310   printf("drag_transform_conv_{x,y}_z = %g,%g\n",
311          drag_transform_conv_x_z, drag_transform_conv_y_z);
312 }
313
314 static void make_z_rotation(double rotz[D3][D3], double cz, double sz) {
315   rotz[0][0]=  cz;    rotz[0][1]=  sz;     rotz[0][2]=  0;
316   rotz[1][0]= -sz;    rotz[1][1]=  cz;     rotz[1][2]=  0;
317   rotz[2][0]=   0;    rotz[2][1]=   0;     rotz[2][2]=  1;
318 }  
319
320 static void drag_rotate_delta(double dx, double dy) {
321   /* We multiple our transformation matrix by a matrix:
322    *
323    * If we just had y movement, we would rotate about x axis:
324    *  rotation X = [  1    0   0 ]
325    *               [  0   cy  sy ]
326    *               [  0  -sy  cy ]
327    *  where cy and sy are sin and cos of y rotation
328    *
329    * But we should pre-rotate this by a rotation about the z axis
330    * to get it to the right angle (to include x rotation).  So
331    * we make cy and sy be cos() and sin(hypot(x,y)) and use
332    * with cr,sr as cos() and sin(atan2(y,y)):
333    *
334    * Ie we would do  T' = Z R^T X R T   where
335    *             or  T' = Z    C    T   where  C = R^T X R  and
336    *
337    *  adjustment R = [  cr  sr  0 ]
338    *                 [ -sr  cr  0 ]
339    *                 [  0    0  1 ]    or make_z_rotation(cr,sr)
340    *
341    *  rotation Z   = [  cz  sz  0 ]
342    *                 [ -sz  cz  0 ]
343    *                 [  0    0  1 ]    or make_z_rotation(cz,sz)
344    */
345
346   double rotx[D3][D3], adjr[D3][D3], rotz[D3][D3];
347   GSL_MATRIX(rotx);
348   GSL_MATRIX(adjr);
349   GSL_MATRIX(rotz);
350
351   static double temp[D3][D3], change[D3][D3];
352   static GSL_MATRIX(temp);
353   static GSL_MATRIX(change);
354
355   printf("\nTRANSFORM %g, %g\n", dx,dy);
356
357   double dz= -drag_transform_conv_x_z * dx +
358               drag_transform_conv_y_z * dy;
359           
360   dx *= (1 - fabs(drag_transform_conv_x_z));
361   dy *= (1 - fabs(drag_transform_conv_y_z));
362
363   double d= hypot(dx,dy);
364
365   printf(" dx,dy,dz = %g, %g, %g    d = %g\n", dx,dy,dz, d);
366
367   if (hypot(d,dz) < 1e-6) return;
368
369   printf(" no xy rotation\n");
370
371   double ang= d*2.0;
372   
373   double cy= cos(ang);
374   double sy= sin(ang);
375
376   double cr, sr;
377   if (d > 1e-6) {
378     cr= -dy / d;
379     sr=  dx / d;
380   } else {
381     cr= 1;
382     sr= 0;
383   }
384   printf("\n d=%g cy,sy=%g,%g cr,sr=%g,%g\n\n", d,cy,sy,cr,sr);
385   
386   rotx[0][0]=   1;    rotx[0][1]=   0;     rotx[0][2]=  0;
387   rotx[1][0]=   0;    rotx[1][1]=  cy;     rotx[1][2]= sy;
388   rotx[2][0]=   0;    rotx[2][1]= -sy;     rotx[2][2]= cy;
389   PMATRIX(rotx);
390
391   make_z_rotation(adjr,cr,sr);
392   PMATRIX(adjr);
393
394   GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
395                      &rotx_gsl,&adjr_gsl,
396                      0.0, &temp_gsl) );
397   pmatrix("X R", temp);
398   
399   GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
400                      &adjr_gsl,&temp_gsl,
401                      0.0, &change_gsl) );
402   PMATRIX(change);
403
404   double angz= dz*2.0;
405   make_z_rotation(rotz,cos(angz),sin(angz));
406   PMATRIX(rotz);
407
408   GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
409                      &rotz_gsl,&change_gsl,
410                      0.0, &temp_gsl) );
411   pmatrix("Z C", temp);
412
413   static double skew[D3][D3];
414   static GSL_MATRIX(skew);
415   
416   GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
417                      &temp_gsl, &transform_gsl,
418                      0.0, &skew_gsl) );
419   PMATRIX(skew);
420
421   memcpy(&transform,&skew,sizeof(transform));
422   show();
423   return;
424
425   /* Now we want to normalise skew, the result becomes new transform */
426   double svd_v[D3][D3];
427   GSL_MATRIX(svd_v);
428
429   double sigma[D3], tau[D3];
430   GSL_VECTOR(sigma);
431   GSL_VECTOR(tau);
432   
433   /* We use notation from Wikipedia Polar_decomposition
434    *     Wikipedia's  W      is GSL's U
435    *     Wikipedia's  Sigma  is GSL's S
436    *     Wikipedia's  V      is GSL's V
437    *     Wikipedia's  U      is our desired result
438    * Wikipedia which says if the SVD is    A = W Sigma V*
439    * then the polar decomposition is       A = U P
440    *   where                               P = V Sigma V*
441    *   and                                 U = W V*
442    */
443   
444   GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
445   pmatrix("W",    skew);
446   pvector("Sigma",sigma);
447   pmatrix("V",    svd_v);
448
449   /* We only need U, not P. */
450   GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
451                      &skew_gsl,&svd_v_gsl,
452                      0.0,&transform_gsl) );
453
454   pmatrix("U", transform);
455
456   printf("drag_rotate_delta...\n");
457   show();
458 }
459 DRAG_SAVING(rotate, transform, drag_transform_prep(e));
460
461 static void drag_sizeadj_delta(double dx, double dy) {
462   sizeadj_scale *= pow(3.0, -dy);
463   show();
464 }
465 DRAG_SAVING(sizeadj, sizeadj_scale, );
466
467 static void drag_3d_delta(double dx, double dy) {
468   eyes_apart += dx * 0.1;
469   if (eyes_apart < eyes_apart_min) eyes_apart= eyes_apart_min;
470   printf("sizeadj eyes_apart %g\n", eyes_apart);
471   show();
472 }
473 DRAG_SAVING(3d, eyes_apart, );
474
475 static const Drag *drag= &drag_none;
476
477 static int drag_last_x, drag_last_y;
478
479 static void drag_position(int x, int y) {
480   drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
481               (y - drag_last_y) * 1.0 / wmaxdim);
482   drag_last_x= x;
483   drag_last_y= y;
484 }
485
486 static void event_button(XButtonEvent *e) {
487   if (e->window != window || !e->same_screen) return;
488   if (e->type == ButtonPress) {
489     if (e->state || drag != &drag_none) {
490       printf("drag=%s press state=0x%lx abandon\n",
491              drag->name, (unsigned long)e->state);
492       drag->abandon();
493       drag= &drag_none;
494       return;
495     }
496     switch (e->button) {
497     case Button1: drag= &drag_rotate;  break;
498     case Button2: drag= &drag_sizeadj; break;
499     case Button3: drag= &drag_3d;      break;
500     default: printf("unknown drag start %d\n", e->button);
501     }
502     printf("drag=%s press button=%lu start %d,%d\n",
503            drag->name, (unsigned long)e->button, e->x, e->y);
504     drag_last_x= e->x;
505     drag_last_y= e->y;
506     drag->start(e);
507   }
508   if (e->type == ButtonRelease) {
509     printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
510     drag_position(e->x, e->y);
511     drag->conclude();
512     drag= &drag_none;
513   }
514 }
515
516 static void event_motion(int x, int y) {
517   printf("drag=%s motion %d,%d\n", drag->name, x, y);
518   drag_position(x,y);
519 }
520
521 static void event_key(XKeyEvent *e) {
522   KeySym ks;
523   char buf[10];
524   int r;
525
526   r= XLookupString(e,buf,sizeof(buf)-1,&ks,0);
527   if (!r) {
528     printf("XLookupString keycode=%u state=0x%x gave %d\n",
529            e->keycode, e->state, r);
530     return;
531   }
532
533   if (!strcmp(buf,"q")) exit(0);
534   if (!strcmp(buf,"w")) {
535     wireframe= !wireframe;
536     show();
537     return;
538   }
539   if (!strcmp(buf,"d")) {
540     eyes_apart= eyes_apart>0 ? eyes_apart_min : eyes_apart_preferred;
541     show();
542     return;
543   }
544 }
545
546 static void event_config(XConfigureEvent *e) {
547   if (e->width == wwidth && e->height == wheight)
548     return;
549
550   wwidth= e->width;  wheight= e->height;
551   wmaxdim= wwidth > wheight ? wwidth : wheight;
552   wmindim= wwidth < wheight ? wwidth : wheight;
553   
554   XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
555   for (currentbuffer=0; currentbuffer<2; currentbuffer++)
556     XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
557
558   mkpixmaps();
559   show();
560 }
561
562 static void check_input(void) {
563   struct stat newstab;
564   int r;
565
566   r= stat(input_filename, &newstab);
567   if (r<0) diee("could not check input");
568
569 #define CI(x) if (newstab.st_##x == input_stab.st_##x) ; else goto changed
570   CI(dev);
571   CI(ino);
572   CI(size);
573   CI(mtime);
574 #undef CI
575   return;
576
577  changed:
578   show();
579 }
580
581 int main(int argc, const char *const *argv) {
582   static const int wantedevents= POLLIN|POLLPRI|POLLERR|POLLHUP;
583
584   XEvent event;
585   int k, i, r, *xfds, nxfds, polls_alloc=0;
586   struct pollfd *polls=0;
587   int motion_deferred=0, motion_x=-1, motion_y=-1;
588   
589   if (argc != 2 || argv[1][0]=='-') {
590     fputs("need filename\n",stderr); exit(8);
591   }
592   input_filename= argv[1];
593
594   read_input();
595   K transform[k][k]= 1.0;
596   display_prepare();
597   show();
598
599   XMapWindow(display,window);
600   for (;;) {
601
602     XA( XInternalConnectionNumbers(display, &xfds, &nxfds) );
603     if (polls_alloc <= nxfds) {
604       polls_alloc= nxfds + polls_alloc + 1;
605       polls= realloc(polls, sizeof(*polls) * polls_alloc);
606       if (!polls) diee("realloc for pollfds");
607     }
608     for (i=0; i<nxfds; i++) {
609       polls[i].fd= xfds[i];
610       polls[i].events= wantedevents;
611       polls[i].revents= 0;
612     }
613     XFree(xfds);
614
615     polls[i].fd= ConnectionNumber(display);
616     polls[i].events= wantedevents;
617
618     r= poll(polls, nxfds+1, motion_deferred ? 0 : 200);
619     if (r<0) diee("poll");
620
621     for (i=0; i<nxfds; i++)
622       if (polls[i].revents)
623         XProcessInternalConnection(display, polls[i].fd);
624
625     r= XCheckMaskEvent(display,~0UL,&event);
626     if (!r) {
627       if (motion_deferred) {
628         event_motion(motion_x, motion_y);
629         motion_deferred=0;
630       }
631       check_input();
632       continue;
633     }
634     
635     switch (event.type) {
636
637     case ButtonPress:
638     case ButtonRelease:     event_button(&event.xbutton);     break;
639       
640     case KeyPress:          event_key(&event.xkey);           break;
641
642     case ConfigureNotify:   event_config(&event.xconfigure);  break;
643
644     case MotionNotify:
645       motion_x= event.xmotion.x;
646       motion_y= event.xmotion.y;
647       motion_deferred= 1;
648       continue;
649       
650     default:
651       printf("unknown event type %u 0x%x\n", event.type,event.type);
652     }
653   }
654 }
655