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 static Display *display;
96 static Pixmap pixmap, doublebuffers[2];
98 static GC linegc, fillgc;
99 static int wwidth=WSZ, wheight=WSZ, wmindim=WSZ, wmaxdim=WSZ;
100 static int ncut, currentbuffer, x11depth, x11screen;
102 static double scale= 0.3;
103 static double eye_z= -10, eye_x= 0;
104 static double cut_z= -9;
106 static void drawtriangle(const Triangle *t) {
110 for (i=0; i<3; i++) {
111 double *v= t->vertex[i];
116 if (z < cut_z) { ncut++; return; }
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;
122 points[3]= points[0];
124 XA( XFillPolygon(display,pixmap,fillgc, points,3,Convex,CoordModeOrigin) );
125 XA( XDrawLines(display,pixmap,linegc,points, 4,CoordModeOrigin) );
128 static const unsigned long core_event_mask=
129 ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask;
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;
139 static void display_prepare(void) {
141 XSetWindowAttributes wa;
145 XA( display= XOpenDisplay(0) );
146 x11screen= DefaultScreen(display);
147 x11depth= DefaultDepth(display,x11screen);
148 XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&vinfo) );
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,
156 hints.flags= USPosition;
159 XSetWMNormalHints(display,window,&hints);
163 gcv.function= GXcopy;
164 gcv.plane_mask= ~0UL;
165 gcv.foreground= WhitePixel(display,x11screen);
166 linegc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask|GCForeground, &gcv);
168 gcv.function= GXclear;
169 gcv.plane_mask= ~0UL;
170 fillgc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask, &gcv);
173 static void display_conformation(void) {
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++)
181 printf("shown, %d/%d triangles cut\n", ncut, ntris);
183 XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
184 XA( XClearWindow(display,window) );
185 currentbuffer= !currentbuffer;
188 static void show(void) {
190 transform_coordinates();
191 generate_display_list();
193 display_conformation();
199 void (*delta)(double dx, double dy);
200 void (*conclude)(void);
201 void (*abandon)(void);
205 static const Drag drag_##x= { \
206 #x, drag_##x##_start, drag_##x##_delta, \
207 drag_##x##_conclude, drag_##x##_abandon \
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)); \
215 static void drag_##x##_conclude(void) { } \
216 static void drag_##x##_abandon(void) { \
217 memcpy(&thing, &original_##thing, sizeof(thing)); \
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) { }
228 static void pvectorcore(const char *n, double v[D3]) {
231 K printf("%# 10.10f ",v[k]);
234 static void pvector(const char *n, double v[D3]) {
238 static void pmatrix(const char *n, double m[D3][D3]) {
240 for (j=0; j<D3; j++) { pvectorcore(n,m[j]); n=""; }
243 #define PMATRIX(x) pmatrix(#x,x);
245 static void drag_rotate_delta(double dx, double dy) {
246 /* We multiple our transformation matrix by a matrix:
248 * If we just had y movement, we would rotate about x axis:
249 * rotation X = [ 1 0 0 ]
252 * where cy and sy are sin and cos of y rotation
254 * But we should pre-rotate this by a rotation about the z axis
255 * to get it to the right angle (to include x rotation). So
256 * we make cy and sy be cos() and sin(hypot(x,y)) and use
257 * with cr,sr as cos() and sin(atan2(y,y)):
259 * Ie we would do T' = R^T X R T where
260 * or T' = C T where C = R^T X R and
262 * adjustment R = [ cr sr 0 ]
267 double rotx[D3][D3], adjr[D3][D3];
271 static double temp[D3][D3], change[D3][D3];
272 static GSL_MATRIX(temp);
273 static GSL_MATRIX(change);
275 double d= hypot(dx,dy);
276 if (d < 1e-6) return;
284 printf("\n d=%g cy,sy=%g,%g cr,sr=%g,%g\n\n", d,cy,sy,cr,sr);
286 rotx[0][0]= 1; rotx[0][1]= 0; rotx[0][2]= 0;
287 rotx[1][0]= 0; rotx[1][1]= cy; rotx[1][2]= sy;
288 rotx[2][0]= 0; rotx[2][1]= -sy; rotx[2][2]= cy;
291 adjr[0][0]= cr; adjr[0][1]= sr; adjr[0][2]= 0;
292 adjr[1][0]= -sr; adjr[1][1]= cr; adjr[1][2]= 0;
293 adjr[2][0]= 0; adjr[2][1]= 0; adjr[2][2]= 1;
296 GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
301 GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
306 static double skew[D3][D3];
307 static GSL_MATRIX(skew);
309 GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
310 &change_gsl,&transform_gsl,
314 memcpy(&transform,&skew,sizeof(transform));
318 /* Now we want to normalise skew, the result becomes new transform */
319 double svd_v[D3][D3];
322 double sigma[D3], tau[D3];
326 /* We use notation from Wikipedia Polar_decomposition
327 * Wikipedia's W is GSL's U
328 * Wikipedia's Sigma is GSL's S
329 * Wikipedia's V is GSL's V
330 * Wikipedia's U is our desired result
331 * Wikipedia which says if the SVD is A = W Sigma V*
332 * then the polar decomposition is A = U P
333 * where P = V Sigma V*
337 GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
339 pvector("Sigma",sigma);
342 /* We only need U, not P. */
343 GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
344 &skew_gsl,&svd_v_gsl,
345 0.0,&transform_gsl) );
347 pmatrix("U", transform);
349 printf("drag_rotate_delta...\n");
352 DRAG_SAVING(rotate, transform);
354 static void drag_scale_delta(double dx, double dy) {
355 scale *= pow(3.0, dy);
358 DRAG_SAVING(scale, scale);
360 static const Drag *drag= &drag_none;
362 static int drag_last_x, drag_last_y;
364 static void drag_position(int x, int y) {
365 drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
366 (y - drag_last_y) * 1.0 / wmaxdim);
371 static void event_button(XButtonEvent *e) {
372 if (e->window != window || !e->same_screen) return;
373 if (e->type == ButtonPress) {
374 if (e->state || drag != &drag_none) {
375 printf("drag=%s press state=0x%lx abandon\n",
376 drag->name, (unsigned long)e->state);
382 case Button1: drag= &drag_rotate; break;
383 case Button2: drag= &drag_scale; break;
384 default: printf("unknown drag start %d\n", e->button);
386 printf("drag=%s press button=%lu start %d,%d\n",
387 drag->name, (unsigned long)e->button, e->x, e->y);
392 if (e->type == ButtonRelease) {
393 printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
394 drag_position(e->x, e->y);
400 static void event_motion(int x, int y) {
401 printf("drag=%s motion %d,%d\n", drag->name, x, y);
405 static void event_config(XConfigureEvent *e) {
406 if (e->width == wwidth && e->height == wheight)
409 wwidth= e->width; wheight= e->height;
410 wmaxdim= wwidth > wheight ? wwidth : wheight;
411 wmindim= wwidth < wheight ? wwidth : wheight;
413 XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
414 for (currentbuffer=0; currentbuffer<2; currentbuffer++)
415 XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
421 int main(int argc, const char *const *argv) {
424 int motion_deferred=0, motion_x=-1, motion_y=-1;
426 if (argc != 2 || argv[1][0]=='-') {
427 fputs("need filename\n",stderr); exit(8);
429 input_filename= argv[1];
432 K transform[k][k]= 1.0;
436 XMapWindow(display,window);
438 if (motion_deferred) {
439 int r= XCheckMaskEvent(display,~0UL,&event);
441 event_motion(motion_x, motion_y);
446 XNextEvent(display,&event);
448 switch (event.type) {
451 case ButtonRelease: event_button(&event.xbutton); break;
453 case ConfigureNotify: event_config(&event.xconfigure); break;
456 motion_x= event.xmotion.x;
457 motion_y= event.xmotion.y;
462 printf("unknown event type %u 0x%x\n", event.type,event.type);