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) ||
  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);
  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     if (oldx*oldy == oldsz) goto found;
  32   }
  33   fail("input file size cannot be interpolated to target file size\n");
  34
  35  found:
  36   inc= 1<<shift;
  37 }
  38
  39 static void read_input(void) {
  40   int x,y;
  41
  42   for (y=0; y<Y; y+=inc)
  43     for (x=0; x<X; x+=inc) {
  44       errno= 0;
  45       if (fread(all.a[(y << YSHIFT) | x], sizeof(double), D3, stdin) != D3)
  46         diee("read input");
  47     }
  48 }
  49
  50   /* We use GSL's interpolation functions.  Each xa array is simple
  51    * integers, and we do the interpolation three times for each line
  52    * of topologically colinear vertices (once, separately, for each of
  53    * the three spatial coordinates in the ya array).  The `horizontal'
  54    * lines are done with periodic boundary conditions, obviously.
  55    *
  56    * This directly gives the required additional vertices:
  57    *
  58    *         ,O.                   O  original vertex
  59    *        /   \                  *  new vertex
  60    *      ,* - - *.                dashed lines are new edges,
  61    *     /  `   '  \                not directly represented
  62    *    O----`*'----O               or manipulated by this program
  63    *
  64    */
  65
  66 #define NEXTV (v= EDGE_END2(v,direction))
  67
  68 static void interpolate_line(int startvertex,
  69                              int direction /* edge number */,
  70                              int norig) {
  71   double xa[norig+1], ya[D3][norig+1], n;
  72   const gsl_interp_type *it;
  73   gsl_interp *interp;
  74   gsl_interp_accel *accel;
  75   int i, v, k;
  76
  77   for (i=0, v=startvertex;
  78        ;
  79        i++, NEXTV, assert(v>=0), NEXTV) {
  80     if (v<0)
  81       break;
  82     assert(i <= norig);
  83     xa[i]= i;
  84     K ya[k][i]= all.a[v][k];
  85     if (v==startvertex)
  86       break;
  87   }
  88   n= i;
  89   it= v>=0 ? gsl_interp_cspline_periodic : gsl_interp_cspline;
  90
  91   K {
  92     interp= gsl_interp_alloc(it,n); if (!interp) diee("gsl_interp_alloc");
  93     accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
  94     GA( gsl_interp_init(interp,xa,ya[k],n) );
  95
  96     for (i=0, v=startvertex;
  97          i<n-1;
  98          i++, NEXTV) {
  99       assert(v>=0); NEXTV; assert(v>=0);
 100       GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel, &all.a[v][k]) );
 101     }
 102
 103     gsl_interp_accel_free(accel);
 104     gsl_interp_free(interp);
 105   }
 106 }
 107
 108 static void interpolate(void) {
 109   int x,y;
 110
 111   if (shift!=1) fail("only interpolation by factor of 2 supported\n");
 112
 113   for (y=0; y<Y; y+=inc) {
 114     interpolate_line(y<<YSHIFT, 0, oldx);
 115   }
 116   for (x=0; x<X; x+=inc) {
 117     interpolate_line(x, 5, oldy);
 118     interpolate_line(x, 4, oldy);
 119   }
 120 }
 121
 122 static void write_output(void) {
 123   if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
 124 }
 125
 126 int main(int argc, const char **argv) {
 127   if (argc!=1 || isatty(0) || isatty(1))
 128     { fputs("usage: interpolate <INPUT.cfm >OUTPUT.cfm\n",stderr); exit(8); }
 129
 130   characterise_input();
 131   read_input();
 132   interpolate();
 133   write_output();
 134   if (fclose(stdout)) diee("fclose stdout");
 135   return 0;
 136 }