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