chiark / gitweb /
before abolish oee_edgemap in favour of something weirder
[moebius2.git] / project.c
index 0921edf..266237f 100644 (file)
--- a/project.c
+++ b/project.c
@@ -15,11 +15,8 @@ static Triangle trisbuffer[MAXTRIS], *displaylist[MAXTRIS];
 static int ntris;
 static Vertices conformation;
 
-#define STATIC_MATRIX(x)                                       \
-  static double x[D3][D3];                                     \
-  static gsl_matrix x##_gsl= { D3,D3,D3,&x[0][0] };
-
-STATIC_MATRIX(transform);
+static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
+static GSL_MATRIX(transform);
 
 const char *input_filename;
 
@@ -34,9 +31,17 @@ static void read_input(void) {
 }
 
 static void transform_coordinates(void) {
-  int v;
+  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];
   }
 }
 
@@ -71,10 +76,14 @@ static void generate_display_list(void) {
 }    
 
 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= ta->vertex[0][2];
-  double zb= tb->vertex[0][2];
+  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;
 }
@@ -87,22 +96,21 @@ static void sort_display_list(void) {
 
 #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 int wwidth=WSZ, wheight=WSZ, wmaxdim=WSZ, currentbuffer;
-static int ncut;
-
-static double scale= WSZ * 0.3;
-static double eye_z= -10, eye_x= 0;
-static double cut_z= -9;
 
-static void xdie(int l, const char *str) {
-  fprintf(stderr,"X library call failed, line %d: %s\n", l, str);
-}
+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;
 
-#define XA(w) ((w) ? (void)0 : xdie(__LINE__, #w))
+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];
@@ -117,67 +125,99 @@ static void drawtriangle(const Triangle *t) {
     if (z < cut_z) { ncut++; return; }
     
     double zezezp= eye_z / (eye_z - z);
-    points[i].x= scale * (zezezp * (x - eye_x) + eye_x) + WSZ/2;
-    points[i].y= scale * (zezezp *  y                 ) + WSZ/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=
   ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask;
 
-static void display_prepare(void) {
+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;
-  int screen, depth;
-  XVisualInfo vinfo;
   XSizeHints hints;
   
   XA( display= XOpenDisplay(0) );
-  screen= DefaultScreen(display);
-  depth= DefaultDepth(display,screen);
-  XA( XMatchVisualInfo(display,screen,depth, TrueColor,&vinfo) );
+  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,depth,
-                           InputOutput, vinfo.visual,
+                           0,0, wwidth,wheight, 0,x11depth,
+                           InputOutput, visinfo.visual,
                            CWEventMask, &wa) );
 
   hints.flags= USPosition;
   hints.x= 10;
   hints.y= 10;
   XSetWMNormalHints(display,window,&hints);
-  
-  for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
-    XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,depth) );
-    doublebuffers[currentbuffer]= pixmap;
-  }
 
-  gcv.function= GXcopy;
-  gcv.plane_mask= ~0UL;
-  gcv.foreground= WhitePixel(display,screen);
-  linegc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask|GCForeground, &gcv);
-  
-  gcv.function= GXclear;
-  gcv.plane_mask= ~0UL;
-  fillgc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask, &gcv);
+  mkpixmaps();
 
-  currentbuffer= 0;
+  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 (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) );
@@ -185,6 +225,7 @@ static void display_conformation(void) {
 }
 
 static void show(void) {
+  scale_wmindim= sizeadj_scale * wmindim;
   read_input();
   transform_coordinates();
   generate_display_list();
@@ -224,41 +265,146 @@ 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 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_MATRIX(rotateby);
-  STATIC_MATRIX(qr);
-  static double tau[D3];
-  static gsl_vector tau_gsl= { D3,1,&tau[0] };
+  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;
   
-  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);
   
-  gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
-                &rotateby_gsl,&transform_gsl, 0.0,&qr_gsl);
-  gsl_linalg_QR_decomp(&qr_gsl, &tau_gsl);
-  gsl_linalg_QR_unpack(&qr_gsl, &tau_gsl,
-                      &transform_gsl, &rotateby_gsl /*dummy*/);
+  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_scale_delta(double dx, double dy) {
-  scale *= pow(3.0, dy);
+static void drag_sizeadj_delta(double dx, double dy) {
+  sizeadj_scale *= pow(3.0, -dy);
   show();
 }
-DRAG_SAVING(scale, scale);
+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;
 
@@ -283,7 +429,8 @@ static void event_button(XButtonEvent *e) {
     }
     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",
@@ -305,10 +452,25 @@ static void event_motion(int x, int y) {
   drag_position(x,y);
 }
 
-static void event_config(XConfigureEvent *e) { printf("configure\n"); }
+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]=='-') {
@@ -317,6 +479,7 @@ int main(int argc, const char *const *argv) {
   input_filename= argv[1];
 
   read_input();
+  K transform[k][k]= 1.0;
   display_prepare();
   show();