printf("]\n");
}
static void pvector(const char *n, double v[D3]) {
- pvectorcore(m,v);
+ pvectorcore(n,v);
putchar('\n');
}
static void pmatrix(const char *n, double m[D3][D3]) {
int j;
- for (j=0; j<D3; j++) { pvector(n,m[j]); n=""; }
+ 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);
+
+ 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);
GSL_VECTOR(sigma);
GSL_VECTOR(tau);
- rotateby[2][0]= dx;
- rotateby[2][1]= dy;
- PMATRIX(rotateby);
-
- GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
- &rotateby_gsl,&transform_gsl,
- 0.0, &skew_gsl) );
-
- /* Now we want to normalise skew, the result becomes new transform */
-
- pmatrix("skew", skew);
-
/* We use notation from Wikipedia Polar_decomposition
* Wikipedia's W is GSL's U
* Wikipedia's Sigma is GSL's S
case ConfigureNotify: event_config(&event.xconfigure); break;
case MotionNotify:
-#if 0
motion_x= event.xmotion.x;
motion_y= event.xmotion.y;
motion_deferred= 1;
-#endif
continue;
default: