chiark / gitweb /
f5d7c1bb9c727517e9bb3784f3cc5f64ac41d9b0
[moebius2.git] / 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 characterise_input(void) {
18   struct stat stab;
19   int r;
20   
21   r= fstat(0,&stab);  if (r) diee("fstat input to find length");
22
23   if (!stab.st_size || stab.st_size % (sizeof(double)*D3) ||
24       stab.st_size > INT_MAX)
25     fail("input file is not reasonable whole number of vertices\n");
26
27   fprintf(stderr,"XBITS=%d XY=%d*%d=%d\n", XBITS, X,Y,N);
28   fprintf(stderr,"OXBITS=%d OXY=%d*%d=%d\n", OXBITS, OX,OY,ON);
29   fprintf(stderr,"sizeof(double)=%d st_size=%lu\n",
30           (int)sizeof(double), (unsigned long)stab.st_size);
31
32   if (stab.st_size != (ON*sizeof(double)*D3))
33     fail("input file wrong size\n");
34 }
35
36 static void read_input(void) {
37   int x,y, ox,oy, v,ov;
38   
39   for (oy=y=0; y<Y; oy++, y+=INC) {
40     fprintf(stderr, "y=%2d>%2d", oy,y);
41     for (ox=x=0; x<X; ox++, x+=INC) {
42       errno= 0;
43       ov= (oy << OXBITS) | ox;
44       v= (y << YSHIFT) | x;
45       fprintf(stderr, " 0%02o->0x%02x", ov, v);
46       if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
47         diee("\nread input");
48       computed_count[v]= -1;
49     }
50     fputc('\n',stderr);
51   }
52 }
53
54   /* We use GSL's interpolation functions.  Each xa array is simple
55    * integers, and we do the interpolation three times for each line
56    * of topologically colinear vertices (once, separately, for each of
57    * the three spatial coordinates in the ya array).  The `horizontal'
58    * lines are done with periodic boundary conditions, obviously.
59    *
60    * This directly gives the required additional vertices:
61    *
62    *         ,O.                   O  original vertex
63    *        /   \                  *  new vertex
64    *      ,* - - *.                dashed lines are new edges,
65    *     /  `   '  \                not directly represented
66    *    O----`*'----O               or manipulated by this program
67    *
68    */
69
70 typedef struct {
71   int v, e;
72 } Traverse;
73
74 static void traverse_next(Traverse *t) {
75   int v2;
76   v2= EDGE_END2(t->v, t->e);
77   if (v2>=0) {
78     int e2= edge_reverse(t->v, t->e);
79     assert(EDGE_END2(v2,e2) == t->v);
80     t->e= (e2+3) % V6;
81   }
82   t->v= v2;
83 }
84      
85 static void interpolate_line(int startvertex,
86                              int direction /* edge number */,
87                              int nmax) {
88   double xa[nmax], ya[D3][nmax];
89   const gsl_interp_type *it;
90   gsl_interp *interp;
91   gsl_interp_accel *accel;
92   int i, k, nold, nnew;
93   Traverse traverse;
94
95 #define STARTV (traverse.v=startvertex, traverse.e=direction)
96 #define NEXTV (traverse_next(&traverse), fprintf(stderr,"-%02x",traverse.v))
97
98   fprintf(stderr,"interpolate_line 0x%02x,%d nmax=%d ",
99           startvertex,direction,nmax);
100
101   for (i=0, STARTV;
102        ;
103        ) {
104     assert(traverse.v>=0);
105     assert(i < nmax);
106     xa[i]= i;
107     K ya[k][i]= all.a[traverse.v][k];
108     fputc('*',stderr);
109     if (i && traverse.v==startvertex) { fputc('@',stderr); break; }
110     i++;
111     NEXTV;
112     if (traverse.v<0) break;
113     NEXTV;
114   }
115   if (traverse.v>=0) {
116     it= gsl_interp_akima_periodic;
117     nold= i+1;
118     nnew= i;
119   } else {
120     it= gsl_interp_akima;
121     nold= i;
122     nnew= i-1;
123   }
124
125   fprintf(stderr,"  n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
126
127   K {
128     fprintf(stderr,"\n  k=%d ", k);
129
130     interp= gsl_interp_alloc(it,nold); if (!interp) diee("gsl_interp_alloc");
131     accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
132     GA( gsl_interp_init(interp,xa,ya[k],nold) );
133
134     for (i=0, STARTV;
135          i<nnew;
136          i++, NEXTV) {
137       assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
138       fputc('#',stderr);
139       assert(computed_count[traverse.v] == k);
140       GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
141                             &all.a[traverse.v][k]) );
142       computed_count[traverse.v]++;
143     }
144     
145     gsl_interp_accel_free(accel);
146     gsl_interp_free(interp);
147   }
148   fprintf(stderr,"    done\n");
149 }
150         
151 static void interpolate(void) {
152   int x,y;
153   
154   for (y=0; y<(Y+1)/2; y+=INC) {
155     interpolate_line(y<<YSHIFT, 0, OX*2+1);
156   }
157   for (x=0; x<X; x+=INC) {
158     interpolate_line(                x, 5, OY);
159     interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
160   }
161 }
162
163 static void write_output(void) {
164   int x, y, v, c, bad=0;
165   for (y=0; y<Y; y++) {
166     fprintf(stderr,"checking y=%d ", y);
167     for (x=0; x<X; x++) {
168       v= y<<YSHIFT | x;
169       fprintf(stderr," %02x",v);
170       c= computed_count[v];
171       if (c==D3) fputc('#',stderr);
172       else if (c==-1) fputc('*',stderr);
173       else { fputc('!',stderr); bad++; }
174     }
175     fputc('\n',stderr);
176   }
177   assert(!bad);
178   if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
179   fprintf(stderr,"interpolated\n");
180 }
181
182 int main(int argc, const char **argv) {
183   if (argc!=1 || isatty(0) || isatty(1))
184     { fputs("usage: interpolate <INPUT.cfm >OUTPUT.cfm\n",stderr); exit(8); }
185
186   characterise_input();
187   read_input();
188   interpolate();
189   write_output();
190   if (fclose(stdout)) diee("fclose stdout");
191   return 0;
192 }