--- /dev/null
+/*
+ * Displays a conformation
+ */
+
+#include <X11/Xlib.h>
+#include <X11/Xutil.h>
+
+#include "mgraph.h"
+
+#define MAXTRIS (N*2)
+
+typedef struct { double vertex[3][D3]; } Triangle;
+
+static Triangle trisbuffer[MAXTRIS], *displaylist[MAXTRIS];
+static int ntris;
+static Vertices conformation;
+
+static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
+static GSL_MATRIX(transform);
+
+const char *input_filename;
+
+static void read_input(void) {
+ FILE *f;
+ int r;
+
+ f= fopen(input_filename, "rb"); if (!f) diee("input file");
+ errno= 0;
+ r= fread(&conformation,sizeof(conformation),1,f); if (r!=1) diee("fread");
+ fclose(f);
+}
+
+static void transform_coordinates(void) {
+ double result[D3];
+ GSL_VECTOR(result);
+ gsl_vector input_gsl= { D3,1 };
+
+ int v, k;
+
+ FOR_VERTEX(v) {
+ input_gsl.data= &conformation[v][0];
+ GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
+ 0.0, &result_gsl) );
+ K conformation[v][k]= result[k];
+ }
+}
+
+static void addtriangle(int va, int vb, int vc) {
+ Triangle *t= &trisbuffer[ntris];
+ int k;
+
+ assert(ntris < MAXTRIS);
+ K {
+ t->vertex[0][k]= conformation[va][k];
+ t->vertex[1][k]= conformation[vb][k];
+ t->vertex[2][k]= conformation[vc][k];
+ }
+ displaylist[ntris++]= t;
+}
+
+static void generate_display_list(void) {
+ int vb, ve[3], e;
+
+ ntris= 0;
+ FOR_VERTEX(vb) {
+ /* We use the two triangles in the parallelogram vb, vb+e0, vb+e1, vb+e2.
+ * We go round each triangle clockwise (although our surface is non-
+ * orientable so it shouldn't matter).
+ */
+ for (e=0; e<3; e++) ve[e]= EDGE_END2(vb,e);
+ if (ve[1]>=0) {
+ if (ve[0]>=0) addtriangle(vb,ve[0],ve[1]);
+ if (ve[2]>=0) addtriangle(vb,ve[1],ve[2]);
+ }
+ }
+}
+
+static int dl_compare(const void *tav, const void *tbv) {
+ int i;
+ const Triangle *const *tap= tav, *ta= *tap;
+ const Triangle *const *tbp= tbp, *tb= *tbp;
+ double za=0, zb=0;
+ for (i=0; i<3; i++) {
+ za += ta->vertex[i][2];
+ zb += tb->vertex[i][2];
+ }
+ return za > zb ? -1 :
+ za < zb ? +1 : 0;
+}
+
+static void sort_display_list(void) {
+ qsort(displaylist, ntris, sizeof(*displaylist), dl_compare);
+}
+
+/*---------- X stuff ----------*/
+
+#define WSZ 400
+
+typedef struct { GC fillgc, linegc; } DrawingMode;
+
+static Display *display;
+static Pixmap pixmap, doublebuffers[2];
+static Window window;
+
+static DrawingMode dmred, dmblue, dmwhite;
+static const DrawingMode *dmcurrent;
+static int wwidth=WSZ, wheight=WSZ, wmindim=WSZ, wmaxdim=WSZ;
+static int ncut, currentbuffer, x11depth, x11screen;
+XVisualInfo visinfo;
+
+static double sizeadj_scale= 0.3, eyes_apart, scale_wmindim;
+static double eye_z= -10, eye_x=0;
+static double cut_z= -9;
+
+static void drawtriangle(const Triangle *t) {
+ XPoint points[4];
+ int i;
+
+ for (i=0; i<3; i++) {
+ double *v= t->vertex[i];
+ double x= v[0];
+ double y= v[1];
+ double z= v[2];
+
+ if (z < cut_z) { ncut++; return; }
+
+ double zezezp= eye_z / (eye_z - z);
+ points[i].x= scale_wmindim * (zezezp * (x - eye_x) + eye_x) + wwidth/2;
+ points[i].y= scale_wmindim * (zezezp * y ) + wheight/2;
+ }
+ points[3]= points[0];
+
+ XA( XFillPolygon(display,pixmap, dmcurrent->fillgc,
+ points,3,Convex,CoordModeOrigin) );
+ XA( XDrawLines(display,pixmap, dmcurrent->linegc,
+ points, 4,CoordModeOrigin) );
+}
+
+static const unsigned long core_event_mask=
+ ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask;
+
+static void mkpixmaps(void) {
+ for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
+ XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,x11depth) );
+ doublebuffers[currentbuffer]= pixmap;
+ }
+ currentbuffer= 0;
+}
+
+static void mkgcs(DrawingMode *dm, unsigned long planes) {
+ XGCValues gcv;
+
+ gcv.function= GXcopy;
+ gcv.foreground= WhitePixel(display,x11screen);
+ gcv.plane_mask= planes;
+ dm->linegc= XCreateGC(display,pixmap,
+ GCFunction|GCForeground|GCPlaneMask,
+ &gcv);
+
+ gcv.function= GXclear;
+ dm->fillgc= XCreateGC(display,pixmap,
+ GCFunction|GCPlaneMask,
+ &gcv);
+}
+
+static void display_prepare(void) {
+ XSetWindowAttributes wa;
+ XSizeHints hints;
+
+ XA( display= XOpenDisplay(0) );
+ x11screen= DefaultScreen(display);
+ x11depth= DefaultDepth(display,x11screen);
+ XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&visinfo) );
+
+ wa.event_mask= core_event_mask;
+ XA( window= XCreateWindow(display, DefaultRootWindow(display),
+ 0,0, wwidth,wheight, 0,x11depth,
+ InputOutput, visinfo.visual,
+ CWEventMask, &wa) );
+
+ hints.flags= USPosition;
+ hints.x= 10;
+ hints.y= 10;
+ XSetWMNormalHints(display,window,&hints);
+
+ mkpixmaps();
+
+ mkgcs(&dmwhite, AllPlanes);
+ mkgcs(&dmblue, visinfo.blue_mask);
+ mkgcs(&dmred, visinfo.red_mask);
+}
+
+static void drawtriangles(const DrawingMode *dm) {
+ Triangle *const *t;
+ int i;
+
+ dmcurrent= dm;
+ for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
+ drawtriangle(*t);
+}
+
+static void display_conformation(void) {
+ pixmap= doublebuffers[currentbuffer];
+
+ XA( XFillRectangle(display,pixmap,dmwhite.fillgc,0,0,wwidth,wheight) );
+
+ if (eyes_apart > 0) {
+ const double preferred=0.05, beyond=0.07;
+
+ eye_x= eyes_apart < preferred ? eyes_apart :
+ eyes_apart < beyond ? preferred :
+ eyes_apart - (beyond - preferred);
+ eye_x /= sizeadj_scale;
+ drawtriangles(&dmblue);
+ eye_x= -eye_x;
+ drawtriangles(&dmred);
+ } else {
+ drawtriangles(&dmwhite);
+ printf("shown, %d/%d triangles cut\n", ncut, ntris);
+ }
+
+ XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
+ XA( XClearWindow(display,window) );
+ currentbuffer= !currentbuffer;
+}
+
+static void show(void) {
+ scale_wmindim= sizeadj_scale * wmindim;
+ read_input();
+ transform_coordinates();
+ generate_display_list();
+ sort_display_list();
+ display_conformation();
+}
+
+typedef struct {
+ const char *name;
+ void (*start)(void);
+ void (*delta)(double dx, double dy);
+ void (*conclude)(void);
+ void (*abandon)(void);
+} Drag;
+
+#define DRAG(x) \
+ static const Drag drag_##x= { \
+ #x, drag_##x##_start, drag_##x##_delta, \
+ drag_##x##_conclude, drag_##x##_abandon \
+ }
+
+#define DRAG_SAVING(x, thing) \
+ static typeof(thing) original_##thing; \
+ static void drag_##x##_start(void) { \
+ memcpy(&original_##thing, &thing, sizeof(thing)); \
+ } \
+ static void drag_##x##_conclude(void) { } \
+ static void drag_##x##_abandon(void) { \
+ memcpy(&thing, &original_##thing, sizeof(thing)); \
+ show(); \
+ } \
+ DRAG(x)
+
+static void drag_none_start(void) { }
+static void drag_none_delta(double dx, double dy) { }
+static void drag_none_conclude(void) { }
+static void drag_none_abandon(void) { }
+DRAG(none);
+
+static void pvectorcore(const char *n, double v[D3]) {
+ int k;
+ printf("%10s [ ",n);
+ K printf("%# 10.10f ",v[k]);
+ printf("]\n");
+}
+static void pvector(const char *n, double v[D3]) {
+ pvectorcore(n,v);
+ putchar('\n');
+}
+static void pmatrix(const char *n, double m[D3][D3]) {
+ int j;
+ for (j=0; j<D3; j++) { pvectorcore(n,m[j]); n=""; }
+ putchar('\n');
+}
+#define PMATRIX(x) pmatrix(#x,x);
+
+static void drag_rotate_delta(double dx, double dy) {
+ /* We multiple our transformation matrix by a matrix:
+ *
+ * If we just had y movement, we would rotate about x axis:
+ * rotation X = [ 1 0 0 ]
+ * [ 0 cy sy ]
+ * [ 0 -sy cy ]
+ * where cy and sy are sin and cos of y rotation
+ *
+ * But we should pre-rotate this by a rotation about the z axis
+ * to get it to the right angle (to include x rotation). So
+ * we make cy and sy be cos() and sin(hypot(x,y)) and use
+ * with cr,sr as cos() and sin(atan2(y,y)):
+ *
+ * Ie we would do T' = R^T X R T where
+ * or T' = C T where C = R^T X R and
+ *
+ * adjustment R = [ cr sr 0 ]
+ * [ -sr cr 0 ]
+ * [ 0 0 1 ]
+ */
+
+ double rotx[D3][D3], adjr[D3][D3];
+ GSL_MATRIX(rotx);
+ GSL_MATRIX(adjr);
+
+ static double temp[D3][D3], change[D3][D3];
+ static GSL_MATRIX(temp);
+ static GSL_MATRIX(change);
+
+ double d= hypot(dx,dy);
+ if (d < 1e-6) return;
+
+ double ang= d*2.0;
+
+ double cy= cos(ang);
+ double sy= sin(ang);
+ double cr= -dy / d;
+ double sr= dx / d;
+ printf("\n d=%g cy,sy=%g,%g cr,sr=%g,%g\n\n", d,cy,sy,cr,sr);
+
+ rotx[0][0]= 1; rotx[0][1]= 0; rotx[0][2]= 0;
+ rotx[1][0]= 0; rotx[1][1]= cy; rotx[1][2]= sy;
+ rotx[2][0]= 0; rotx[2][1]= -sy; rotx[2][2]= cy;
+ PMATRIX(rotx);
+
+ adjr[0][0]= cr; adjr[0][1]= sr; adjr[0][2]= 0;
+ adjr[1][0]= -sr; adjr[1][1]= cr; adjr[1][2]= 0;
+ adjr[2][0]= 0; adjr[2][1]= 0; adjr[2][2]= 1;
+ PMATRIX(adjr);
+
+ GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
+ &rotx_gsl,&adjr_gsl,
+ 0.0, &temp_gsl) );
+ PMATRIX(temp);
+
+ GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
+ &adjr_gsl,&temp_gsl,
+ 0.0, &change_gsl) );
+ PMATRIX(change);
+
+ static double skew[D3][D3];
+ static GSL_MATRIX(skew);
+
+ GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
+ &change_gsl,&transform_gsl,
+ 0.0, &skew_gsl) );
+ PMATRIX(skew);
+
+ memcpy(&transform,&skew,sizeof(transform));
+ show();
+ return;
+
+ /* Now we want to normalise skew, the result becomes new transform */
+ double svd_v[D3][D3];
+ GSL_MATRIX(svd_v);
+
+ double sigma[D3], tau[D3];
+ GSL_VECTOR(sigma);
+ GSL_VECTOR(tau);
+
+ /* We use notation from Wikipedia Polar_decomposition
+ * Wikipedia's W is GSL's U
+ * Wikipedia's Sigma is GSL's S
+ * Wikipedia's V is GSL's V
+ * Wikipedia's U is our desired result
+ * Wikipedia which says if the SVD is A = W Sigma V*
+ * then the polar decomposition is A = U P
+ * where P = V Sigma V*
+ * and U = W V*
+ */
+
+ GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
+ pmatrix("W", skew);
+ pvector("Sigma",sigma);
+ pmatrix("V", svd_v);
+
+ /* We only need U, not P. */
+ GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
+ &skew_gsl,&svd_v_gsl,
+ 0.0,&transform_gsl) );
+
+ pmatrix("U", transform);
+
+ printf("drag_rotate_delta...\n");
+ show();
+}
+DRAG_SAVING(rotate, transform);
+
+static void drag_sizeadj_delta(double dx, double dy) {
+ sizeadj_scale *= pow(3.0, -dy);
+ show();
+}
+DRAG_SAVING(sizeadj, sizeadj_scale);
+
+static void drag_3d_delta(double dx, double dy) {
+ const double min_eyes_apart= -0.02;
+ eyes_apart += dx * 0.1;
+ if (eyes_apart < min_eyes_apart) eyes_apart= min_eyes_apart;
+ printf("sizeadj eyes_apart %g\n", eyes_apart);
+ show();
+}
+DRAG_SAVING(3d, eyes_apart);
+
+static const Drag *drag= &drag_none;
+
+static int drag_last_x, drag_last_y;
+
+static void drag_position(int x, int y) {
+ drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
+ (y - drag_last_y) * 1.0 / wmaxdim);
+ drag_last_x= x;
+ drag_last_y= y;
+}
+
+static void event_button(XButtonEvent *e) {
+ if (e->window != window || !e->same_screen) return;
+ if (e->type == ButtonPress) {
+ if (e->state || drag != &drag_none) {
+ printf("drag=%s press state=0x%lx abandon\n",
+ drag->name, (unsigned long)e->state);
+ drag->abandon();
+ drag= &drag_none;
+ return;
+ }
+ switch (e->button) {
+ case Button1: drag= &drag_rotate; break;
+ case Button2: drag= &drag_sizeadj; break;
+ case Button3: drag= &drag_3d; break;
+ default: printf("unknown drag start %d\n", e->button);
+ }
+ printf("drag=%s press button=%lu start %d,%d\n",
+ drag->name, (unsigned long)e->button, e->x, e->y);
+ drag_last_x= e->x;
+ drag_last_y= e->y;
+ drag->start();
+ }
+ if (e->type == ButtonRelease) {
+ printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
+ drag_position(e->x, e->y);
+ drag->conclude();
+ drag= &drag_none;
+ }
+}
+
+static void event_motion(int x, int y) {
+ printf("drag=%s motion %d,%d\n", drag->name, x, y);
+ drag_position(x,y);
+}
+
+static void event_config(XConfigureEvent *e) {
+ if (e->width == wwidth && e->height == wheight)
+ return;
+
+ wwidth= e->width; wheight= e->height;
+ wmaxdim= wwidth > wheight ? wwidth : wheight;
+ wmindim= wwidth < wheight ? wwidth : wheight;
+
+ XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
+ for (currentbuffer=0; currentbuffer<2; currentbuffer++)
+ XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
+
+ mkpixmaps();
+ show();
+}
+
+int main(int argc, const char *const *argv) {
+ XEvent event;
+ int k;
+ int motion_deferred=0, motion_x=-1, motion_y=-1;
+
+ if (argc != 2 || argv[1][0]=='-') {
+ fputs("need filename\n",stderr); exit(8);
+ }
+ input_filename= argv[1];
+
+ read_input();
+ K transform[k][k]= 1.0;
+ display_prepare();
+ show();
+
+ XMapWindow(display,window);
+ for (;;) {
+ if (motion_deferred) {
+ int r= XCheckMaskEvent(display,~0UL,&event);
+ if (!r) {
+ event_motion(motion_x, motion_y);
+ motion_deferred=0;
+ continue;
+ }
+ } else {
+ XNextEvent(display,&event);
+ }
+ switch (event.type) {
+
+ case ButtonPress:
+ case ButtonRelease: event_button(&event.xbutton); break;
+
+ case ConfigureNotify: event_config(&event.xconfigure); break;
+
+ case MotionNotify:
+ motion_x= event.xmotion.x;
+ motion_y= event.xmotion.y;
+ motion_deferred= 1;
+ continue;
+
+ default:
+ printf("unknown event type %u 0x%x\n", event.type,event.type);
+ }
+ }
+}
+
~ian / moebius2.git / blobdiff