chiark / gitweb /
interpolation more flexible
authorIan Jackson <ian@davenant.relativity.greenend.org.uk>
Sat, 19 Jan 2008 11:10:15 +0000 (11:10 +0000)
committerIan Jackson <ian@davenant.relativity.greenend.org.uk>
Sat, 19 Jan 2008 11:10:15 +0000 (11:10 +0000)
Makefile
interpolate.c

index 79e8edf..bb869a1 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,8 +1,8 @@
 
 
-VIEWDIMS=33 64 125
+VIEWDIMS=33 64 125 246
 TARGETS= minimise primer lumpy.cfm sgtatham.cfm ring.cfm \
-       interpolate-64 \
+       interpolate-64 interpolate-125 interpolate-246 \
        $(addprefix view-, $(VIEWDIMS))
 SGTATHAM=sgtatham
 
index 778e37c..9225b8f 100644 (file)
@@ -4,11 +4,8 @@
 
 #include "mgraph.h"
 
-#define OXBITS (XBITS-1)
-#define OX (1<<OXBITS)
-#define OY (Y/2 + 1)
-#define ON (OX*OY)
-#define INC 2
+/* filled in by characterise_input: */
+static int oXBITS,oX,oY,oN,inc;
 
 /* filled in by read_input: */
 static struct Vertices all;
@@ -21,7 +18,7 @@ static void note_computed(int v, int k) {
 
 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");
 
@@ -29,33 +26,50 @@ static void characterise_input(void) {
       stab.st_size > INT_MAX)
     fail("input file is not reasonable whole number of vertices\n");
 
-  fprintf(stderr,"XBITS=%d XY=%d*%d=%d\n", XBITS, X,Y,N);
-  fprintf(stderr,"OXBITS=%d OXY=%d*%d=%d\n", OXBITS, OX,OY,ON);
-  fprintf(stderr,"sizeof(double)=%d st_size=%lu\n",
-         (int)sizeof(double), (unsigned long)stab.st_size);
+  printf("XBITS=%d XY=%d*%d=%d", XBITS, X,Y,N);
 
-  if (stab.st_size != (ON*sizeof(double)*D3))
-    fail("input file wrong size\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("\ninput size cannot be reconciled\n");
+
+ found:
+  oXBITS= XBITS-shift;
+  printf("oXBITS=%d\n", oXBITS);
 }
 
 static void read_input(void) {
   int x,y, ox,oy, v,ov;
   
-  for (oy=0; oy<OY; oy++) {
+  for (oy=0; oy<oY; oy++) {
     y= oy*2;
-    fprintf(stderr, "y=%2d>%2d", oy,y);
-    for (ox=0; ox<OX; ox++) {
+    printf("y=%2d>%2d", oy,y);
+    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;
+      ov= (oy << oXBITS) | ox;
       v= (y << YSHIFT) | x;
-      fprintf(stderr, " 0%02o->0x%02x", ov, v);
+      printf(" 0%02o->0x%02x", ov, v);
       if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
        diee("\nread input");
       computed_count[v]= -1;
     }
-    fputc('\n',stderr);
+    putchar('\n');
   }
 }
 
@@ -94,19 +108,19 @@ static void traverse_next(Traverse *t) {
   t->v= v2;
 }
 
-#if 0
-
-static void interpolate(void) {
+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 */
-    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) {
+    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;
 
@@ -127,7 +141,7 @@ static void interpolate(void) {
        traverse_next(&trav);
        vp= trav.v;
        
-       fprintf(stderr," 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
+       printf(" 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
                vp&0xff,vq,vm,vr,vs&0xff);
 
        if (vp>=0)
@@ -141,7 +155,7 @@ static void interpolate(void) {
          K srtarg[k]= all.a[vq][k];
 
        K {
-         const double alpha= 3.0;
+         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]);
@@ -150,18 +164,14 @@ static void interpolate(void) {
          note_computed(vm,k);
        }
       }
-      fputc('\n',stderr);
+      putchar('\n');
     }
   }
 }
 
-#endif
-
-#if 1
-
-static void interpolate_line(int startvertex,
-                            int direction /* edge number */,
-                            int nmax) {
+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];
   const gsl_interp_type *it;
   gsl_interp *interp;
@@ -170,9 +180,9 @@ static void interpolate_line(int startvertex,
   Traverse traverse;
 
 #define STARTV (traverse.v=startvertex, traverse.e=direction)
-#define NEXTV (traverse_next(&traverse), fprintf(stderr,"-%02x",traverse.v))
+#define NEXTV (traverse_next(&traverse), printf("-%02x",traverse.v))
 
-  fprintf(stderr,"interpolate_line 0x%02x,%d nmax=%d ",
+  printf("interpolate_line 0x%02x,%d nmax=%d ",
          startvertex,direction,nmax);
 
   for (i=0, STARTV;
@@ -182,27 +192,27 @@ static void interpolate_line(int startvertex,
     assert(i < nmax);
     xa[i]= i;
     K ya[k][i]= all.a[traverse.v][k];
-    fputc('*',stderr);
-    if (i && traverse.v==startvertex) { fputc('@',stderr); break; }
+    putchar('*');
+    if (i && traverse.v==startvertex) { putchar('@'); break; }
     i++;
     NEXTV;
     if (traverse.v<0) break;
     NEXTV;
   }
   if (traverse.v>=0) {
-    it= gsl_interp_akima_periodic;
+    it= periodic;
     nold= i+1;
     nnew= i;
   } else {
-    it= gsl_interp_akima;
+    it= aperiodic;
     nold= i;
     nnew= i-1;
   }
 
-  fprintf(stderr,"  n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
+  printf("  n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
 
   K {
-    fprintf(stderr,"\n  k=%d ", k);
+    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");
@@ -212,7 +222,7 @@ static void interpolate_line(int startvertex,
         i<nnew;
         i++, NEXTV) {
       assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
-      fputc('#',stderr);
+      putchar('#');
       GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
                            &all.a[traverse.v][k]) );
       note_computed(traverse.v, k);
@@ -221,49 +231,73 @@ static void interpolate_line(int startvertex,
     gsl_interp_accel_free(accel);
     gsl_interp_free(interp);
   }
-  fprintf(stderr,"    done\n");
+  printf("    done\n");
 }
        
-static void interpolate(void) {
+static void interpolate_gsl(const gsl_interp_type *aperiodic,
+                           const gsl_interp_type *periodic) {
   int x,y;
   
-  for (y=0; y<(Y+1)/2; y+=INC) {
-    interpolate_line(y<<YSHIFT | ((y>>1)&1), 0, OX*2+1);
+  for (y=0; y<(Y+1)/2; y+=inc) {
+    interpolate_gsl_line(y<<YSHIFT | ((y>>1)&1), 0, oX*2+1,
+                        aperiodic, periodic);
   }
-  for (x=0; x<X; x+=INC) {
-    interpolate_line(                x, 5, OY);
-    interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
+  for (x=0; x<X; x+=inc) {
+    interpolate_gsl_line(                x, 5, oY, aperiodic, periodic);
+    interpolate_gsl_line((Y-1)<<YSHIFT | x, 1, oY, aperiodic, periodic);
   }
 }
-#endif
 
-static void write_output(void) {
+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++) {
-    fprintf(stderr,"checking y=%d ", y);
+    printf("checking y=%d ", y);
     for (x=0; x<X; x++) {
       v= y<<YSHIFT | x;
-      fprintf(stderr," %02x",v);
+      printf(" %02x",v);
       c= computed_count[v];
-      if (c==D3) fputc('#',stderr);
-      else if (c==-1) fputc('*',stderr);
-      else { fprintf(stderr,"!%d",c); bad++; }
+      if (c==D3) putchar('#');
+      else if (c==-1) putchar('*');
+      else { printf("!%d",c); bad++; }
     }
-    fputc('\n',stderr);
+    putchar('\n');
   }
   assert(!bad);
-  if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
-  fprintf(stderr,"interpolated\n");
+
+  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");
   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;
 }