project.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+e1, vb+e0, 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[0]>=0) {
  72       if (ve[1]>=0) addtriangle(vb,ve[0],ve[1]);
  73       if (ve[2]>=0) addtriangle(vb,ve[2],ve[0]);
  74     }
  75   }
  76 }
  77
  78 static int dl_compare(const void *tav, const void *tbv) {
  79   const Triangle *const *tap= tav, *ta= *tap;
  80   const Triangle *const *tbp= tbp, *tb= *tbp;
  81   double za= ta->vertex[0][2];
  82   double zb= tb->vertex[0][2];
  83   return za > zb ? -1 :
  84          za < zb ? +1 : 0;
  85 }
  86
  87 static void sort_display_list(void) {
  88   qsort(displaylist, ntris, sizeof(*displaylist), dl_compare);
  89 }
  90
  91 /*---------- X stuff ----------*/
  92
  93 #define WSZ 400
  94
  95 static Display *display;
  96 static Pixmap pixmap, doublebuffers[2];
  97 static Window window;
  98 static GC linegc, fillgc;
  99 static int wwidth=WSZ, wheight=WSZ, wmindim=WSZ, wmaxdim=WSZ;
 100 static int ncut, currentbuffer, x11depth, x11screen;
 101
 102 static double scale= 0.3;
 103 static double eye_z= -10, eye_x= 0;
 104 static double cut_z= -9;
 105
 106 static void drawtriangle(const Triangle *t) {
 107   XPoint points[4];
 108   int i;
 109
 110   for (i=0; i<3; i++) {
 111     double *v= t->vertex[i];
 112     double x= v[0];
 113     double y= v[1];
 114     double z= v[2];
 115
 116     if (z < cut_z) { ncut++; return; }
 117
 118     double zezezp= eye_z / (eye_z - z);
 119     points[i].x= scale * wmindim * (zezezp * (x - eye_x) + eye_x) + wwidth/2;
 120     points[i].y= scale * wmindim * (zezezp *  y                 ) + wheight/2;
 121   }
 122   points[3]= points[0];
 123
 124   XA( XFillPolygon(display,pixmap,fillgc, points,3,Convex,CoordModeOrigin) );
 125   XA( XDrawLines(display,pixmap,linegc,points, 4,CoordModeOrigin) );
 126 }
 127
 128 static const unsigned long core_event_mask=
 129   ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask;
 130
 131 static void mkpixmaps(void) {
 132   for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
 133     XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,x11depth) );
 134     doublebuffers[currentbuffer]= pixmap;
 135   }
 136   currentbuffer= 0;
 137 }
 138
 139 static void display_prepare(void) {
 140   XGCValues gcv;
 141   XSetWindowAttributes wa;
 142   XVisualInfo vinfo;
 143   XSizeHints hints;
 144
 145   XA( display= XOpenDisplay(0) );
 146   x11screen= DefaultScreen(display);
 147   x11depth= DefaultDepth(display,x11screen);
 148   XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&vinfo) );
 149
 150   wa.event_mask= core_event_mask;
 151   XA( window= XCreateWindow(display, DefaultRootWindow(display),
 152                             0,0, wwidth,wheight, 0,x11depth,
 153                             InputOutput, vinfo.visual,
 154                             CWEventMask, &wa) );
 155
 156   hints.flags= USPosition;
 157   hints.x= 10;
 158   hints.y= 10;
 159   XSetWMNormalHints(display,window,&hints);
 160
 161   mkpixmaps();
 162
 163   gcv.function= GXcopy;
 164   gcv.plane_mask= ~0UL;
 165   gcv.foreground= WhitePixel(display,x11screen);
 166   linegc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask|GCForeground, &gcv);
 167
 168   gcv.function= GXclear;
 169   gcv.plane_mask= ~0UL;
 170   fillgc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask, &gcv);
 171 }
 172
 173 static void display_conformation(void) {
 174   Triangle *const *t;
 175   int i;
 176
 177   pixmap= doublebuffers[currentbuffer];
 178   XA( XFillRectangle(display,pixmap,fillgc,0,0,wwidth,wheight) );
 179   for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
 180     drawtriangle(*t);
 181   printf("shown, %d/%d triangles cut\n", ncut, ntris);
 182
 183   XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
 184   XA( XClearWindow(display,window) );
 185   currentbuffer= !currentbuffer;
 186 }
 187
 188 static void show(void) {
 189   read_input();
 190   transform_coordinates();
 191   generate_display_list();
 192   sort_display_list();
 193   display_conformation();
 194 }
 195
 196 typedef struct {
 197   const char *name;
 198   void (*start)(void);
 199   void (*delta)(double dx, double dy);
 200   void (*conclude)(void);
 201   void (*abandon)(void);
 202 } Drag;
 203
 204 #define DRAG(x)                                 \
 205   static const Drag drag_##x= {                 \
 206     #x, drag_##x##_start, drag_##x##_delta,     \
 207     drag_##x##_conclude, drag_##x##_abandon     \
 208   }
 209
 210 #define DRAG_SAVING(x, thing)                           \
 211   static typeof(thing) original_##thing;                \
 212   static void drag_##x##_start(void) {                  \
 213     memcpy(&original_##thing, &thing, sizeof(thing));   \
 214   }                                                     \
 215   static void drag_##x##_conclude(void) { }             \
 216   static void drag_##x##_abandon(void) {                \
 217     memcpy(&thing, &original_##thing, sizeof(thing));   \
 218     show();                                             \
 219   }                                                     \
 220   DRAG(x)
 221
 222 static void drag_none_start(void) { }
 223 static void drag_none_delta(double dx, double dy) { }
 224 static void drag_none_conclude(void) { }
 225 static void drag_none_abandon(void) { }
 226 DRAG(none);
 227
 228 static void pmatrix(const char *n, double m[D3][D3]) {
 229   int j,k;
 230   for (j=0; j<D3; j++) {
 231     printf("%10s [ ",n);
 232     K printf("%# f ",m[j][k]);
 233     printf("]\n");
 234     n="";
 235   }
 236 }
 237
 238 #define PMATRIX(x) pmatrix(#x,x);
 239
 240 static void drag_rotate_delta(double dx, double dy) {
 241   /* We multiple our transformation matrix by this matrix:
 242    *     [  1   0   0 ]
 243    *     [  0   1   0 ]
 244    *     [ dx  dy   1 ]
 245    * and then renormalise.
 246    */
 247
 248   static double rotateby[D3][D3]= {{1,0,0},{0,1,0},{-20,-20,1}};
 249   static GSL_MATRIX(rotateby);
 250
 251   double qr[D3][D3];
 252   gsl_matrix qr_gsl= { D3,D3,D3,&qr[0][0] };
 253
 254   double tau[D3];
 255   GSL_VECTOR(tau);
 256
 257   rotateby[2][0]= dx;
 258   rotateby[2][1]= dy;
 259
 260   PMATRIX(rotateby);
 261
 262   GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
 263                      &rotateby_gsl,&transform_gsl, 0.0,&qr_gsl) );
 264
 265   PMATRIX(rotateby);
 266
 267   GA( gsl_linalg_QR_decomp(&qr_gsl, &tau_gsl) );
 268   GA( gsl_linalg_QR_unpack(&qr_gsl, &tau_gsl,
 269                            &transform_gsl, &rotateby_gsl /*dummy*/) );
 270
 271   PMATRIX(transform);
 272   PMATRIX(rotateby);
 273
 274   printf("drag_rotate_delta...\n");
 275   show();
 276 }
 277 DRAG_SAVING(rotate, transform);
 278
 279 static void drag_scale_delta(double dx, double dy) {
 280   scale *= pow(3.0, dy);
 281   show();
 282 }
 283 DRAG_SAVING(scale, scale);
 284
 285 static const Drag *drag= &drag_none;
 286
 287 static int drag_last_x, drag_last_y;
 288
 289 static void drag_position(int x, int y) {
 290   drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
 291               (y - drag_last_y) * 1.0 / wmaxdim);
 292   drag_last_x= x;
 293   drag_last_y= y;
 294 }
 295
 296 static void event_button(XButtonEvent *e) {
 297   if (e->window != window || !e->same_screen) return;
 298   if (e->type == ButtonPress) {
 299     if (e->state || drag != &drag_none) {
 300       printf("drag=%s press state=0x%lx abandon\n",
 301              drag->name, (unsigned long)e->state);
 302       drag->abandon();
 303       drag= &drag_none;
 304       return;
 305     }
 306     switch (e->button) {
 307     case Button1: drag= &drag_rotate;  break;
 308     case Button2: drag= &drag_scale;   break;
 309     default: printf("unknown drag start %d\n", e->button);
 310     }
 311     printf("drag=%s press button=%lu start %d,%d\n",
 312            drag->name, (unsigned long)e->button, e->x, e->y);
 313     drag_last_x= e->x;
 314     drag_last_y= e->y;
 315     drag->start();
 316   }
 317   if (e->type == ButtonRelease) {
 318     printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
 319     drag_position(e->x, e->y);
 320     drag->conclude();
 321     drag= &drag_none;
 322   }
 323 }
 324
 325 static void event_motion(int x, int y) {
 326   printf("drag=%s motion %d,%d\n", drag->name, x, y);
 327   drag_position(x,y);
 328 }
 329
 330 static void event_config(XConfigureEvent *e) {
 331   if (e->width == wwidth && e->height == wheight)
 332     return;
 333
 334   wwidth= e->width;  wheight= e->height;
 335   wmaxdim= wwidth > wheight ? wwidth : wheight;
 336   wmindim= wwidth < wheight ? wwidth : wheight;
 337
 338   XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
 339   for (currentbuffer=0; currentbuffer<2; currentbuffer++)
 340     XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
 341
 342   mkpixmaps();
 343   show();
 344 }
 345
 346 int main(int argc, const char *const *argv) {
 347   XEvent event;
 348   int k;
 349   int motion_deferred=0, motion_x=-1, motion_y=-1;
 350
 351   if (argc != 2 || argv[1][0]=='-') {
 352     fputs("need filename\n",stderr); exit(8);
 353   }
 354   input_filename= argv[1];
 355
 356   read_input();
 357   K transform[k][k]= 1.0;
 358   display_prepare();
 359   show();
 360
 361   XMapWindow(display,window);
 362   for (;;) {
 363     if (motion_deferred) {
 364       int r= XCheckMaskEvent(display,~0UL,&event);
 365       if (!r) {
 366         event_motion(motion_x, motion_y);
 367         motion_deferred=0;
 368         continue;
 369       }
 370     } else {
 371       XNextEvent(display,&event);
 372     }
 373     switch (event.type) {
 374
 375     case ButtonPress:
 376     case ButtonRelease:     event_button(&event.xbutton);     break;
 377
 378     case ConfigureNotify:   event_config(&event.xconfigure);  break;
 379
 380     case MotionNotify:
 381 #if 0
 382       motion_x= event.xmotion.x;
 383       motion_y= event.xmotion.y;
 384       motion_deferred= 1;
 385 #endif
 386       continue;
 387
 388     default:
 389       printf("unknown event type %u 0x%x\n", event.type,event.type);
 390     }
 391   }
 392 }
 393