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