2 * increases the resolution of a conformation by interpolation
7 /* filled in by characterise_input: */
8 static int shift, oldxbits, oldybits, oldx, oldy, oldsz, inc;
10 /* filled in by read_input: */
11 static struct Vertices all;
13 static void characterise_input(void) {
17 r= fstat(0,&stab); if (r) diee("fstat input to find length");
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");
23 oldsz= stab.st_size / sizeof(double);
25 shift > XBITS+1 && shift > YBITS+1;
30 oldy= (1<<oldybits)-1;
31 if (oldx*oldy == oldsz) goto found;
33 fail("input file size cannot be interpolated to target file size\n");
39 static void read_input(void) {
42 for (y=0; y<Y; y+=inc)
43 for (x=0; x<X; x+=inc) {
45 if (fread(all.a[(y << YSHIFT) | x], sizeof(double), D3, stdin) != D3)
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.
56 * This directly gives the required additional vertices:
58 * ,O. O original vertex
60 * ,* - - *. dashed lines are new edges,
61 * / ` ' \ not directly represented
62 * O----`*'----O or manipulated by this program
66 #define NEXTV (v= EDGE_END2(v,direction))
68 static void interpolate_line(int startvertex,
69 int direction /* edge number */,
71 double xa[norig+1], ya[D3][norig+1], n;
72 const gsl_interp_type *it;
74 gsl_interp_accel *accel;
77 for (i=0, v=startvertex;
79 i++, NEXTV, assert(v>=0), NEXTV) {
84 K ya[k][i]= all.a[v][k];
89 it= v>=0 ? gsl_interp_cspline_periodic : gsl_interp_cspline;
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) );
96 for (i=0, v=startvertex;
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]) );
103 gsl_interp_accel_free(accel);
104 gsl_interp_free(interp);
108 static void interpolate(void) {
111 if (shift!=1) fail("only interpolation by factor of 2 supported\n");
113 for (y=0; y<Y; y+=inc) {
114 interpolate_line(y<<YSHIFT, 0, oldx);
116 for (x=0; x<X; x+=inc) {
117 interpolate_line(x, 5, oldy);
118 interpolate_line(x, 4, oldy);
122 static void write_output(void) {
123 if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
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); }
130 characterise_input();
134 if (fclose(stdout)) diee("fclose stdout");