interpolate.c

   1 /*
   2  * increases the resolution of a conformation by interpolation
   3  */
   4
   5 #include "mgraph.h"
   6
   7 #define OXBITS (XBITS-1)
   8 #define OX (1<<OXBITS)
   9 #define OY (Y/2 + 1)
  10 #define ON (OX*OY)
  11 #define INC 2
  12
  13 /* filled in by read_input: */
  14 static struct Vertices all;
  15 static int computed_count[N];
  16
  17 static void characterise_input(void) {
  18   struct stat stab;
  19   int r;
  20
  21   r= fstat(0,&stab);  if (r) diee("fstat input to find length");
  22
  23   if (!stab.st_size || stab.st_size % (sizeof(double)*D3) ||
  24       stab.st_size > INT_MAX)
  25     fail("input file is not reasonable whole number of vertices\n");
  26
  27   fprintf(stderr,"XBITS=%d XY=%d*%d=%d\n", XBITS, X,Y,N);
  28   fprintf(stderr,"OXBITS=%d OXY=%d*%d=%d\n", OXBITS, OX,OY,ON);
  29   fprintf(stderr,"sizeof(double)=%d st_size=%lu\n",
  30           (int)sizeof(double), (unsigned long)stab.st_size);
  31
  32   if (stab.st_size != (ON*sizeof(double)*D3))
  33     fail("input file wrong size\n");
  34 }
  35
  36 static void read_input(void) {
  37   int x,y, ox,oy, v,ov;
  38
  39   for (oy=y=0; y<Y; oy++, y+=INC) {
  40     fprintf(stderr, "y=%2d>%2d", oy,y);
  41     for (ox=x=0; x<X; ox++, x+=INC) {
  42       errno= 0;
  43       ov= (oy << OXBITS) | ox;
  44       v= (y << YSHIFT) | x;
  45       fprintf(stderr, " 0%02o->0x%02x", ov, v);
  46       if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
  47         diee("\nread input");
  48       computed_count[v]= -1;
  49     }
  50     fputc('\n',stderr);
  51   }
  52 }
  53
  54   /* We use GSL's interpolation functions.  Each xa array is simple
  55    * integers, and we do the interpolation three times for each line
  56    * of topologically colinear vertices (once, separately, for each of
  57    * the three spatial coordinates in the ya array).  The `horizontal'
  58    * lines are done with periodic boundary conditions, obviously.
  59    *
  60    * This directly gives the required additional vertices:
  61    *
  62    *         ,O.                   O  original vertex
  63    *        /   \                  *  new vertex
  64    *      ,* - - *.                dashed lines are new edges,
  65    *     /  `   '  \                not directly represented
  66    *    O----`*'----O               or manipulated by this program
  67    *
  68    */
  69
  70 typedef struct {
  71   int v, e;
  72 } Traverse;
  73
  74 static void traverse_next(Traverse *t) {
  75   int v2;
  76   v2= EDGE_END2(t->v, t->e);
  77   if (v2>=0) {
  78     int e2= edge_reverse(t->v, t->e);
  79     assert(EDGE_END2(v2,e2) == t->v);
  80     t->e= (e2+3) % V6;
  81   }
  82   t->v= v2;
  83 }
  84
  85 static void interpolate_line(int startvertex,
  86                              int direction /* edge number */,
  87                              int nmax) {
  88   double xa[nmax], ya[D3][nmax];
  89   const gsl_interp_type *it;
  90   gsl_interp *interp;
  91   gsl_interp_accel *accel;
  92   int i, k, nold, nnew;
  93   Traverse traverse;
  94
  95 #define STARTV (traverse.v=startvertex, traverse.e=direction)
  96 #define NEXTV (traverse_next(&traverse), fprintf(stderr,"-%02x",traverse.v))
  97
  98   fprintf(stderr,"interpolate_line 0x%02x,%d nmax=%d ",
  99           startvertex,direction,nmax);
 100
 101   for (i=0, STARTV;
 102        ;
 103        ) {
 104     assert(traverse.v>=0);
 105     assert(i < nmax);
 106     xa[i]= i;
 107     K ya[k][i]= all.a[traverse.v][k];
 108     fputc('*',stderr);
 109     if (i && traverse.v==startvertex) { fputc('@',stderr); break; }
 110     i++;
 111     NEXTV;
 112     if (traverse.v<0) break;
 113     NEXTV;
 114   }
 115   if (traverse.v>=0) {
 116     it= gsl_interp_akima_periodic;
 117     nold= i+1;
 118     nnew= i;
 119   } else {
 120     it= gsl_interp_akima;
 121     nold= i;
 122     nnew= i-1;
 123   }
 124
 125   fprintf(stderr,"  n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
 126
 127   K {
 128     fprintf(stderr,"\n  k=%d ", k);
 129
 130     interp= gsl_interp_alloc(it,nold); if (!interp) diee("gsl_interp_alloc");
 131     accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
 132     GA( gsl_interp_init(interp,xa,ya[k],nold) );
 133
 134     for (i=0, STARTV;
 135          i<nnew;
 136          i++, NEXTV) {
 137       assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
 138       fputc('#',stderr);
 139       assert(computed_count[traverse.v] == k);
 140       GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
 141                             &all.a[traverse.v][k]) );
 142       computed_count[traverse.v]++;
 143     }
 144
 145     gsl_interp_accel_free(accel);
 146     gsl_interp_free(interp);
 147   }
 148   fprintf(stderr,"    done\n");
 149 }
 150
 151 static void interpolate(void) {
 152   int x,y;
 153
 154   for (y=0; y<(Y+1)/2; y+=INC) {
 155     interpolate_line(y<<YSHIFT, 0, OX*2+1);
 156   }
 157   for (x=0; x<X; x+=INC) {
 158     interpolate_line(                x, 5, OY);
 159     interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
 160   }
 161 }
 162
 163 static void write_output(void) {
 164   int x, y, v, c, bad=0;
 165   for (y=0; y<Y; y++) {
 166     fprintf(stderr,"checking y=%d ", y);
 167     for (x=0; x<X; x++) {
 168       v= y<<YSHIFT | x;
 169       fprintf(stderr," %02x",v);
 170       c= computed_count[v];
 171       if (c==D3) fputc('#',stderr);
 172       else if (c==-1) fputc('*',stderr);
 173       else { fputc('!',stderr); bad++; }
 174     }
 175     fputc('\n',stderr);
 176   }
 177   assert(!bad);
 178   if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
 179   fprintf(stderr,"interpolated\n");
 180 }
 181
 182 int main(int argc, const char **argv) {
 183   if (argc!=1 || isatty(0) || isatty(1))
 184     { fputs("usage: interpolate <INPUT.cfm >OUTPUT.cfm\n",stderr); exit(8); }
 185
 186   characterise_input();
 187   read_input();
 188   interpolate();
 189   write_output();
 190   if (fclose(stdout)) diee("fclose stdout");
 191   return 0;
 192 }