interpolate and various sizes build

[moebius2.git] / interpolate.c

diff --git a/interpolate.c b/interpolate.c
new file mode 100644 (file)
index 0000000..6f519e4
--- /dev/null
+++ b/interpolate.c
@@ -0,0 +1,136 @@
+/*
+ * increases the resolution of a conformation by interpolation
+ */
+
+#include "mgraph.h"
+
+/* filled in by characterise_input: */
+static int shift, oldxbits, oldybits, oldx, oldy, oldsz, inc;
+
+/* filled in by read_input: */
+static struct Vertices all;
+
+static void characterise_input(void) {
+  struct stat stab;
+  int r;
+  
+  r= fstat(0,&stab);  if (r) diee("fstat input to find length");
+
+  if (!stab.st_size || stab.st_size % sizeof(double) ||
+      stab.st_size > INT_MAX)
+    fail("input file is not reasonable whole number of doubles\n");
+
+  oldsz= stab.st_size / sizeof(double);
+  for (shift=1;
+       shift > XBITS+1 && shift > YBITS+1;
+       shift++) {
+    oldxbits= XBITS-1;
+    oldybits= YBITS-1;
+    oldx=  1<<oldxbits;
+    oldy= (1<<oldybits)-1;
+    if (oldx*oldy == oldsz) goto found;
+  }
+  fail("input file size cannot be interpolated to target file size\n");
+
+ found:
+  inc= 1<<shift;
+}
+
+static void read_input(void) {
+  int x,y;
+  
+  for (y=0; y<Y; y+=inc)
+    for (x=0; x<X; x+=inc) {
+      errno= 0;
+      if (fread(all.a[(y << YSHIFT) | x], sizeof(double), D3, stdin) != D3)
+       diee("read input");
+    }
+}
+
+  /* We use GSL's interpolation functions.  Each xa array is simple
+   * integers, and we do the interpolation three times for each line
+   * of topologically colinear vertices (once, separately, for each of
+   * the three spatial coordinates in the ya array).  The `horizontal'
+   * lines are done with periodic boundary conditions, obviously.
+   *
+   * This directly gives the required additional vertices:
+   *
+   *         ,O.                   O  original vertex
+   *        /   \                  *  new vertex
+   *      ,* - - *.                dashed lines are new edges,
+   *     /  `   '  \                not directly represented
+   *    O----`*'----O               or manipulated by this program
+   *
+   */
+
+#define NEXTV (v= EDGE_END2(v,direction))
+
+static void interpolate_line(int startvertex,
+                            int direction /* edge number */,
+                            int norig) {
+  double xa[norig+1], ya[D3][norig+1], n;
+  const gsl_interp_type *it;
+  gsl_interp *interp;
+  gsl_interp_accel *accel;
+  int i, v, k;
+
+  for (i=0, v=startvertex;
+       ;
+       i++, NEXTV, assert(v>=0), NEXTV) {
+    if (v<0)
+      break;
+    assert(i <= norig);
+    xa[i]= i;
+    K ya[k][i]= all.a[v][k];
+    if (v==startvertex)
+      break;
+  }
+  n= i;
+  it= v>=0 ? gsl_interp_cspline_periodic : gsl_interp_cspline;
+
+  K {
+    interp= gsl_interp_alloc(it,n); if (!interp) diee("gsl_interp_alloc");
+    accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
+    GA( gsl_interp_init(interp,xa,ya[k],n) );
+
+    for (i=0, v=startvertex;
+        i<n-1;
+        i++, NEXTV) {
+      assert(v>=0); NEXTV; assert(v>=0);
+      GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel, &all.a[v][k]) );
+    }
+    
+    gsl_interp_accel_free(accel);
+    gsl_interp_free(interp);
+  }
+}
+       
+static void interpolate(void) {
+  int x,y;
+  
+  if (shift!=1) fail("only interpolation by factor of 2 supported\n");
+
+  for (y=0; y<Y; y+=inc) {
+    interpolate_line(y<<YSHIFT, 0, oldx);
+  }
+  for (x=0; x<X; x+=inc) {
+    interpolate_line(x, 5, oldy);
+    interpolate_line(x, 4, oldy);
+  }
+}
+
+static void write_output(void) {
+  if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
+}
+
+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); }
+
+  characterise_input();
+  read_input();
+  interpolate();
+  write_output();
+  if (fclose(stdout)) diee("fclose stdout");
+  return 0;
+}