2 * Displays a conformation
12 typedef struct { double vertex[3][D3]; } Triangle;
14 static Triangle trisbuffer[MAXTRIS], *displaylist[MAXTRIS];
16 static Vertices conformation;
18 static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
19 static GSL_MATRIX(transform);
21 const char *input_filename;
23 static void read_input(void) {
27 f= fopen(input_filename, "rb"); if (!f) diee("input file");
29 r= fread(&conformation,sizeof(conformation),1,f); if (r!=1) diee("fread");
33 static void transform_coordinates(void) {
36 gsl_vector input_gsl= { D3,1 };
41 input_gsl.data= &conformation[v][0];
42 GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
44 K conformation[v][k]= result[k];
48 static void addtriangle(int va, int vb, int vc) {
49 Triangle *t= &trisbuffer[ntris];
52 assert(ntris < MAXTRIS);
54 t->vertex[0][k]= conformation[va][k];
55 t->vertex[1][k]= conformation[vb][k];
56 t->vertex[2][k]= conformation[vc][k];
58 displaylist[ntris++]= t;
61 static void generate_display_list(void) {
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).
70 for (e=0; e<3; e++) ve[e]= EDGE_END2(vb,e);
72 if (ve[1]>=0) addtriangle(vb,ve[0],ve[1]);
73 if (ve[2]>=0) addtriangle(vb,ve[2],ve[0]);
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];
87 static void sort_display_list(void) {
88 qsort(displaylist, ntris, sizeof(*displaylist), dl_compare);
91 /*---------- X stuff ----------*/
95 typedef struct { GC fillgc, linegc; } DrawingMode;
97 static Display *display;
98 static Pixmap pixmap, doublebuffers[2];
101 static DrawingMode dmred, dmblue, dmwhite;
102 static const DrawingMode *dmcurrent;
103 static int wwidth=WSZ, wheight=WSZ, wmindim=WSZ, wmaxdim=WSZ;
104 static int ncut, currentbuffer, x11depth, x11screen;
107 static struct { double scale, eyes_apart; } sizeadj= { 0.3, 0.0 };
108 static double scale_wmindim;
109 static double eye_z= -10, eye_x=0;
110 static double cut_z= -9;
112 static void drawtriangle(const Triangle *t) {
116 for (i=0; i<3; i++) {
117 double *v= t->vertex[i];
122 if (z < cut_z) { ncut++; return; }
124 double zezezp= eye_z / (eye_z - z);
125 points[i].x= scale_wmindim * (zezezp * (x - eye_x) + eye_x) + wwidth/2;
126 points[i].y= scale_wmindim * (zezezp * y ) + wheight/2;
128 points[3]= points[0];
130 XA( XFillPolygon(display,pixmap, dmcurrent->fillgc,
131 points,3,Convex,CoordModeOrigin) );
132 XA( XDrawLines(display,pixmap, dmcurrent->linegc,
133 points, 4,CoordModeOrigin) );
136 static const unsigned long core_event_mask=
137 ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask;
139 static void mkpixmaps(void) {
140 for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
141 XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,x11depth) );
142 doublebuffers[currentbuffer]= pixmap;
147 static void mkgcs(DrawingMode *dm, unsigned long planes) {
150 gcv.function= GXcopy;
151 gcv.foreground= WhitePixel(display,x11screen);
152 gcv.plane_mask= planes;
153 dm->linegc= XCreateGC(display,pixmap,
154 GCFunction|GCForeground|GCPlaneMask,
157 gcv.function= GXclear;
158 dm->fillgc= XCreateGC(display,pixmap,
159 GCFunction|GCPlaneMask,
163 static void display_prepare(void) {
164 XSetWindowAttributes wa;
167 XA( display= XOpenDisplay(0) );
168 x11screen= DefaultScreen(display);
169 x11depth= DefaultDepth(display,x11screen);
170 XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&visinfo) );
172 wa.event_mask= core_event_mask;
173 XA( window= XCreateWindow(display, DefaultRootWindow(display),
174 0,0, wwidth,wheight, 0,x11depth,
175 InputOutput, visinfo.visual,
178 hints.flags= USPosition;
181 XSetWMNormalHints(display,window,&hints);
185 mkgcs(&dmwhite, AllPlanes);
186 mkgcs(&dmblue, visinfo.blue_mask);
187 mkgcs(&dmred, visinfo.red_mask);
190 static void drawtriangles(const DrawingMode *dm) {
195 for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
199 static void display_conformation(void) {
200 pixmap= doublebuffers[currentbuffer];
202 XA( XFillRectangle(display,pixmap,dmwhite.fillgc,0,0,wwidth,wheight) );
204 if (sizeadj.eyes_apart > 0) {
205 const double preferred=100, beyond=200;
207 eye_x= sizeadj.eyes_apart < preferred ? sizeadj.eyes_apart :
208 sizeadj.eyes_apart < beyond ? preferred :
209 sizeadj.eyes_apart - (beyond - preferred);
210 eye_x /= scale_wmindim;
211 drawtriangles(&dmblue);
213 drawtriangles(&dmred);
215 drawtriangles(&dmwhite);
216 printf("shown, %d/%d triangles cut\n", ncut, ntris);
219 XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
220 XA( XClearWindow(display,window) );
221 currentbuffer= !currentbuffer;
224 static void show(void) {
225 scale_wmindim= sizeadj.scale * wmindim;
227 transform_coordinates();
228 generate_display_list();
230 display_conformation();
236 void (*delta)(double dx, double dy);
237 void (*conclude)(void);
238 void (*abandon)(void);
242 static const Drag drag_##x= { \
243 #x, drag_##x##_start, drag_##x##_delta, \
244 drag_##x##_conclude, drag_##x##_abandon \
247 #define DRAG_SAVING(x, thing) \
248 static typeof(thing) original_##thing; \
249 static void drag_##x##_start(void) { \
250 memcpy(&original_##thing, &thing, sizeof(thing)); \
252 static void drag_##x##_conclude(void) { } \
253 static void drag_##x##_abandon(void) { \
254 memcpy(&thing, &original_##thing, sizeof(thing)); \
259 static void drag_none_start(void) { }
260 static void drag_none_delta(double dx, double dy) { }
261 static void drag_none_conclude(void) { }
262 static void drag_none_abandon(void) { }
265 static void pvectorcore(const char *n, double v[D3]) {
268 K printf("%# 10.10f ",v[k]);
271 static void pvector(const char *n, double v[D3]) {
275 static void pmatrix(const char *n, double m[D3][D3]) {
277 for (j=0; j<D3; j++) { pvectorcore(n,m[j]); n=""; }
280 #define PMATRIX(x) pmatrix(#x,x);
282 static void drag_rotate_delta(double dx, double dy) {
283 /* We multiple our transformation matrix by a matrix:
285 * If we just had y movement, we would rotate about x axis:
286 * rotation X = [ 1 0 0 ]
289 * where cy and sy are sin and cos of y rotation
291 * But we should pre-rotate this by a rotation about the z axis
292 * to get it to the right angle (to include x rotation). So
293 * we make cy and sy be cos() and sin(hypot(x,y)) and use
294 * with cr,sr as cos() and sin(atan2(y,y)):
296 * Ie we would do T' = R^T X R T where
297 * or T' = C T where C = R^T X R and
299 * adjustment R = [ cr sr 0 ]
304 double rotx[D3][D3], adjr[D3][D3];
308 static double temp[D3][D3], change[D3][D3];
309 static GSL_MATRIX(temp);
310 static GSL_MATRIX(change);
312 double d= hypot(dx,dy);
313 if (d < 1e-6) return;
321 printf("\n d=%g cy,sy=%g,%g cr,sr=%g,%g\n\n", d,cy,sy,cr,sr);
323 rotx[0][0]= 1; rotx[0][1]= 0; rotx[0][2]= 0;
324 rotx[1][0]= 0; rotx[1][1]= cy; rotx[1][2]= sy;
325 rotx[2][0]= 0; rotx[2][1]= -sy; rotx[2][2]= cy;
328 adjr[0][0]= cr; adjr[0][1]= sr; adjr[0][2]= 0;
329 adjr[1][0]= -sr; adjr[1][1]= cr; adjr[1][2]= 0;
330 adjr[2][0]= 0; adjr[2][1]= 0; adjr[2][2]= 1;
333 GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
338 GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
343 static double skew[D3][D3];
344 static GSL_MATRIX(skew);
346 GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
347 &change_gsl,&transform_gsl,
351 memcpy(&transform,&skew,sizeof(transform));
355 /* Now we want to normalise skew, the result becomes new transform */
356 double svd_v[D3][D3];
359 double sigma[D3], tau[D3];
363 /* We use notation from Wikipedia Polar_decomposition
364 * Wikipedia's W is GSL's U
365 * Wikipedia's Sigma is GSL's S
366 * Wikipedia's V is GSL's V
367 * Wikipedia's U is our desired result
368 * Wikipedia which says if the SVD is A = W Sigma V*
369 * then the polar decomposition is A = U P
370 * where P = V Sigma V*
374 GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
376 pvector("Sigma",sigma);
379 /* We only need U, not P. */
380 GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
381 &skew_gsl,&svd_v_gsl,
382 0.0,&transform_gsl) );
384 pmatrix("U", transform);
386 printf("drag_rotate_delta...\n");
389 DRAG_SAVING(rotate, transform);
391 static void drag_sizeadj_delta(double dx, double dy) {
392 const double min_eyes_apart= -200;
393 sizeadj.scale *= pow(3.0, dy);
395 sizeadj.eyes_apart += dx * wmaxdim;
396 if (sizeadj.eyes_apart < min_eyes_apart) sizeadj.eyes_apart= min_eyes_apart;
397 printf("sizeadj eyes_apart %g\n", sizeadj.eyes_apart);
400 DRAG_SAVING(sizeadj, sizeadj);
402 static const Drag *drag= &drag_none;
404 static int drag_last_x, drag_last_y;
406 static void drag_position(int x, int y) {
407 drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
408 (y - drag_last_y) * 1.0 / wmaxdim);
413 static void event_button(XButtonEvent *e) {
414 if (e->window != window || !e->same_screen) return;
415 if (e->type == ButtonPress) {
416 if (e->state || drag != &drag_none) {
417 printf("drag=%s press state=0x%lx abandon\n",
418 drag->name, (unsigned long)e->state);
424 case Button1: drag= &drag_rotate; break;
425 case Button2: drag= &drag_sizeadj; break;
426 case Button3: drag= &drag_3d; break;
427 default: printf("unknown drag start %d\n", e->button);
429 printf("drag=%s press button=%lu start %d,%d\n",
430 drag->name, (unsigned long)e->button, e->x, e->y);
435 if (e->type == ButtonRelease) {
436 printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
437 drag_position(e->x, e->y);
443 static void event_motion(int x, int y) {
444 printf("drag=%s motion %d,%d\n", drag->name, x, y);
448 static void event_config(XConfigureEvent *e) {
449 if (e->width == wwidth && e->height == wheight)
452 wwidth= e->width; wheight= e->height;
453 wmaxdim= wwidth > wheight ? wwidth : wheight;
454 wmindim= wwidth < wheight ? wwidth : wheight;
456 XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
457 for (currentbuffer=0; currentbuffer<2; currentbuffer++)
458 XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
464 int main(int argc, const char *const *argv) {
467 int motion_deferred=0, motion_x=-1, motion_y=-1;
469 if (argc != 2 || argv[1][0]=='-') {
470 fputs("need filename\n",stderr); exit(8);
472 input_filename= argv[1];
475 K transform[k][k]= 1.0;
479 XMapWindow(display,window);
481 if (motion_deferred) {
482 int r= XCheckMaskEvent(display,~0UL,&event);
484 event_motion(motion_x, motion_y);
489 XNextEvent(display,&event);
491 switch (event.type) {
494 case ButtonRelease: event_button(&event.xbutton); break;
496 case ConfigureNotify: event_config(&event.xconfigure); break;
499 motion_x= event.xmotion.x;
500 motion_y= event.xmotion.y;
505 printf("unknown event type %u 0x%x\n", event.type,event.type);