X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=interpolate.c;h=f162fc040981b12f97864ec8169792db3b64cf3f;hp=6f519e414a2d35664a443babb6b665f1b3cae4cd;hb=701de0c6ae90ac6cd5282ed882ad91b11b098b2c;hpb=97ef96dd318b684d829ac950403582788df427d8

diff --git a/interpolate.c b/interpolate.c
index 6f519e4..f162fc0 100644
--- a/interpolate.c
+++ b/interpolate.c
@@ -5,46 +5,72 @@
 #include "mgraph.h"
 
 /* filled in by characterise_input: */
-static int shift, oldxbits, oldybits, oldx, oldy, oldsz, inc;
+static int oXBITS,oX,oY,oN,inc;
 
 /* filled in by read_input: */
 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;
+  int r, shift, oN_calc;
   
   r= fstat(0,&stab);  if (r) diee("fstat input to find length");
 
-  if (!stab.st_size || stab.st_size % sizeof(double) ||
+  if (!stab.st_size || stab.st_size % (sizeof(double)*D3) ||
       stab.st_size > INT_MAX)
-    fail("input file is not reasonable whole number of doubles\n");
-
-  oldsz= stab.st_size / sizeof(double);
-  for (shift=1;
-       shift > XBITS+1 && shift > YBITS+1;
-       shift++) {
-    oldxbits= XBITS-1;
-    oldybits= YBITS-1;
-    oldx=  1<<oldxbits;
-    oldy= (1<<oldybits)-1;
-    if (oldx*oldy == oldsz) goto found;
+    fail("input file is not reasonable whole number of vertices\n");
+
+  printf("XBITS=%d XY=%d*%d=%d", XBITS, X,Y,N);
+
+  oN= stab.st_size / (sizeof(double)*D3);
+
+  for (shift=1, inc=2;
+       shift<XBITS-1;
+       shift++, inc<<=1) {
+    printf("\n shift=%d inc=%d ",shift,inc);
+
+    if (X % inc) continue;
+    oX= X/inc; printf("oX=%d ",oX);
+
+    if ((Y-1) % inc) continue;
+    oY= (Y-1)/inc + 1; printf("oY=%d ",oY);
+
+    oN_calc= oX*oY;
+    printf("oN=%d",oN_calc);
+    if (oN_calc == oN) goto found;
   }
-  fail("input file size cannot be interpolated to target file size\n");
+  fail("\ninput size cannot be reconciled\n");
 
  found:
-  inc= 1<<shift;
+  oXBITS= XBITS-shift;
+  printf("oXBITS=%d\n", oXBITS);
 }
 
 static void read_input(void) {
-  int x,y;
+  int x,y, ox,oy, v,ov;
   
-  for (y=0; y<Y; y+=inc)
-    for (x=0; x<X; x+=inc) {
+  for (oy=0; oy<oY; oy++) {
+    y= oy*inc;
+    printf("y=%2d>%2d", oy,y);
+    for (ox=0; ox<oX; ox++) {
+      x= (ox - (oy>>1))*inc + (y>>1) + X*2;
+      while (x>=X) { y= (Y-1)-y; x-=X; }
       errno= 0;
-      if (fread(all.a[(y << YSHIFT) | x], sizeof(double), D3, stdin) != D3)
-	diee("read input");
+      ov= (oy << oXBITS) | ox;
+      v= (y << YSHIFT) | x;
+      printf(" 0%02o->0x%02x", ov, v);
+      if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
+	diee("\nread input");
+      computed_count[v]= D3*2;
     }
+    putchar('\n');
+  }
 }
 
   /* We use GSL's interpolation functions.  Each xa array is simple
@@ -63,74 +89,234 @@ static void read_input(void) {
    *
    */
 
-#define NEXTV (v= EDGE_END2(v,direction))
+typedef struct {
+  int v, e;
+} Traverse;
+
+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);
+    assert(EDGE_END2(v2,e2) == t->v);
+    t->e= (e2+3) % V6;
+  }
+  t->v= v2;
+}
+
+static void interpolate_lin4point(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];
+
+  if (inc != 2) fail("cannot do >2x lin4point interpolation\n");
+  
+  for (eqr=1; eqr!=4; eqr+=5, eqr%=V6) { /* each old edge exactly once */
+    printf("eqr=%d\n",eqr);
+    for (yq=0; yq<Y; yq+=inc) {
+      printf(" 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;
+	
+	printf(" 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];
 
-static void interpolate_line(int startvertex,
-			     int direction /* edge number */,
-			     int norig) {
-  double xa[norig+1], ya[D3][norig+1], n;
+	K {
+	  const double alpha= 6;
+	  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]= (pqtarg[k] + srtarg[k] + alpha * all.a[vm][k])
+	    / (2 + alpha);
+	  note_computed(vm,k);
+	}
+      }
+      putchar('\n');
+    }
+  }
+}
+
+static void interpolate_gsl_line
+  (int startvertex, int direction /* edge number */, int nmax,
+   const gsl_interp_type *aperiodic, const gsl_interp_type *periodic) {
+  double xa[nmax], ya[D3][nmax];
+  int vnew[nmax], xnew[nmax];
   const gsl_interp_type *it;
   gsl_interp *interp;
   gsl_interp_accel *accel;
-  int i, v, k;
+  int i, k, nold, nnew, x, skipped, realstart;
+  Traverse traverse;
+
+  printf("interpolate_line 0x%02x,%d nmax=%d ",
+	  startvertex,direction,nmax);
+
+  traverse.v=startvertex, traverse.e=direction;
+  skipped= 0;
+  while (!computed_count[traverse.v]) {
+    printf("-%02x~",traverse.v);
+    traverse_next(&traverse);
+    skipped++;
+    assert(skipped < nmax);
+  }
+  realstart= traverse.v;
 
-  for (i=0, v=startvertex;
+  nold=0, nnew=0;
+  for (x=0;
        ;
-       i++, NEXTV, assert(v>=0), NEXTV) {
-    if (v<0)
+       x++, traverse_next(&traverse)) {
+    printf("-%02x",traverse.v);
+    if (traverse.v<0) {
+      putchar('$');
+      it= aperiodic;
+      assert(!skipped);
       break;
-    assert(i <= norig);
-    xa[i]= i;
-    K ya[k][i]= all.a[v][k];
-    if (v==startvertex)
+    }
+    if (computed_count[traverse.v]) {
+      assert(nold < nmax);
+      xa[nold]= x;
+      K ya[k][nold]= all.a[traverse.v][k];
+      nold++;
+      putchar('*');
+    } else {
+      assert(nnew < nmax);
+      xnew[nnew]= x;
+      vnew[nnew]= traverse.v;
+      nnew++;
+      putchar('#');
+    }
+    if (x && traverse.v==realstart) {
+      putchar('@');
+      it= periodic;
       break;
+    }
   }
-  n= i;
-  it= v>=0 ? gsl_interp_cspline_periodic : gsl_interp_cspline;
+
+  printf("  n=%d->%d", nold,nnew);
+  if (!nnew) { printf(" nothing\n"); return; }
+
+  assert(nold);
+  printf(" x=%g..%g->%d..%d", xa[0],xa[nold-1], xnew[0],xnew[nnew-1]);
 
   K {
-    interp= gsl_interp_alloc(it,n); if (!interp) diee("gsl_interp_alloc");
+    printf("\n  k=%d ", k);
+
+    interp= gsl_interp_alloc(it,nold); if (!interp) diee("gsl_interp_alloc");
     accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
-    GA( gsl_interp_init(interp,xa,ya[k],n) );
+    GA( gsl_interp_init(interp,xa,ya[k],nold) );
 
-    for (i=0, v=startvertex;
-	 i<n-1;
-	 i++, NEXTV) {
-      assert(v>=0); NEXTV; assert(v>=0);
-      GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel, &all.a[v][k]) );
+    for (i=0; i<nnew; i++) {
+      int v= vnew[i];
+      x= xnew[i];
+      printf(" %02x(%d)",v,x);
+      GA( gsl_interp_eval_e(interp,xa,ya[k], x, accel, &all.a[v][k]) );
+      note_computed(v,k);
     }
     
     gsl_interp_accel_free(accel);
     gsl_interp_free(interp);
   }
+  printf("    done\n");
 }
 	
-static void interpolate(void) {
+static void interpolate_gsl(const gsl_interp_type *aperiodic,
+			    const gsl_interp_type *periodic) {
   int x,y;
   
-  if (shift!=1) fail("only interpolation by factor of 2 supported\n");
-
-  for (y=0; y<Y; y+=inc) {
-    interpolate_line(y<<YSHIFT, 0, oldx);
-  }
   for (x=0; x<X; x+=inc) {
-    interpolate_line(x, 5, oldy);
-    interpolate_line(x, 4, oldy);
+    interpolate_gsl_line(                x, 5, Y, aperiodic, periodic);
+    interpolate_gsl_line((Y-1)<<YSHIFT | x, 1, Y, aperiodic, periodic);
+  }
+  for (y=0; y<(Y+1)/2; y++) {
+    interpolate_gsl_line(y<<YSHIFT | ((y>>1)&1), 0, X*2+1,
+			 aperiodic, periodic);
   }
 }
 
-static void write_output(void) {
-  if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
+static void write_output(const char *fn) {
+  FILE *f;
+  char *fn_new;
+  int x, y, v, c, bad=0;
+  
+  for (y=0; y<Y; y++) {
+    printf("checking y=%d ", y);
+    for (x=0; x<X; x++) {
+      v= y<<YSHIFT | x;
+      printf(" %02x",v);
+      c= computed_count[v];
+      if (c==D3) putchar('#');
+      else if (c==D3*2) putchar('*');
+      else { printf("!%d",c); bad++; }
+    }
+    putchar('\n');
+  }
+  assert(!bad);
+
+  if (asprintf(&fn_new,"%s.new",fn) <= 0) diee("asprintf");
+  f= fopen(fn_new,"wb");  if (!f) diee("open new output file");
+  if (fwrite(all.a,sizeof(all.a),1,f) != 1) diee("fwrite");
+  if (fclose(f)) diee("fclose output");
+  if (rename(fn_new,fn)) diee("install new output");
+  printf("interpolated\n");
 }
 
 int main(int argc, const char **argv) {
-  if (argc!=1 || isatty(0) || isatty(1))
-    { fputs("usage: interpolate <INPUT.cfm >OUTPUT.cfm\n",stderr); exit(8); }
+  if (argc!=4 ||
+      strncmp(argv[1],"-a",2) || strlen(argv[1])!=3 ||
+      argv[2][0]=='-' ||
+      strncmp(argv[3],"-o",2)) {
+    fputs("usage: interpolate -aK INPUT.cfm -oOUTPUT.cfm\n",stderr);
+    exit(8);
+  }
+
+  if (!freopen(argv[2],"rb",stdin)) diee("open input");
 
+  mgraph_prepare();
   characterise_input();
   read_input();
-  interpolate();
-  write_output();
-  if (fclose(stdout)) diee("fclose stdout");
+
+  switch (argv[1][2]) {
+  case '4':  interpolate_lin4point();                       break;
+  case 'a':  interpolate_gsl(gsl_interp_akima,
+			     gsl_interp_akima_periodic);    break;
+  case 'c':  interpolate_gsl(gsl_interp_cspline,
+			     gsl_interp_cspline_periodic);  break;
+  default:
+    fail("unknown interpolation algorithm\n");
+  }
+  
+  write_output(argv[3]+2);
   return 0;
 }