2 * Displays a conformation
14 typedef struct { double vertex[3][D3]; } Triangle;
16 static Triangle trisbuffer[MAXTRIS], *displaylist[MAXTRIS];
18 static Vertices conformation;
20 static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
21 static GSL_MATRIX(transform);
24 static struct stat input_stab;
25 static const char *input_filename;
26 static int pause_updates;
29 static char title[TITLE_MAX+1];
31 static void read_input(void) {
32 int r, last_readlink_r, last_readlink_errno=0;
33 char symlink_check[sizeof(title)];
38 r= readlink(input_filename, symlink_check, TITLE_MAX);
42 if (r>=0) symlink_check[r]= 0;
43 else symlink_check[0]= 0;
45 if (r == last_readlink_r &&
47 ? (errno==last_readlink_errno)
48 : !strcmp(symlink_check, title))) {
51 if (last_readlink_errno==EINVAL)
52 snprintf(title,sizeof(title),"%s",input_filename);
54 snprintf(title,sizeof(title),"? %s",strerror(last_readlink_errno));
59 strcpy(title, symlink_check);
61 last_readlink_errno= errno;
63 if (input_f) fclose(input_f);
64 input_f= fopen(input_filename, "rb"); if (!input_f) diee("input file");
67 if (fstat(fileno(input_f), &input_stab)) diee("fstat input file");
70 r= fread(&conformation,sizeof(conformation),1,input_f);
71 if (r!=1) diee("fread");
74 static void transform_coordinates(void) {
77 gsl_vector input_gsl= { D3,1 };
82 input_gsl.data= &conformation[v][0];
83 GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
85 K conformation[v][k]= result[k];
89 static int vertex_in_triangles[N], vertex_in_triangles_checked;
91 static void addtriangle(int va, int vb, int vc) {
92 Triangle *t= &trisbuffer[ntris];
95 assert(ntris < MAXTRIS);
97 t->vertex[0][k]= conformation[va][k];
98 t->vertex[1][k]= conformation[vb][k];
99 t->vertex[2][k]= conformation[vc][k];
101 if (!vertex_in_triangles_checked) {
102 vertex_in_triangles[va]++;
103 vertex_in_triangles[vb]++;
104 vertex_in_triangles[vc]++;
106 displaylist[ntris++]= t;
109 static void generate_display_list(void) {
113 FOR_VERTEX(vb,INNER) {
114 /* We use the two triangles in the parallelogram vb, vb+e5, vb+e0, vb+e1.
115 * We go round each triangle clockwise (although our surface is non-
116 * orientable so it shouldn't matter). Picking the parallelogram
117 * to our right avoids getting it wrong at the join.
119 //if ((vb & YMASK) > Y1) continue;
120 //if ((vb & XMASK) > 2) continue;
121 for (e=0; e<V6; e++) ve[e]= EDGE_END2(vb,e);
123 if (ve[5]>=0) addtriangle(vb,ve[0],ve[5]);
125 if (ve[1]>=0) addtriangle(vb,ve[1],ve[0]);
128 if (!vertex_in_triangles_checked) {
130 FOR_VERTEX(v,INNER) {
131 expd= RIM_VERTEX_P(v) ? 3 : 6;
132 if (vertex_in_triangles[v] != expd) {
133 fprintf(stderr,"vertex %02x used for %d triangles, expected %d\n",
134 v, vertex_in_triangles[v], expd);
138 vertex_in_triangles_checked= 1;
142 static int dl_compare(const void *tav, const void *tbv) {
144 const Triangle *const *tap= tav, *ta= *tap;
145 const Triangle *const *tbp= tbv, *tb= *tbp;
147 for (i=0; i<3; i++) {
148 za += ta->vertex[i][2];
149 zb += tb->vertex[i][2];
151 return za > zb ? -1 :
155 static void sort_display_list(void) {
156 qsort(displaylist, ntris, sizeof(*displaylist), dl_compare);
159 /*---------- X stuff ----------*/
163 typedef struct { GC fillgc, linegc; } DrawingMode;
165 static Display *display;
166 static Pixmap pixmap, doublebuffers[2];
167 static Window window;
169 static DrawingMode dmred, dmblue, dmwhite;
170 static const DrawingMode *dmcurrent;
171 static int wwidth=WSZ, wheight=WSZ, wmindim=WSZ, wmaxdim=WSZ;
172 static int ncut, currentbuffer, x11depth, x11screen, wireframe;
175 static double sizeadj_scale= 0.3, eyes_apart, scale_wmindim;
176 static double eye_z= -10, eye_x=0;
177 static double cut_z= -9;
178 static const double eyes_apart_preferred=0.05, eyes_apart_min= -0.02;
181 static void drawtriangle(const Triangle *t) {
185 for (i=0; i<3; i++) {
186 const double *v= t->vertex[i];
191 if (z < cut_z) { ncut++; return; }
193 double zezezp= eye_z / (eye_z - z);
194 points[i].x= scale_wmindim * (zezezp * (x - eye_x) + eye_x) + wwidth/2;
195 points[i].y= scale_wmindim * (zezezp * y ) + wheight/2;
197 points[3]= points[0];
200 XA( XFillPolygon(display,pixmap, dmcurrent->fillgc,
201 points,3,Convex,CoordModeOrigin) );
202 XA( XDrawLines(display,pixmap, dmcurrent->linegc,
203 points, 4,CoordModeOrigin) );
206 static const unsigned long core_event_mask=
207 ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask|
208 KeyPressMask|SubstructureNotifyMask;
210 static void mkpixmaps(void) {
211 for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
212 XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,x11depth) );
213 doublebuffers[currentbuffer]= pixmap;
218 static void mkgcs(DrawingMode *dm, unsigned long planes) {
221 gcv.function= GXcopy;
222 gcv.foreground= WhitePixel(display,x11screen);
223 gcv.plane_mask= planes;
224 dm->linegc= XCreateGC(display,pixmap,
225 GCFunction|GCForeground|GCPlaneMask,
228 gcv.function= GXclear;
229 dm->fillgc= XCreateGC(display,pixmap,
230 GCFunction|GCPlaneMask,
234 static void display_prepare(void) {
235 XSetWindowAttributes wa;
238 XA( display= XOpenDisplay(0) );
239 x11screen= DefaultScreen(display);
240 x11depth= DefaultDepth(display,x11screen);
241 XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&visinfo) );
243 wa.event_mask= core_event_mask;
244 XA( window= XCreateWindow(display, DefaultRootWindow(display),
245 0,0, wwidth,wheight, 0,x11depth,
246 InputOutput, visinfo.visual,
249 hints.flags= USPosition;
252 XSetWMNormalHints(display,window,&hints);
256 mkgcs(&dmwhite, AllPlanes);
257 mkgcs(&dmblue, visinfo.blue_mask);
258 mkgcs(&dmred, visinfo.red_mask);
261 static void drawtriangles(const DrawingMode *dm) {
266 for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
270 static void display_conformation(void) {
271 pixmap= doublebuffers[currentbuffer];
273 XA( XFillRectangle(display,pixmap,dmwhite.fillgc,0,0,wwidth,wheight) );
275 if (eyes_apart > 0) {
276 const double stationary= 0.07;
278 eye_x= eyes_apart < eyes_apart_preferred
280 eyes_apart < (eyes_apart_preferred + stationary)
281 ? eyes_apart_preferred
282 : eyes_apart - stationary;
283 eye_x /= sizeadj_scale;
284 drawtriangles(&dmblue);
286 drawtriangles(&dmred);
288 drawtriangles(&dmwhite);
289 printf("shown, %d/%d triangles cut\n", ncut, ntris);
292 XA( XStoreName(display,window,title) );
293 XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
294 XA( XClearWindow(display,window) );
295 currentbuffer= !currentbuffer;
298 static void show(void) {
299 scale_wmindim= sizeadj_scale * wmindim;
301 transform_coordinates();
302 generate_display_list();
304 display_conformation();
309 void (*start)(const XButtonEvent *e);
310 void (*delta)(double dx, double dy);
311 void (*conclude)(void);
312 void (*abandon)(void);
316 static const Drag drag_##x= { \
317 #x, drag_##x##_start, drag_##x##_delta, \
318 drag_##x##_conclude, drag_##x##_abandon \
321 #define DRAG_SAVING(x, thing, hook) \
322 static typeof(thing) original_##thing; \
323 static void drag_##x##_start(const XButtonEvent *e) { \
324 memcpy(&original_##thing, &thing, sizeof(thing)); \
327 static void drag_##x##_conclude(void) { } \
328 static void drag_##x##_abandon(void) { \
329 memcpy(&thing, &original_##thing, sizeof(thing)); \
334 static void drag_none_start(const XButtonEvent *e) { }
335 static void drag_none_delta(double dx, double dy) { }
336 static void drag_none_conclude(void) { }
337 static void drag_none_abandon(void) { }
340 static void pvectorcore(const char *n, double v[D3]) {
343 K printf("%# 10.10f ",v[k]);
346 static void pvector(const char *n, double v[D3]) {
350 static void pmatrix(const char *n, double m[D3][D3]) {
352 for (j=0; j<D3; j++) { pvectorcore(n,m[j]); n=""; }
355 #define PMATRIX(x) pmatrix(#x,x);
357 static double drag_transform_conv_x_z= 0;
358 static double drag_transform_conv_y_z= 0;
360 static void drag_transform_prep(const XButtonEvent *e) {
361 static const double factor= 2.5;
362 drag_transform_conv_x_z= MAX( MIN(e->y * factor / wheight - (factor/2),
364 drag_transform_conv_y_z= MAX( MIN(e->x * factor / wwidth - (factor/2),
366 printf("drag_transform_conv_{x,y}_z = %g,%g\n",
367 drag_transform_conv_x_z, drag_transform_conv_y_z);
370 static void make_z_rotation(double rotz[D3][D3], double cz, double sz) {
371 rotz[0][0]= cz; rotz[0][1]= sz; rotz[0][2]= 0;
372 rotz[1][0]= -sz; rotz[1][1]= cz; rotz[1][2]= 0;
373 rotz[2][0]= 0; rotz[2][1]= 0; rotz[2][2]= 1;
376 static void drag_rotate_delta(double dx, double dy) {
377 /* We multiple our transformation matrix by a matrix:
379 * If we just had y movement, we would rotate about x axis:
380 * rotation X = [ 1 0 0 ]
383 * where cy and sy are sin and cos of y rotation
385 * But we should pre-rotate this by a rotation about the z axis
386 * to get it to the right angle (to include x rotation). So
387 * we make cy and sy be cos() and sin(hypot(x,y)) and use
388 * with cr,sr as cos() and sin(atan2(y,y)):
390 * Ie we would do T' = Z R^T X R T where
391 * or T' = Z C T where C = R^T X R and
393 * adjustment R = [ cr sr 0 ]
395 * [ 0 0 1 ] or make_z_rotation(cr,sr)
397 * rotation Z = [ cz sz 0 ]
399 * [ 0 0 1 ] or make_z_rotation(cz,sz)
402 double rotx[D3][D3], adjr[D3][D3], rotz[D3][D3];
407 static double temp[D3][D3], change[D3][D3];
408 static GSL_MATRIX(temp);
409 static GSL_MATRIX(change);
411 printf("\nTRANSFORM %g, %g\n", dx,dy);
413 double dz= -drag_transform_conv_x_z * dx +
414 drag_transform_conv_y_z * dy;
416 dx *= (1 - fabs(drag_transform_conv_x_z));
417 dy *= (1 - fabs(drag_transform_conv_y_z));
419 double d= hypot(dx,dy);
421 printf(" dx,dy,dz = %g, %g, %g d = %g\n", dx,dy,dz, d);
423 if (hypot(d,dz) < 1e-6) return;
425 printf(" no xy rotation\n");
440 printf("\n d=%g cy,sy=%g,%g cr,sr=%g,%g\n\n", d,cy,sy,cr,sr);
442 rotx[0][0]= 1; rotx[0][1]= 0; rotx[0][2]= 0;
443 rotx[1][0]= 0; rotx[1][1]= cy; rotx[1][2]= sy;
444 rotx[2][0]= 0; rotx[2][1]= -sy; rotx[2][2]= cy;
447 make_z_rotation(adjr,cr,sr);
450 GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
453 pmatrix("X R", temp);
455 GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
461 make_z_rotation(rotz,cos(angz),sin(angz));
464 GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
465 &rotz_gsl,&change_gsl,
467 pmatrix("Z C", temp);
469 static double skew[D3][D3];
470 static GSL_MATRIX(skew);
472 GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
473 &temp_gsl, &transform_gsl,
477 memcpy(&transform,&skew,sizeof(transform));
481 /* Now we want to normalise skew, the result becomes new transform */
482 double svd_v[D3][D3];
485 double sigma[D3], tau[D3];
489 /* We use notation from Wikipedia Polar_decomposition
490 * Wikipedia's W is GSL's U
491 * Wikipedia's Sigma is GSL's S
492 * Wikipedia's V is GSL's V
493 * Wikipedia's U is our desired result
494 * Wikipedia which says if the SVD is A = W Sigma V*
495 * then the polar decomposition is A = U P
496 * where P = V Sigma V*
500 GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
502 pvector("Sigma",sigma);
505 /* We only need U, not P. */
506 GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
507 &skew_gsl,&svd_v_gsl,
508 0.0,&transform_gsl) );
510 pmatrix("U", transform);
512 printf("drag_rotate_delta...\n");
515 DRAG_SAVING(rotate, transform, drag_transform_prep(e));
517 static void drag_sizeadj_delta(double dx, double dy) {
518 sizeadj_scale *= pow(3.0, -dy);
521 DRAG_SAVING(sizeadj, sizeadj_scale, );
523 static void drag_3d_delta(double dx, double dy) {
524 eyes_apart += dx * 0.1;
525 if (eyes_apart < eyes_apart_min) eyes_apart= eyes_apart_min;
526 printf("sizeadj eyes_apart %g\n", eyes_apart);
529 DRAG_SAVING(3d, eyes_apart, );
531 static const Drag *drag= &drag_none;
533 static int drag_last_x, drag_last_y;
535 static void drag_position(int x, int y) {
536 drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
537 (y - drag_last_y) * 1.0 / wmaxdim);
542 static void event_button(XButtonEvent *e) {
543 if (e->window != window || !e->same_screen) return;
544 if (e->type == ButtonPress) {
545 if (e->state || drag != &drag_none) {
546 printf("drag=%s press state=0x%lx abandon\n",
547 drag->name, (unsigned long)e->state);
553 case Button1: drag= &drag_rotate; break;
554 case Button2: drag= &drag_sizeadj; break;
555 case Button3: drag= &drag_3d; break;
556 default: printf("unknown drag start %d\n", e->button);
558 printf("drag=%s press button=%lu start %d,%d\n",
559 drag->name, (unsigned long)e->button, e->x, e->y);
564 if (e->type == ButtonRelease) {
565 printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
566 drag_position(e->x, e->y);
572 static void event_motion(int x, int y) {
573 printf("drag=%s motion %d,%d\n", drag->name, x, y);
577 static void transform_preset_record(const char *fn, const char *fn_new) {
579 f= fopen(fn_new,"wb");
580 if (!f) diee("open new transform");
581 if (fwrite(transform,sizeof(transform),1,f) != 1) diee("write transform");
582 if (fclose(f)) diee("fclose new transform");
583 if (rename(fn_new,fn)) diee("install transform");
586 static void transform_preset_playback(const char *fn) {
589 if (!f && errno==ENOENT) {
590 fprintf(stderr,"no preset %s\n",fn);
595 if (fread(transform,sizeof(transform),1,f) != 1) {
596 perror("read preset!");
604 static void event_key(XKeyEvent *e) {
607 char buf[10], buf_nomod[10];
610 r= XLookupString(e,buf,sizeof(buf)-1,&ks,0);
612 printf("XLookupString keycode=%u state=0x%x gave %d\n",
613 e->keycode, e->state, r);
617 if (!strcmp(buf,"q"))
619 else if (!strcmp(buf,"p"))
620 pause_updates= !pause_updates;
621 else if (!strcmp(buf,"w")) {
622 wireframe= !wireframe;
625 } else if (!strcmp(buf,"d")) {
626 eyes_apart= eyes_apart>0 ? eyes_apart_min : eyes_apart_preferred;
634 r_nomod= XLookupString(&e_nomod,buf_nomod,sizeof(buf_nomod)-1,&ks,0);
635 if (r_nomod && !buf_nomod[1] && buf_nomod[0]>='0' && buf_nomod[0]<='9') {
636 char filename[20], filename_new[25];
637 snprintf(filename,sizeof(filename)-1,".view-preset-%s",buf_nomod);
638 snprintf(filename_new,sizeof(filename_new)-1,"%s.new",filename);
639 printf("transform preset %d %s\n", e->state, filename);
640 if (e->state) transform_preset_record(filename,filename_new);
641 else transform_preset_playback(filename);
645 printf("unknown key keycode=%d state=0x%x char=%c 0x%02x "
646 "[rnm=%d bnm[0,1]=0x%02x,%02x]\n",
647 e->keycode, e->state, buf[0]>' ' && buf[0]<127 ? buf[0] : '?',
648 buf[0], r_nomod, buf_nomod[0], buf_nomod[1]);
649 printf("%d %d %d %d\n",
656 static void event_config(XConfigureEvent *e) {
657 if (e->width == wwidth && e->height == wheight)
660 wwidth= e->width; wheight= e->height;
661 wmaxdim= wwidth > wheight ? wwidth : wheight;
662 wmindim= wwidth < wheight ? wwidth : wheight;
664 XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
665 for (currentbuffer=0; currentbuffer<2; currentbuffer++)
666 XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
672 static void check_input(void) {
676 r= stat(input_filename, &newstab);
677 if (r<0) diee("could not check input");
679 #define CI(x) if (newstab.st_##x == input_stab.st_##x) ; else goto changed
691 static void topocheck(void) {
692 int v1,e,v2,eprime,v1prime, count;
693 FOR_EDGE(v1,e,v2, INNER) {
695 FOR_VEDGE(v2,eprime,v1prime)
696 if (v1prime==v1) count++;
698 fprintf(stderr,"%02x -%d-> %02x reverse edge count = %d!\n",
700 FOR_VEDGE(v2,eprime,v1prime)
701 fprintf(stderr,"%02x -%d-> %02x -> %d -> %02x\n",
702 v1,e,v2,eprime,v1prime);
708 int main(int argc, const char *const *argv) {
709 static const int wantedevents= POLLIN|POLLPRI|POLLERR|POLLHUP;
712 int k, i, r, *xfds, nxfds, polls_alloc=0;
713 struct pollfd *polls=0;
714 int motion_deferred=0, motion_x=-1, motion_y=-1;
718 if (argc==1) { printf("topology self-consistent, ok\n"); exit(0); }
720 if (argc != 2 || argv[1][0]=='-') {
721 fputs("need filename\n",stderr); exit(8);
723 input_filename= argv[1];
726 K transform[k][k]= 1.0;
730 XMapWindow(display,window);
733 XA( XInternalConnectionNumbers(display, &xfds, &nxfds) );
734 if (polls_alloc <= nxfds) {
735 polls_alloc= nxfds + polls_alloc + 1;
736 polls= realloc(polls, sizeof(*polls) * polls_alloc);
737 if (!polls) diee("realloc for pollfds");
739 for (i=0; i<nxfds; i++) {
740 polls[i].fd= xfds[i];
741 polls[i].events= wantedevents;
746 polls[i].fd= ConnectionNumber(display);
747 polls[i].events= wantedevents;
749 r= poll(polls, nxfds+1, motion_deferred ? 0 : pause_updates ? -1 : 200);
751 if (errno==EINTR) continue;
755 for (i=0; i<nxfds; i++)
756 if (polls[i].revents)
757 XProcessInternalConnection(display, polls[i].fd);
759 r= XCheckMaskEvent(display,~0UL,&event);
761 if (motion_deferred) {
762 event_motion(motion_x, motion_y);
770 switch (event.type) {
773 case ButtonRelease: event_button(&event.xbutton); break;
775 case KeyPress: event_key(&event.xkey); break;
777 case ConfigureNotify: event_config(&event.xconfigure); break;
780 motion_x= event.xmotion.x;
781 motion_y= event.xmotion.y;
786 printf("unknown event type %u 0x%x\n", event.type,event.type);