chiark / gitweb /
linear interpolation works; view presets
[moebius2.git] / interpolate.c
index f5d7c1b..b705bcf 100644 (file)
 static struct Vertices all;
 static int computed_count[N];
 
 static struct Vertices all;
 static int computed_count[N];
 
+static void note_computed(int v, int k) {
+  assert(computed_count[v] == k);
+  computed_count[v]++;
+}
+
 static void characterise_input(void) {
   struct stat stab;
   int r;
 static void characterise_input(void) {
   struct stat stab;
   int r;
@@ -36,9 +41,12 @@ static void characterise_input(void) {
 static void read_input(void) {
   int x,y, ox,oy, v,ov;
   
 static void read_input(void) {
   int x,y, ox,oy, v,ov;
   
-  for (oy=y=0; y<Y; oy++, y+=INC) {
+  for (oy=0; oy<OY; oy++) {
+    y= oy*2;
     fprintf(stderr, "y=%2d>%2d", oy,y);
     fprintf(stderr, "y=%2d>%2d", oy,y);
-    for (ox=x=0; x<X; ox++, x+=INC) {
+    for (ox=0; ox<OX; ox++) {
+      x= ox*2 + (oy&1);
+      if (x>=X) { y= (Y-1)-y; x-=X; }
       errno= 0;
       ov= (oy << OXBITS) | ox;
       v= (y << YSHIFT) | x;
       errno= 0;
       ov= (oy << OXBITS) | ox;
       v= (y << YSHIFT) | x;
@@ -73,6 +81,10 @@ typedef struct {
 
 static void traverse_next(Traverse *t) {
   int v2;
 
 static void traverse_next(Traverse *t) {
   int v2;
+
+  if (t->v<0)
+    return;
+  
   v2= EDGE_END2(t->v, t->e);
   if (v2>=0) {
     int e2= edge_reverse(t->v, t->e);
   v2= EDGE_END2(t->v, t->e);
   if (v2>=0) {
     int e2= edge_reverse(t->v, t->e);
@@ -81,7 +93,66 @@ static void traverse_next(Traverse *t) {
   }
   t->v= v2;
 }
   }
   t->v= v2;
 }
-     
+
+static void interpolate(void) {
+  /* four points P Q R S, although P and S may be missing
+   * interpolate in QR finding M. */
+  int xq,yq,eqr, vp,vq,vm,vr,vs, k;
+  double pqtarg[D3], srtarg[D3];
+
+  for (eqr=1; eqr!=4; eqr+=5, eqr%=V6) { /* each old edge exactly once */
+    fprintf(stderr,"eqr=%d\n",eqr);
+    for (yq=0; yq<Y; yq+=INC) {
+      fprintf(stderr," yq=%2d ",yq);
+      for (xq=((yq>>1)&1); xq<X; xq+=INC) {
+       vq= yq << YSHIFT | xq;
+       Traverse trav; trav.v=vq; trav.e=eqr;
+
+       traverse_next(&trav);
+       vm= trav.v;
+       if (vm<0) continue;
+
+       traverse_next(&trav);
+       vr= trav.v;
+       assert(vr>=0);
+
+       traverse_next(&trav);
+       traverse_next(&trav);
+       vs= trav.v;
+
+       trav.v= vq; trav.e= EDGE_OPPOSITE(eqr);
+       traverse_next(&trav);
+       traverse_next(&trav);
+       vp= trav.v;
+       
+       fprintf(stderr," 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
+               vp&0xff,vq,vm,vr,vs&0xff);
+
+       if (vp>=0)
+         K pqtarg[k]= all.a[vq][k]*2 - all.a[vp][k];
+       else
+         K pqtarg[k]= all.a[vr][k];
+
+       if (vs>=0)
+         K srtarg[k]= all.a[vr][k]*2 - all.a[vs][k];
+       else
+         K srtarg[k]= all.a[vq][k];
+
+       K {
+         all.a[vm][k]= 0.5 * (all.a[vq][k] + all.a[vr][k]);
+//       pqtarg[k]= 0.5 * (pqtarg[k] + all.a[vm][k]);
+//       srtarg[k]= 0.5 * (srtarg[k] + all.a[vm][k]);
+//       all.a[vm][k]= 0.5 * (pqtarg[k] + srtarg[k]);
+         note_computed(vm,k);
+       }
+      }
+      fputc('\n',stderr);
+    }
+  }
+}
+
+#if 0
+
 static void interpolate_line(int startvertex,
                             int direction /* edge number */,
                             int nmax) {
 static void interpolate_line(int startvertex,
                             int direction /* edge number */,
                             int nmax) {
@@ -136,10 +207,9 @@ static void interpolate_line(int startvertex,
         i++, NEXTV) {
       assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
       fputc('#',stderr);
         i++, NEXTV) {
       assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
       fputc('#',stderr);
-      assert(computed_count[traverse.v] == k);
       GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
                            &all.a[traverse.v][k]) );
       GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
                            &all.a[traverse.v][k]) );
-      computed_count[traverse.v]++;
+      note_computed(traverse.v, k);
     }
     
     gsl_interp_accel_free(accel);
     }
     
     gsl_interp_accel_free(accel);
@@ -159,6 +229,7 @@ static void interpolate(void) {
     interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
   }
 }
     interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
   }
 }
+#endif
 
 static void write_output(void) {
   int x, y, v, c, bad=0;
 
 static void write_output(void) {
   int x, y, v, c, bad=0;
@@ -170,7 +241,7 @@ static void write_output(void) {
       c= computed_count[v];
       if (c==D3) fputc('#',stderr);
       else if (c==-1) fputc('*',stderr);
       c= computed_count[v];
       if (c==D3) fputc('#',stderr);
       else if (c==-1) fputc('*',stderr);
-      else { fputc('!',stderr); bad++; }
+      else { fprintf(stderr,"!%d",c); bad++; }
     }
     fputc('\n',stderr);
   }
     }
     fputc('\n',stderr);
   }