2 * Displays a conformation
16 typedef struct { double vertex[3][D3]; } Triangle;
18 static Triangle trisbuffer[MAXTRIS], *displaylist[MAXTRIS];
20 static Vertices conformation;
22 static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
23 static GSL_MATRIX(transform);
26 static struct stat input_stab;
27 static const char *input_filename;
29 static void read_input(void) {
32 if (input_f) fclose(input_f);
33 input_f= fopen(input_filename, "rb"); if (!input_f) diee("input file");
35 if (fstat(fileno(input_f), &input_stab)) diee("fstat input file");
38 r= fread(&conformation,sizeof(conformation),1,input_f);
39 if (r!=1) diee("fread");
42 static void transform_coordinates(void) {
45 gsl_vector input_gsl= { D3,1 };
50 input_gsl.data= &conformation[v][0];
51 GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
53 K conformation[v][k]= result[k];
57 static void addtriangle(int va, int vb, int vc) {
58 Triangle *t= &trisbuffer[ntris];
61 assert(ntris < MAXTRIS);
63 t->vertex[0][k]= conformation[va][k];
64 t->vertex[1][k]= conformation[vb][k];
65 t->vertex[2][k]= conformation[vc][k];
67 displaylist[ntris++]= t;
70 static void generate_display_list(void) {
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).
79 for (e=0; e<3; e++) ve[e]= EDGE_END2(vb,e);
81 if (ve[0]>=0) addtriangle(vb,ve[0],ve[1]);
82 if (ve[2]>=0) addtriangle(vb,ve[1],ve[2]);
87 static int dl_compare(const void *tav, const void *tbv) {
89 const Triangle *const *tap= tav, *ta= *tap;
90 const Triangle *const *tbp= tbp, *tb= *tbp;
93 za += ta->vertex[i][2];
94 zb += tb->vertex[i][2];
100 static void sort_display_list(void) {
101 qsort(displaylist, ntris, sizeof(*displaylist), dl_compare);
104 /*---------- X stuff ----------*/
108 typedef struct { GC fillgc, linegc; } DrawingMode;
110 static Display *display;
111 static Pixmap pixmap, doublebuffers[2];
112 static Window window;
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;
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;
126 static void drawtriangle(const Triangle *t) {
130 for (i=0; i<3; i++) {
131 double *v= t->vertex[i];
136 if (z < cut_z) { ncut++; return; }
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;
142 points[3]= points[0];
144 XA( XFillPolygon(display,pixmap, dmcurrent->fillgc,
145 points,3,Convex,CoordModeOrigin) );
146 XA( XDrawLines(display,pixmap, dmcurrent->linegc,
147 points, 4,CoordModeOrigin) );
150 static const unsigned long core_event_mask=
151 ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask|
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;
162 static void mkgcs(DrawingMode *dm, unsigned long planes) {
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,
172 gcv.function= GXclear;
173 dm->fillgc= XCreateGC(display,pixmap,
174 GCFunction|GCPlaneMask,
178 static void display_prepare(void) {
179 XSetWindowAttributes wa;
182 XA( display= XOpenDisplay(0) );
183 x11screen= DefaultScreen(display);
184 x11depth= DefaultDepth(display,x11screen);
185 XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&visinfo) );
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,
193 hints.flags= USPosition;
196 XSetWMNormalHints(display,window,&hints);
200 mkgcs(&dmwhite, AllPlanes);
201 mkgcs(&dmblue, visinfo.blue_mask);
202 mkgcs(&dmred, visinfo.red_mask);
205 static void drawtriangles(const DrawingMode *dm) {
210 for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
214 static void display_conformation(void) {
215 pixmap= doublebuffers[currentbuffer];
217 XA( XFillRectangle(display,pixmap,dmwhite.fillgc,0,0,wwidth,wheight) );
219 if (eyes_apart > 0) {
220 const double stationary= 0.07;
222 eye_x= eyes_apart < eyes_apart_preferred
224 eyes_apart < (eyes_apart_preferred + stationary)
225 ? eyes_apart_preferred
226 : eyes_apart - stationary;
227 eye_x /= sizeadj_scale;
228 drawtriangles(&dmblue);
230 drawtriangles(&dmred);
232 drawtriangles(&dmwhite);
233 printf("shown, %d/%d triangles cut\n", ncut, ntris);
236 XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
237 XA( XClearWindow(display,window) );
238 currentbuffer= !currentbuffer;
241 static void show(void) {
242 scale_wmindim= sizeadj_scale * wmindim;
244 transform_coordinates();
245 generate_display_list();
247 display_conformation();
252 void (*start)(const XButtonEvent *e);
253 void (*delta)(double dx, double dy);
254 void (*conclude)(void);
255 void (*abandon)(void);
259 static const Drag drag_##x= { \
260 #x, drag_##x##_start, drag_##x##_delta, \
261 drag_##x##_conclude, drag_##x##_abandon \
264 #define DRAG_SAVING(x, thing, hook) \
265 static typeof(thing) original_##thing; \
266 static void drag_##x##_start(const XButtonEvent *e) { \
267 memcpy(&original_##thing, &thing, sizeof(thing)); \
270 static void drag_##x##_conclude(void) { } \
271 static void drag_##x##_abandon(void) { \
272 memcpy(&thing, &original_##thing, sizeof(thing)); \
277 static void drag_none_start(const XButtonEvent *e) { }
278 static void drag_none_delta(double dx, double dy) { }
279 static void drag_none_conclude(void) { }
280 static void drag_none_abandon(void) { }
283 static void pvectorcore(const char *n, double v[D3]) {
286 K printf("%# 10.10f ",v[k]);
289 static void pvector(const char *n, double v[D3]) {
293 static void pmatrix(const char *n, double m[D3][D3]) {
295 for (j=0; j<D3; j++) { pvectorcore(n,m[j]); n=""; }
298 #define PMATRIX(x) pmatrix(#x,x);
300 static double drag_transform_conv_x_z= 0;
301 static double drag_transform_conv_y_z= 0;
303 static void drag_transform_prep(const XButtonEvent *e) {
304 static const double factor= 2.5;
305 drag_transform_conv_x_z= MAX( MIN(e->y * factor / wheight - (factor/2),
307 drag_transform_conv_y_z= MAX( MIN(e->x * factor / wwidth - (factor/2),
309 printf("drag_transform_conv_{x,y}_z = %g,%g\n",
310 drag_transform_conv_x_z, drag_transform_conv_y_z);
313 static void make_z_rotation(double rotz[D3][D3], double cz, double sz) {
314 rotz[0][0]= cz; rotz[0][1]= sz; rotz[0][2]= 0;
315 rotz[1][0]= -sz; rotz[1][1]= cz; rotz[1][2]= 0;
316 rotz[2][0]= 0; rotz[2][1]= 0; rotz[2][2]= 1;
319 static void drag_rotate_delta(double dx, double dy) {
320 /* We multiple our transformation matrix by a matrix:
322 * If we just had y movement, we would rotate about x axis:
323 * rotation X = [ 1 0 0 ]
326 * where cy and sy are sin and cos of y rotation
328 * But we should pre-rotate this by a rotation about the z axis
329 * to get it to the right angle (to include x rotation). So
330 * we make cy and sy be cos() and sin(hypot(x,y)) and use
331 * with cr,sr as cos() and sin(atan2(y,y)):
333 * Ie we would do T' = Z R^T X R T where
334 * or T' = Z C T where C = R^T X R and
336 * adjustment R = [ cr sr 0 ]
338 * [ 0 0 1 ] or make_z_rotation(cr,sr)
340 * rotation Z = [ cz sz 0 ]
342 * [ 0 0 1 ] or make_z_rotation(cz,sz)
345 double rotx[D3][D3], adjr[D3][D3], rotz[D3][D3];
350 static double temp[D3][D3], change[D3][D3];
351 static GSL_MATRIX(temp);
352 static GSL_MATRIX(change);
354 printf("\nTRANSFORM %g, %g\n", dx,dy);
356 double dz= -drag_transform_conv_x_z * dx +
357 drag_transform_conv_y_z * dy;
359 dx *= (1 - fabs(drag_transform_conv_x_z));
360 dy *= (1 - fabs(drag_transform_conv_y_z));
362 double d= hypot(dx,dy);
364 printf(" dx,dy,dz = %g, %g, %g d = %g\n", dx,dy,dz, d);
366 if (hypot(d,dz) < 1e-6) return;
368 printf(" no xy rotation\n");
383 printf("\n d=%g cy,sy=%g,%g cr,sr=%g,%g\n\n", d,cy,sy,cr,sr);
385 rotx[0][0]= 1; rotx[0][1]= 0; rotx[0][2]= 0;
386 rotx[1][0]= 0; rotx[1][1]= cy; rotx[1][2]= sy;
387 rotx[2][0]= 0; rotx[2][1]= -sy; rotx[2][2]= cy;
390 make_z_rotation(adjr,cr,sr);
393 GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
396 pmatrix("X R", temp);
398 GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
404 make_z_rotation(rotz,cos(angz),sin(angz));
407 GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
408 &rotz_gsl,&change_gsl,
410 pmatrix("Z C", temp);
412 static double skew[D3][D3];
413 static GSL_MATRIX(skew);
415 GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
416 &temp_gsl, &transform_gsl,
420 memcpy(&transform,&skew,sizeof(transform));
424 /* Now we want to normalise skew, the result becomes new transform */
425 double svd_v[D3][D3];
428 double sigma[D3], tau[D3];
432 /* We use notation from Wikipedia Polar_decomposition
433 * Wikipedia's W is GSL's U
434 * Wikipedia's Sigma is GSL's S
435 * Wikipedia's V is GSL's V
436 * Wikipedia's U is our desired result
437 * Wikipedia which says if the SVD is A = W Sigma V*
438 * then the polar decomposition is A = U P
439 * where P = V Sigma V*
443 GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
445 pvector("Sigma",sigma);
448 /* We only need U, not P. */
449 GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
450 &skew_gsl,&svd_v_gsl,
451 0.0,&transform_gsl) );
453 pmatrix("U", transform);
455 printf("drag_rotate_delta...\n");
458 DRAG_SAVING(rotate, transform, drag_transform_prep(e));
460 static void drag_sizeadj_delta(double dx, double dy) {
461 sizeadj_scale *= pow(3.0, -dy);
464 DRAG_SAVING(sizeadj, sizeadj_scale, );
466 static void drag_3d_delta(double dx, double dy) {
467 eyes_apart += dx * 0.1;
468 if (eyes_apart < eyes_apart_min) eyes_apart= eyes_apart_min;
469 printf("sizeadj eyes_apart %g\n", eyes_apart);
472 DRAG_SAVING(3d, eyes_apart, );
474 static const Drag *drag= &drag_none;
476 static int drag_last_x, drag_last_y;
478 static void drag_position(int x, int y) {
479 drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
480 (y - drag_last_y) * 1.0 / wmaxdim);
485 static void event_button(XButtonEvent *e) {
486 if (e->window != window || !e->same_screen) return;
487 if (e->type == ButtonPress) {
488 if (e->state || drag != &drag_none) {
489 printf("drag=%s press state=0x%lx abandon\n",
490 drag->name, (unsigned long)e->state);
496 case Button1: drag= &drag_rotate; break;
497 case Button2: drag= &drag_sizeadj; break;
498 case Button3: drag= &drag_3d; break;
499 default: printf("unknown drag start %d\n", e->button);
501 printf("drag=%s press button=%lu start %d,%d\n",
502 drag->name, (unsigned long)e->button, e->x, e->y);
507 if (e->type == ButtonRelease) {
508 printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
509 drag_position(e->x, e->y);
515 static void event_motion(int x, int y) {
516 printf("drag=%s motion %d,%d\n", drag->name, x, y);
520 static void event_key(XKeyEvent *e) {
525 r= XLookupString(e,buf,sizeof(buf)-1,&ks,0);
527 printf("XLookupString keycode=%u state=0x%x gave %d\n",
528 e->keycode, e->state, r);
532 if (!strcmp(buf,"q")) exit(0);
533 if (!strcmp(buf,"d")) {
534 eyes_apart= eyes_apart>0 ? eyes_apart_min : eyes_apart_preferred;
540 static void event_config(XConfigureEvent *e) {
541 if (e->width == wwidth && e->height == wheight)
544 wwidth= e->width; wheight= e->height;
545 wmaxdim= wwidth > wheight ? wwidth : wheight;
546 wmindim= wwidth < wheight ? wwidth : wheight;
548 XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
549 for (currentbuffer=0; currentbuffer<2; currentbuffer++)
550 XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
556 static void check_input(void) {
560 r= stat(input_filename, &newstab);
561 if (r<0) diee("could not check input");
563 #define CI(x) if (newstab.st_##x == input_stab.st_##x) ; else goto changed
575 int main(int argc, const char *const *argv) {
576 static const int wantedevents= POLLIN|POLLPRI|POLLERR|POLLHUP;
579 int k, i, r, *xfds, nxfds, polls_alloc=0;
580 struct pollfd *polls=0;
581 int motion_deferred=0, motion_x=-1, motion_y=-1;
583 if (argc != 2 || argv[1][0]=='-') {
584 fputs("need filename\n",stderr); exit(8);
586 input_filename= argv[1];
589 K transform[k][k]= 1.0;
593 XMapWindow(display,window);
596 XA( XInternalConnectionNumbers(display, &xfds, &nxfds) );
597 if (polls_alloc <= nxfds) {
598 polls_alloc= nxfds + polls_alloc + 1;
599 polls= realloc(polls, sizeof(*polls) * polls_alloc);
600 if (!polls) diee("realloc for pollfds");
602 for (i=0; i<nxfds; i++) {
603 polls[i].fd= xfds[i];
604 polls[i].events= wantedevents;
609 polls[i].fd= ConnectionNumber(display);
610 polls[i].events= wantedevents;
612 r= poll(polls, nxfds+1, motion_deferred ? 0 : 200);
613 if (r<0) diee("poll");
615 for (i=0; i<nxfds; i++)
616 if (polls[i].revents)
617 XProcessInternalConnection(display, polls[i].fd);
619 r= XCheckMaskEvent(display,~0UL,&event);
621 if (motion_deferred) {
622 event_motion(motion_x, motion_y);
629 switch (event.type) {
632 case ButtonRelease: event_button(&event.xbutton); break;
634 case KeyPress: event_key(&event.xkey); break;
636 case ConfigureNotify: event_config(&event.xconfigure); break;
639 motion_x= event.xmotion.x;
640 motion_y= event.xmotion.y;
645 printf("unknown event type %u 0x%x\n", event.type,event.type);