static Vertices conformation;
static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
-STATIC_GSL_MATRIX(transform);
+static GSL_MATRIX(transform);
const char *input_filename;
}
static void transform_coordinates(void) {
- /*
- static double result[D3];
- static const gsl_vector result_gsl= { D3,1,&result[0]; }
+ double result[D3];
+ GSL_VECTOR(result);
+ gsl_vector input_gsl= { D3,1 };
- int v;
+ int v, k;
-
FOR_VERTEX(v) {
- GA( gsl_blas_dgemv(CblasNoTrans,
- }*/
+ 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) {
#define WSZ 400
+typedef struct { GC fillgc, linegc; } DrawingMode;
+
static Display *display;
static Pixmap pixmap, doublebuffers[2];
static Window window;
-static GC linegc, fillgc;
+
+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 scale= 0.3;
-static double eye_z= -10, eye_x= 0;
+static struct { double scale, eyes_apart; } sizeadj= { 0.3, 0.0 };
+static double scale_wmindim;
+static double eye_z= -10, eye_x=0;
static double cut_z= -9;
static void drawtriangle(const Triangle *t) {
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[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,fillgc, points,3,Convex,CoordModeOrigin) );
- XA( XDrawLines(display,pixmap,linegc,points, 4,CoordModeOrigin) );
+ 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=
currentbuffer= 0;
}
-static void display_prepare(void) {
+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;
- XVisualInfo vinfo;
XSizeHints hints;
XA( display= XOpenDisplay(0) );
x11screen= DefaultScreen(display);
x11depth= DefaultDepth(display,x11screen);
- XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&vinfo) );
+ 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, vinfo.visual,
+ InputOutput, visinfo.visual,
CWEventMask, &wa) );
hints.flags= USPosition;
mkpixmaps();
- gcv.function= GXcopy;
- gcv.plane_mask= ~0UL;
- gcv.foreground= WhitePixel(display,x11screen);
- linegc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask|GCForeground, &gcv);
-
- gcv.function= GXclear;
- gcv.plane_mask= ~0UL;
- fillgc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask, &gcv);
+ mkgcs(&dmwhite, AllPlanes);
+ mkgcs(&dmblue, visinfo.blue_mask);
+ mkgcs(&dmred, visinfo.red_mask);
}
-static void display_conformation(void) {
+static void drawtriangles(const DrawingMode *dm) {
Triangle *const *t;
int i;
-
- pixmap= doublebuffers[currentbuffer];
- XA( XFillRectangle(display,pixmap,fillgc,0,0,wwidth,wheight) );
+
+ dmcurrent= dm;
for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
drawtriangle(*t);
- printf("shown, %d/%d triangles cut\n", ncut, ntris);
+}
+
+static void display_conformation(void) {
+ pixmap= doublebuffers[currentbuffer];
+
+ XA( XFillRectangle(display,pixmap,dmwhite.fillgc,0,0,wwidth,wheight) );
+
+ if (sizeadj.eyes_apart > 0) {
+ const double preferred=100, beyond=200;
+
+ eye_x= sizeadj.eyes_apart < preferred ? sizeadj.eyes_apart :
+ sizeadj.eyes_apart < beyond ? preferred :
+ sizeadj.eyes_apart - (beyond - preferred);
+ eye_x /= scale_wmindim;
+ 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) );
}
static void show(void) {
+ scale_wmindim= sizeadj.scale * wmindim;
read_input();
transform_coordinates();
generate_display_list();
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 this matrix:
- * [ 1 0 0 ]
- * [ 0 1 0 ]
- * [ dx dy 1 ]
- * and then renormalise.
+ /* 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 ]
*/
- static double rotateby[D3][D3]= {{1,0,0},{0,1,0},{-20,-20,1}};
- STATIC_GSL_MATRIX(rotateby);
+ 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);
- static double qr[D3][D3];
- STATIC_GSL_MATRIX(qr);
+ double d= hypot(dx,dy);
+ if (d < 1e-6) return;
- static double tau[D3];
- STATIC_GSL_VECTOR(tau);
+ double ang= d*2.0;
- int k;
+ 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);
- K rotateby[k][k]= 1;
- rotateby[0][2]= dx;
- rotateby[1][2]= dy;
+ 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,
- &rotateby_gsl,&transform_gsl, 0.0,&qr_gsl) );
- GA( gsl_linalg_QR_decomp(&qr_gsl, &tau_gsl) );
- GA( gsl_linalg_QR_unpack(&qr_gsl, &tau_gsl,
- &transform_gsl, &rotateby_gsl /*dummy*/) );
+ &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_scale_delta(double dx, double dy) {
- scale *= pow(3.0, dy);
+static void drag_sizeadj_delta(double dx, double dy) {
+ const double min_eyes_apart= -200;
+ sizeadj.scale *= pow(3.0, dy);
+
+ sizeadj.eyes_apart += dx * wmaxdim;
+ if (sizeadj.eyes_apart < min_eyes_apart) sizeadj.eyes_apart= min_eyes_apart;
+ printf("sizeadj eyes_apart %g\n", sizeadj.eyes_apart);
show();
}
-DRAG_SAVING(scale, scale);
+DRAG_SAVING(sizeadj, sizeadj);
static const Drag *drag= &drag_none;
}
switch (e->button) {
case Button1: drag= &drag_rotate; break;
- case Button2: drag= &drag_scale; 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",