view.c

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