before check svd

[moebius2.git] / project.c

diff --git a/project.c b/project.c
index b68e47393bb1b966579879cf0e54c0dd2513c764..764a4dd8fd095d544446c2b357d6dc87a0673d5c 100644 (file)
--- a/project.c
+++ b/project.c
@@ -16,7 +16,7 @@ static int ntris;
 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;
 
@@ -31,16 +31,18 @@ static void read_input(void) {
 }
 
 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) {
@@ -223,6 +225,23 @@ 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(m,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=""; }
+  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 ]
@@ -232,25 +251,52 @@ static void drag_rotate_delta(double dx, double dy) {
    */
 
   static double rotateby[D3][D3]= {{1,0,0},{0,1,0},{-20,-20,1}};
-  STATIC_GSL_MATRIX(rotateby);
+  static GSL_MATRIX(rotateby);
 
-  static double qr[D3][D3];
-  STATIC_GSL_MATRIX(qr);
+  static double skew[D3][D3];
+  static GSL_MATRIX(skew);
 
-  static double tau[D3];
-  STATIC_GSL_VECTOR(tau);
-  
-  int k;
+  double svd_v[D3][D3];
+  GSL_MATRIX(svd_v);
+
+  double sigma[D3], tau[D3];
+  GSL_VECTOR(sigma);
+  GSL_VECTOR(tau);
   
-  K rotateby[k][k]= 1;
-  rotateby[0][2]= dx;
-  rotateby[1][2]= dy;
+  rotateby[2][0]= dx;
+  rotateby[2][1]= dy;
+  PMATRIX(rotateby);
   
   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*/) );
+                    &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
+   *     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();
@@ -359,9 +405,11 @@ int main(int argc, const char *const *argv) {
     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: