interpolate.c

   1 /*
   2  * increases the resolution of a conformation by interpolation
   3  */
   4
   5 #include "mgraph.h"
   6
   7 /* filled in by characterise_input: */
   8 static int shift, oldxbits, oldybits, oldx, oldy, oldsz, inc;
   9
  10 /* filled in by read_input: */
  11 static struct Vertices all;
  12
  13 static void characterise_input(void) {
  14   struct stat stab;
  15   int r;
  16
  17   r= fstat(0,&stab);  if (r) diee("fstat input to find length");
  18
  19   if (!stab.st_size || stab.st_size % (sizeof(double)*D3) ||
  20       stab.st_size > INT_MAX)
  21     fail("input file is not reasonable whole number of doubles\n");
  22
  23   oldsz= stab.st_size / (sizeof(double)*D3);
  24   for (shift=1;
  25        shift < XBITS+1 && shift < YBITS+1;
  26        shift++) {
  27     oldxbits= XBITS-1;
  28     oldybits= YBITS-1;
  29     oldx=  1<<oldxbits;
  30     oldy= (1<<oldybits)-1;
  31     fprintf(stderr,"sizeof(double)=%d XYBITS=%d,%d, XY=%d*%d=%d"
  32             " oldsz=%d shift=%d oldxybits=%d,%d oldxy=%d*%d=%d\n",
  33             (int)sizeof(double), XBITS,YBITS, X,Y,N,
  34             oldsz,shift,oldxbits,oldybits,oldx,oldy,oldx*oldy);
  35     if (oldx*oldy == oldsz) goto found;
  36   }
  37   fail("input file size cannot be interpolated to target file size\n");
  38
  39  found:
  40   inc= 1<<shift;
  41   fprintf(stderr,"inc=%d\n",inc);
  42 }
  43
  44 static void read_input(void) {
  45   int x,y, ox,oy, v,ov;
  46
  47   for (oy=y=0; y<Y; oy++, y+=inc) {
  48     fprintf(stderr, "y=%2d>%2d", oy,y);
  49     for (ox=x=0; x<X; ox++, x+=inc) {
  50       errno= 0;
  51       ov= (oy << oldxbits) | ox;
  52       v= (y << YSHIFT) | x;
  53       fprintf(stderr, " 0%02o->0x%02x", ov, v);
  54       if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
  55         diee("\nread input");
  56     }
  57     fputc('\n',stderr);
  58   }
  59 }
  60
  61   /* We use GSL's interpolation functions.  Each xa array is simple
  62    * integers, and we do the interpolation three times for each line
  63    * of topologically colinear vertices (once, separately, for each of
  64    * the three spatial coordinates in the ya array).  The `horizontal'
  65    * lines are done with periodic boundary conditions, obviously.
  66    *
  67    * This directly gives the required additional vertices:
  68    *
  69    *         ,O.                   O  original vertex
  70    *        /   \                  *  new vertex
  71    *      ,* - - *.                dashed lines are new edges,
  72    *     /  `   '  \                not directly represented
  73    *    O----`*'----O               or manipulated by this program
  74    *
  75    */
  76
  77 #define NEXTV (v= EDGE_END2(v,direction))
  78
  79 static void interpolate_line(int startvertex,
  80                              int direction /* edge number */,
  81                              int norig) {
  82   double xa[norig+1], ya[D3][norig+1], n;
  83   const gsl_interp_type *it;
  84   gsl_interp *interp;
  85   gsl_interp_accel *accel;
  86   int i, v, k;
  87
  88   for (i=0, v=startvertex;
  89        ;
  90        i++, NEXTV, assert(v>=0), NEXTV) {
  91     if (v<0)
  92       break;
  93     assert(i <= norig);
  94     xa[i]= i;
  95     K ya[k][i]= all.a[v][k];
  96     if (v==startvertex)
  97       break;
  98   }
  99   n= i;
 100   it= v>=0 ? gsl_interp_cspline_periodic : gsl_interp_cspline;
 101
 102   K {
 103     interp= gsl_interp_alloc(it,n); if (!interp) diee("gsl_interp_alloc");
 104     accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
 105     GA( gsl_interp_init(interp,xa,ya[k],n) );
 106
 107     for (i=0, v=startvertex;
 108          i<n-1;
 109          i++, NEXTV) {
 110       assert(v>=0); NEXTV; assert(v>=0);
 111       GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel, &all.a[v][k]) );
 112     }
 113
 114     gsl_interp_accel_free(accel);
 115     gsl_interp_free(interp);
 116   }
 117 }
 118
 119 static void interpolate(void) {
 120   int x,y;
 121
 122   if (shift!=1) fail("only interpolation by factor of 2 supported\n");
 123
 124   for (y=0; y<Y; y+=inc) {
 125     interpolate_line(y<<YSHIFT, 0, oldx);
 126   }
 127   for (x=0; x<X; x+=inc) {
 128     interpolate_line(x, 5, oldy);
 129     interpolate_line(x, 4, oldy);
 130   }
 131 }
 132
 133 static void write_output(void) {
 134   if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
 135 }
 136
 137 int main(int argc, const char **argv) {
 138   if (argc!=1 || isatty(0) || isatty(1))
 139     { fputs("usage: interpolate <INPUT.cfm >OUTPUT.cfm\n",stderr); exit(8); }
 140
 141   characterise_input();
 142   read_input();
 143   interpolate();
 144   write_output();
 145   if (fclose(stdout)) diee("fclose stdout");
 146   return 0;
 147 }