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;
 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   XA( XFillPolygon(display,pixmap, dmcurrent->fillgc,
 145                    points,3,Convex,CoordModeOrigin) );
 146   XA( XDrawLines(display,pixmap, dmcurrent->linegc,
 147                  points, 4,CoordModeOrigin) );
 148 }
 149
 150 static const unsigned long core_event_mask=
 151   ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask|
 152   KeyPressMask;
 153
 154 static void mkpixmaps(void) {
 155   for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
 156     XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,x11depth) );
 157     doublebuffers[currentbuffer]= pixmap;
 158   }
 159   currentbuffer= 0;
 160 }
 161
 162 static void mkgcs(DrawingMode *dm, unsigned long planes) {
 163   XGCValues gcv;
 164
 165   gcv.function= GXcopy;
 166   gcv.foreground= WhitePixel(display,x11screen);
 167   gcv.plane_mask= planes;
 168   dm->linegc= XCreateGC(display,pixmap,
 169                         GCFunction|GCForeground|GCPlaneMask,
 170                         &gcv);
 171
 172   gcv.function= GXclear;
 173   dm->fillgc= XCreateGC(display,pixmap,
 174                         GCFunction|GCPlaneMask,
 175                         &gcv);
 176 }
 177
 178 static void display_prepare(void) {
 179   XSetWindowAttributes wa;
 180   XSizeHints hints;
 181
 182   XA( display= XOpenDisplay(0) );
 183   x11screen= DefaultScreen(display);
 184   x11depth= DefaultDepth(display,x11screen);
 185   XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&visinfo) );
 186
 187   wa.event_mask= core_event_mask;
 188   XA( window= XCreateWindow(display, DefaultRootWindow(display),
 189                             0,0, wwidth,wheight, 0,x11depth,
 190                             InputOutput, visinfo.visual,
 191                             CWEventMask, &wa) );
 192
 193   hints.flags= USPosition;
 194   hints.x= 10;
 195   hints.y= 10;
 196   XSetWMNormalHints(display,window,&hints);
 197
 198   mkpixmaps();
 199
 200   mkgcs(&dmwhite, AllPlanes);
 201   mkgcs(&dmblue, visinfo.blue_mask);
 202   mkgcs(&dmred, visinfo.red_mask);
 203 }
 204
 205 static void drawtriangles(const DrawingMode *dm) {
 206   Triangle *const *t;
 207   int i;
 208
 209   dmcurrent= dm;
 210   for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
 211     drawtriangle(*t);
 212 }
 213
 214 static void display_conformation(void) {
 215   pixmap= doublebuffers[currentbuffer];
 216
 217   XA( XFillRectangle(display,pixmap,dmwhite.fillgc,0,0,wwidth,wheight) );
 218
 219   if (eyes_apart > 0) {
 220     const double stationary= 0.07;
 221
 222     eye_x= eyes_apart < eyes_apart_preferred
 223               ? eyes_apart :
 224            eyes_apart < (eyes_apart_preferred + stationary)
 225               ? eyes_apart_preferred
 226               : eyes_apart - stationary;
 227     eye_x /= sizeadj_scale;
 228     drawtriangles(&dmblue);
 229     eye_x= -eye_x;
 230     drawtriangles(&dmred);
 231   } else {
 232     drawtriangles(&dmwhite);
 233     printf("shown, %d/%d triangles cut\n", ncut, ntris);
 234   }
 235
 236   XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
 237   XA( XClearWindow(display,window) );
 238   currentbuffer= !currentbuffer;
 239 }
 240
 241 static void show(void) {
 242   scale_wmindim= sizeadj_scale * wmindim;
 243   read_input();
 244   transform_coordinates();
 245   generate_display_list();
 246   sort_display_list();
 247   display_conformation();
 248 }
 249
 250 typedef struct {
 251   const char *name;
 252   void (*start)(void);
 253   void (*delta)(double dx, double dy);
 254   void (*conclude)(void);
 255   void (*abandon)(void);
 256 } Drag;
 257
 258 #define DRAG(x)                                 \
 259   static const Drag drag_##x= {                 \
 260     #x, drag_##x##_start, drag_##x##_delta,     \
 261     drag_##x##_conclude, drag_##x##_abandon     \
 262   }
 263
 264 #define DRAG_SAVING(x, thing)                           \
 265   static typeof(thing) original_##thing;                \
 266   static void drag_##x##_start(void) {                  \
 267     memcpy(&original_##thing, &thing, sizeof(thing));   \
 268   }                                                     \
 269   static void drag_##x##_conclude(void) { }             \
 270   static void drag_##x##_abandon(void) {                \
 271     memcpy(&thing, &original_##thing, sizeof(thing));   \
 272     show();                                             \
 273   }                                                     \
 274   DRAG(x)
 275
 276 static void drag_none_start(void) { }
 277 static void drag_none_delta(double dx, double dy) { }
 278 static void drag_none_conclude(void) { }
 279 static void drag_none_abandon(void) { }
 280 DRAG(none);
 281
 282 static void pvectorcore(const char *n, double v[D3]) {
 283   int k;
 284   printf("%10s [ ",n);
 285   K printf("%# 10.10f ",v[k]);
 286   printf("]\n");
 287 }
 288 static void pvector(const char *n, double v[D3]) {
 289   pvectorcore(n,v);
 290   putchar('\n');
 291 }
 292 static void pmatrix(const char *n, double m[D3][D3]) {
 293   int j;
 294   for (j=0; j<D3; j++) { pvectorcore(n,m[j]); n=""; }
 295   putchar('\n');
 296 }
 297 #define PMATRIX(x) pmatrix(#x,x);
 298
 299 static void drag_rotate_delta(double dx, double dy) {
 300   /* We multiple our transformation matrix by a matrix:
 301    *
 302    * If we just had y movement, we would rotate about x axis:
 303    *  rotation X = [  1    0   0 ]
 304    *               [  0   cy  sy ]
 305    *               [  0  -sy  cy ]
 306    *  where cy and sy are sin and cos of y rotation
 307    *
 308    * But we should pre-rotate this by a rotation about the z axis
 309    * to get it to the right angle (to include x rotation).  So
 310    * we make cy and sy be cos() and sin(hypot(x,y)) and use
 311    * with cr,sr as cos() and sin(atan2(y,y)):
 312    *
 313    * Ie we would do  T' = R^T X R T   where
 314    *             or  T' =    C    T   where  C = R^T X R  and
 315    *
 316    *  adjustment R = [  cr  sr  0 ]
 317    *                 [ -sr  cr  0 ]
 318    *                 [  0    0  1 ]
 319    */
 320
 321   double rotx[D3][D3], adjr[D3][D3];
 322   GSL_MATRIX(rotx);
 323   GSL_MATRIX(adjr);
 324
 325   static double temp[D3][D3], change[D3][D3];
 326   static GSL_MATRIX(temp);
 327   static GSL_MATRIX(change);
 328
 329   double d= hypot(dx,dy);
 330   if (d < 1e-6) return;
 331
 332   double ang= d*2.0;
 333
 334   double cy= cos(ang);
 335   double sy= sin(ang);
 336   double cr= -dy / d;
 337   double sr=  dx / d;
 338   printf("\n d=%g cy,sy=%g,%g cr,sr=%g,%g\n\n", d,cy,sy,cr,sr);
 339
 340   rotx[0][0]=   1;    rotx[0][1]=   0;     rotx[0][2]=  0;
 341   rotx[1][0]=   0;    rotx[1][1]=  cy;     rotx[1][2]= sy;
 342   rotx[2][0]=   0;    rotx[2][1]= -sy;     rotx[2][2]= cy;
 343   PMATRIX(rotx);
 344
 345   adjr[0][0]=  cr;    adjr[0][1]=  sr;     adjr[0][2]=  0;
 346   adjr[1][0]= -sr;    adjr[1][1]=  cr;     adjr[1][2]=  0;
 347   adjr[2][0]=   0;    adjr[2][1]=   0;     adjr[2][2]=  1;
 348   PMATRIX(adjr);
 349
 350   GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
 351                      &rotx_gsl,&adjr_gsl,
 352                      0.0, &temp_gsl) );
 353   PMATRIX(temp);
 354
 355   GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
 356                      &adjr_gsl,&temp_gsl,
 357                      0.0, &change_gsl) );
 358   PMATRIX(change);
 359
 360   static double skew[D3][D3];
 361   static GSL_MATRIX(skew);
 362
 363   GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
 364                      &change_gsl,&transform_gsl,
 365                      0.0, &skew_gsl) );
 366   PMATRIX(skew);
 367
 368   memcpy(&transform,&skew,sizeof(transform));
 369   show();
 370   return;
 371
 372   /* Now we want to normalise skew, the result becomes new transform */
 373   double svd_v[D3][D3];
 374   GSL_MATRIX(svd_v);
 375
 376   double sigma[D3], tau[D3];
 377   GSL_VECTOR(sigma);
 378   GSL_VECTOR(tau);
 379
 380   /* We use notation from Wikipedia Polar_decomposition
 381    *     Wikipedia's  W      is GSL's U
 382    *     Wikipedia's  Sigma  is GSL's S
 383    *     Wikipedia's  V      is GSL's V
 384    *     Wikipedia's  U      is our desired result
 385    * Wikipedia which says if the SVD is    A = W Sigma V*
 386    * then the polar decomposition is       A = U P
 387    *   where                               P = V Sigma V*
 388    *   and                                 U = W V*
 389    */
 390
 391   GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
 392   pmatrix("W",    skew);
 393   pvector("Sigma",sigma);
 394   pmatrix("V",    svd_v);
 395
 396   /* We only need U, not P. */
 397   GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
 398                      &skew_gsl,&svd_v_gsl,
 399                      0.0,&transform_gsl) );
 400
 401   pmatrix("U", transform);
 402
 403   printf("drag_rotate_delta...\n");
 404   show();
 405 }
 406 DRAG_SAVING(rotate, transform);
 407
 408 static void drag_sizeadj_delta(double dx, double dy) {
 409   sizeadj_scale *= pow(3.0, -dy);
 410   show();
 411 }
 412 DRAG_SAVING(sizeadj, sizeadj_scale);
 413
 414 static void drag_3d_delta(double dx, double dy) {
 415   eyes_apart += dx * 0.1;
 416   if (eyes_apart < eyes_apart_min) eyes_apart= eyes_apart_min;
 417   printf("sizeadj eyes_apart %g\n", eyes_apart);
 418   show();
 419 }
 420 DRAG_SAVING(3d, eyes_apart);
 421
 422 static const Drag *drag= &drag_none;
 423
 424 static int drag_last_x, drag_last_y;
 425
 426 static void drag_position(int x, int y) {
 427   drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
 428               (y - drag_last_y) * 1.0 / wmaxdim);
 429   drag_last_x= x;
 430   drag_last_y= y;
 431 }
 432
 433 static void event_button(XButtonEvent *e) {
 434   if (e->window != window || !e->same_screen) return;
 435   if (e->type == ButtonPress) {
 436     if (e->state || drag != &drag_none) {
 437       printf("drag=%s press state=0x%lx abandon\n",
 438              drag->name, (unsigned long)e->state);
 439       drag->abandon();
 440       drag= &drag_none;
 441       return;
 442     }
 443     switch (e->button) {
 444     case Button1: drag= &drag_rotate;  break;
 445     case Button2: drag= &drag_sizeadj; break;
 446     case Button3: drag= &drag_3d;      break;
 447     default: printf("unknown drag start %d\n", e->button);
 448     }
 449     printf("drag=%s press button=%lu start %d,%d\n",
 450            drag->name, (unsigned long)e->button, e->x, e->y);
 451     drag_last_x= e->x;
 452     drag_last_y= e->y;
 453     drag->start();
 454   }
 455   if (e->type == ButtonRelease) {
 456     printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
 457     drag_position(e->x, e->y);
 458     drag->conclude();
 459     drag= &drag_none;
 460   }
 461 }
 462
 463 static void event_motion(int x, int y) {
 464   printf("drag=%s motion %d,%d\n", drag->name, x, y);
 465   drag_position(x,y);
 466 }
 467
 468 static void event_key(XKeyEvent *e) {
 469   KeySym ks;
 470   char buf[10];
 471   int r;
 472
 473   r= XLookupString(e,buf,sizeof(buf)-1,&ks,0);
 474   if (!r) {
 475     printf("XLookupString keycode=%u state=0x%x gave %d\n",
 476            e->keycode, e->state, r);
 477     return;
 478   }
 479
 480   if (!strcmp(buf,"q")) exit(0);
 481   if (!strcmp(buf,"d")) {
 482     eyes_apart= eyes_apart>0 ? eyes_apart_min : eyes_apart_preferred;
 483     show();
 484     return;
 485   }
 486 }
 487
 488 static void event_config(XConfigureEvent *e) {
 489   if (e->width == wwidth && e->height == wheight)
 490     return;
 491
 492   wwidth= e->width;  wheight= e->height;
 493   wmaxdim= wwidth > wheight ? wwidth : wheight;
 494   wmindim= wwidth < wheight ? wwidth : wheight;
 495
 496   XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
 497   for (currentbuffer=0; currentbuffer<2; currentbuffer++)
 498     XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
 499
 500   mkpixmaps();
 501   show();
 502 }
 503
 504 static void check_input(void) {
 505   struct stat newstab;
 506   int r;
 507
 508   r= stat(input_filename, &newstab);
 509   if (r<0) diee("could not check input");
 510
 511 #define CI(x) if (newstab.st_##x == input_stab.st_##x) ; else goto changed
 512   CI(dev);
 513   CI(ino);
 514   CI(size);
 515   CI(mtime);
 516 #undef CI
 517   return;
 518
 519  changed:
 520   show();
 521 }
 522
 523 int main(int argc, const char *const *argv) {
 524   static const int wantedevents= POLLIN|POLLPRI|POLLERR|POLLHUP;
 525
 526   XEvent event;
 527   int k, i, r, *xfds, nxfds, polls_alloc=0;
 528   struct pollfd *polls=0;
 529   int motion_deferred=0, motion_x=-1, motion_y=-1;
 530
 531   if (argc != 2 || argv[1][0]=='-') {
 532     fputs("need filename\n",stderr); exit(8);
 533   }
 534   input_filename= argv[1];
 535
 536   read_input();
 537   K transform[k][k]= 1.0;
 538   display_prepare();
 539   show();
 540
 541   XMapWindow(display,window);
 542   for (;;) {
 543
 544     XA( XInternalConnectionNumbers(display, &xfds, &nxfds) );
 545     if (polls_alloc <= nxfds) {
 546       polls_alloc= nxfds + polls_alloc + 1;
 547       polls= realloc(polls, sizeof(*polls) * polls_alloc);
 548       if (!polls) diee("realloc for pollfds");
 549     }
 550     for (i=0; i<nxfds; i++) {
 551       polls[i].fd= xfds[i];
 552       polls[i].events= wantedevents;
 553       polls[i].revents= 0;
 554     }
 555     XFree(xfds);
 556
 557     polls[i].fd= ConnectionNumber(display);
 558     polls[i].events= wantedevents;
 559
 560     r= poll(polls, nxfds+1, motion_deferred ? 0 : 200);
 561     if (r<0) diee("poll");
 562
 563     for (i=0; i<nxfds; i++)
 564       if (polls[i].revents)
 565         XProcessInternalConnection(display, polls[i].fd);
 566
 567     r= XCheckMaskEvent(display,~0UL,&event);
 568     if (!r) {
 569       if (motion_deferred) {
 570         event_motion(motion_x, motion_y);
 571         motion_deferred=0;
 572       }
 573       check_input();
 574       continue;
 575     }
 576
 577     switch (event.type) {
 578
 579     case ButtonPress:
 580     case ButtonRelease:     event_button(&event.xbutton);     break;
 581
 582     case KeyPress:          event_key(&event.xkey);           break;
 583
 584     case ConfigureNotify:   event_config(&event.xconfigure);  break;
 585
 586     case MotionNotify:
 587       motion_x= event.xmotion.x;
 588       motion_y= event.xmotion.y;
 589       motion_deferred= 1;
 590       continue;
 591
 592     default:
 593       printf("unknown event type %u 0x%x\n", event.type,event.type);
 594     }
 595   }
 596 }
 597