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 note_computed(int v, int k) {
  18   assert(computed_count[v] == k);
  19   computed_count[v]++;
  20 }
  21
  22 static void characterise_input(void) {
  23   struct stat stab;
  24   int r;
  25
  26   r= fstat(0,&stab);  if (r) diee("fstat input to find length");
  27
  28   if (!stab.st_size || stab.st_size % (sizeof(double)*D3) ||
  29       stab.st_size > INT_MAX)
  30     fail("input file is not reasonable whole number of vertices\n");
  31
  32   fprintf(stderr,"XBITS=%d XY=%d*%d=%d\n", XBITS, X,Y,N);
  33   fprintf(stderr,"OXBITS=%d OXY=%d*%d=%d\n", OXBITS, OX,OY,ON);
  34   fprintf(stderr,"sizeof(double)=%d st_size=%lu\n",
  35           (int)sizeof(double), (unsigned long)stab.st_size);
  36
  37   if (stab.st_size != (ON*sizeof(double)*D3))
  38     fail("input file wrong size\n");
  39 }
  40
  41 static void read_input(void) {
  42   int x,y, ox,oy, v,ov;
  43
  44   for (oy=0; oy<OY; oy++) {
  45     y= oy*2;
  46     fprintf(stderr, "y=%2d>%2d", oy,y);
  47     for (ox=0; ox<OX; ox++) {
  48       x= ox*2 + (oy&1);
  49       if (x>=X) { y= (Y-1)-y; x-=X; }
  50       errno= 0;
  51       ov= (oy << OXBITS) | 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       computed_count[v]= -1;
  57     }
  58     fputc('\n',stderr);
  59   }
  60 }
  61
  62   /* We use GSL's interpolation functions.  Each xa array is simple
  63    * integers, and we do the interpolation three times for each line
  64    * of topologically colinear vertices (once, separately, for each of
  65    * the three spatial coordinates in the ya array).  The `horizontal'
  66    * lines are done with periodic boundary conditions, obviously.
  67    *
  68    * This directly gives the required additional vertices:
  69    *
  70    *         ,O.                   O  original vertex
  71    *        /   \                  *  new vertex
  72    *      ,* - - *.                dashed lines are new edges,
  73    *     /  `   '  \                not directly represented
  74    *    O----`*'----O               or manipulated by this program
  75    *
  76    */
  77
  78 typedef struct {
  79   int v, e;
  80 } Traverse;
  81
  82 static void traverse_next(Traverse *t) {
  83   int v2;
  84
  85   if (t->v<0)
  86     return;
  87
  88   v2= EDGE_END2(t->v, t->e);
  89   if (v2>=0) {
  90     int e2= edge_reverse(t->v, t->e);
  91     assert(EDGE_END2(v2,e2) == t->v);
  92     t->e= (e2+3) % V6;
  93   }
  94   t->v= v2;
  95 }
  96
  97 #if 0
  98
  99 static void interpolate(void) {
 100   /* four points P Q R S, although P and S may be missing
 101    * interpolate in QR finding M. */
 102   int xq,yq,eqr, vp,vq,vm,vr,vs, k;
 103   double pqtarg[D3], srtarg[D3];
 104
 105   for (eqr=1; eqr!=4; eqr+=5, eqr%=V6) { /* each old edge exactly once */
 106     fprintf(stderr,"eqr=%d\n",eqr);
 107     for (yq=0; yq<Y; yq+=INC) {
 108       fprintf(stderr," yq=%2d ",yq);
 109       for (xq=((yq>>1)&1); xq<X; xq+=INC) {
 110         vq= yq << YSHIFT | xq;
 111         Traverse trav; trav.v=vq; trav.e=eqr;
 112
 113         traverse_next(&trav);
 114         vm= trav.v;
 115         if (vm<0) continue;
 116
 117         traverse_next(&trav);
 118         vr= trav.v;
 119         assert(vr>=0);
 120
 121         traverse_next(&trav);
 122         traverse_next(&trav);
 123         vs= trav.v;
 124
 125         trav.v= vq; trav.e= EDGE_OPPOSITE(eqr);
 126         traverse_next(&trav);
 127         traverse_next(&trav);
 128         vp= trav.v;
 129
 130         fprintf(stderr," 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
 131                 vp&0xff,vq,vm,vr,vs&0xff);
 132
 133         if (vp>=0)
 134           K pqtarg[k]= all.a[vq][k]*2 - all.a[vp][k];
 135         else
 136           K pqtarg[k]= all.a[vr][k];
 137
 138         if (vs>=0)
 139           K srtarg[k]= all.a[vr][k]*2 - all.a[vs][k];
 140         else
 141           K srtarg[k]= all.a[vq][k];
 142
 143         K {
 144           const double alpha= 3.0;
 145           all.a[vm][k]= 0.5 * (all.a[vq][k] + all.a[vr][k]);
 146           pqtarg[k]= 0.5 * (pqtarg[k] + all.a[vm][k]);
 147           srtarg[k]= 0.5 * (srtarg[k] + all.a[vm][k]);
 148           all.a[vm][k]= (pqtarg[k] + srtarg[k] + alpha * all.a[vm][k])
 149             / (2 + alpha);
 150           note_computed(vm,k);
 151         }
 152       }
 153       fputc('\n',stderr);
 154     }
 155   }
 156 }
 157
 158 #endif
 159
 160 #if 1
 161
 162 static void interpolate_line(int startvertex,
 163                              int direction /* edge number */,
 164                              int nmax) {
 165   double xa[nmax], ya[D3][nmax];
 166   const gsl_interp_type *it;
 167   gsl_interp *interp;
 168   gsl_interp_accel *accel;
 169   int i, k, nold, nnew;
 170   Traverse traverse;
 171
 172 #define STARTV (traverse.v=startvertex, traverse.e=direction)
 173 #define NEXTV (traverse_next(&traverse), fprintf(stderr,"-%02x",traverse.v))
 174
 175   fprintf(stderr,"interpolate_line 0x%02x,%d nmax=%d ",
 176           startvertex,direction,nmax);
 177
 178   for (i=0, STARTV;
 179        ;
 180        ) {
 181     assert(traverse.v>=0);
 182     assert(i < nmax);
 183     xa[i]= i;
 184     K ya[k][i]= all.a[traverse.v][k];
 185     fputc('*',stderr);
 186     if (i && traverse.v==startvertex) { fputc('@',stderr); break; }
 187     i++;
 188     NEXTV;
 189     if (traverse.v<0) break;
 190     NEXTV;
 191   }
 192   if (traverse.v>=0) {
 193     it= gsl_interp_akima_periodic;
 194     nold= i+1;
 195     nnew= i;
 196   } else {
 197     it= gsl_interp_akima;
 198     nold= i;
 199     nnew= i-1;
 200   }
 201
 202   fprintf(stderr,"  n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
 203
 204   K {
 205     fprintf(stderr,"\n  k=%d ", k);
 206
 207     interp= gsl_interp_alloc(it,nold); if (!interp) diee("gsl_interp_alloc");
 208     accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
 209     GA( gsl_interp_init(interp,xa,ya[k],nold) );
 210
 211     for (i=0, STARTV;
 212          i<nnew;
 213          i++, NEXTV) {
 214       assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
 215       fputc('#',stderr);
 216       GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
 217                             &all.a[traverse.v][k]) );
 218       note_computed(traverse.v, k);
 219     }
 220
 221     gsl_interp_accel_free(accel);
 222     gsl_interp_free(interp);
 223   }
 224   fprintf(stderr,"    done\n");
 225 }
 226
 227 static void interpolate(void) {
 228   int x,y;
 229
 230   for (y=0; y<(Y+1)/2; y+=INC) {
 231     interpolate_line(y<<YSHIFT | ((y>>1)&1), 0, OX*2+1);
 232   }
 233   for (x=0; x<X; x+=INC) {
 234     interpolate_line(                x, 5, OY);
 235     interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
 236   }
 237 }
 238 #endif
 239
 240 static void write_output(void) {
 241   int x, y, v, c, bad=0;
 242   for (y=0; y<Y; y++) {
 243     fprintf(stderr,"checking y=%d ", y);
 244     for (x=0; x<X; x++) {
 245       v= y<<YSHIFT | x;
 246       fprintf(stderr," %02x",v);
 247       c= computed_count[v];
 248       if (c==D3) fputc('#',stderr);
 249       else if (c==-1) fputc('*',stderr);
 250       else { fprintf(stderr,"!%d",c); bad++; }
 251     }
 252     fputc('\n',stderr);
 253   }
 254   assert(!bad);
 255   if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
 256   fprintf(stderr,"interpolated\n");
 257 }
 258
 259 int main(int argc, const char **argv) {
 260   if (argc!=1 || isatty(0) || isatty(1))
 261     { fputs("usage: interpolate <INPUT.cfm >OUTPUT.cfm\n",stderr); exit(8); }
 262
 263   characterise_input();
 264   read_input();
 265   interpolate();
 266   write_output();
 267   if (fclose(stdout)) diee("fclose stdout");
 268   return 0;
 269 }