2 * increases the resolution of a conformation by interpolation
7 #define OXBITS (XBITS-1)
13 /* filled in by read_input: */
14 static struct Vertices all;
15 static int computed_count[N];
17 static void note_computed(int v, int k) {
18 assert(computed_count[v] == k);
22 static void characterise_input(void) {
26 r= fstat(0,&stab); if (r) diee("fstat input to find length");
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");
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);
37 if (stab.st_size != (ON*sizeof(double)*D3))
38 fail("input file wrong size\n");
41 static void read_input(void) {
44 for (oy=0; oy<OY; oy++) {
46 fprintf(stderr, "y=%2d>%2d", oy,y);
47 for (ox=0; ox<OX; ox++) {
49 if (x>=X) { y= (Y-1)-y; x-=X; }
51 ov= (oy << OXBITS) | ox;
53 fprintf(stderr, " 0%02o->0x%02x", ov, v);
54 if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
56 computed_count[v]= -1;
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.
68 * This directly gives the required additional vertices:
70 * ,O. O original vertex
72 * ,* - - *. dashed lines are new edges,
73 * / ` ' \ not directly represented
74 * O----`*'----O or manipulated by this program
82 static void traverse_next(Traverse *t) {
88 v2= EDGE_END2(t->v, t->e);
90 int e2= edge_reverse(t->v, t->e);
91 assert(EDGE_END2(v2,e2) == t->v);
97 static void interpolate(void) {
98 /* four points P Q R S, although P and S may be missing
99 * interpolate in QR finding M. */
100 int xq,yq,eqr, vp,vq,vm,vr,vs, k;
101 double pqtarg[D3], srtarg[D3];
103 for (eqr=1; eqr!=4; eqr+=5, eqr%=V6) { /* each old edge exactly once */
104 fprintf(stderr,"eqr=%d\n",eqr);
105 for (yq=0; yq<Y; yq+=INC) {
106 fprintf(stderr," yq=%2d ",yq);
107 for (xq=((yq>>1)&1); xq<X; xq+=INC) {
108 vq= yq << YSHIFT | xq;
109 Traverse trav; trav.v=vq; trav.e=eqr;
111 traverse_next(&trav);
115 traverse_next(&trav);
119 traverse_next(&trav);
120 traverse_next(&trav);
123 trav.v= vq; trav.e= EDGE_OPPOSITE(eqr);
124 traverse_next(&trav);
125 traverse_next(&trav);
128 fprintf(stderr," 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
129 vp&0xff,vq,vm,vr,vs&0xff);
132 K pqtarg[k]= all.a[vq][k]*2 - all.a[vp][k];
134 K pqtarg[k]= all.a[vr][k];
137 K srtarg[k]= all.a[vr][k]*2 - all.a[vs][k];
139 K srtarg[k]= all.a[vq][k];
142 all.a[vm][k]= 0.5 * (all.a[vq][k] + all.a[vr][k]);
143 // pqtarg[k]= 0.5 * (pqtarg[k] + all.a[vm][k]);
144 // srtarg[k]= 0.5 * (srtarg[k] + all.a[vm][k]);
145 // all.a[vm][k]= 0.5 * (pqtarg[k] + srtarg[k]);
156 static void interpolate_line(int startvertex,
157 int direction /* edge number */,
159 double xa[nmax], ya[D3][nmax];
160 const gsl_interp_type *it;
162 gsl_interp_accel *accel;
163 int i, k, nold, nnew;
166 #define STARTV (traverse.v=startvertex, traverse.e=direction)
167 #define NEXTV (traverse_next(&traverse), fprintf(stderr,"-%02x",traverse.v))
169 fprintf(stderr,"interpolate_line 0x%02x,%d nmax=%d ",
170 startvertex,direction,nmax);
175 assert(traverse.v>=0);
178 K ya[k][i]= all.a[traverse.v][k];
180 if (i && traverse.v==startvertex) { fputc('@',stderr); break; }
183 if (traverse.v<0) break;
187 it= gsl_interp_akima_periodic;
191 it= gsl_interp_akima;
196 fprintf(stderr," n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
199 fprintf(stderr,"\n k=%d ", k);
201 interp= gsl_interp_alloc(it,nold); if (!interp) diee("gsl_interp_alloc");
202 accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
203 GA( gsl_interp_init(interp,xa,ya[k],nold) );
208 assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
210 GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
211 &all.a[traverse.v][k]) );
212 note_computed(traverse.v, k);
215 gsl_interp_accel_free(accel);
216 gsl_interp_free(interp);
218 fprintf(stderr," done\n");
221 static void interpolate(void) {
224 for (y=0; y<(Y+1)/2; y+=INC) {
225 interpolate_line(y<<YSHIFT, 0, OX*2+1);
227 for (x=0; x<X; x+=INC) {
228 interpolate_line( x, 5, OY);
229 interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
234 static void write_output(void) {
235 int x, y, v, c, bad=0;
236 for (y=0; y<Y; y++) {
237 fprintf(stderr,"checking y=%d ", y);
238 for (x=0; x<X; x++) {
240 fprintf(stderr," %02x",v);
241 c= computed_count[v];
242 if (c==D3) fputc('#',stderr);
243 else if (c==-1) fputc('*',stderr);
244 else { fprintf(stderr,"!%d",c); bad++; }
249 if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
250 fprintf(stderr,"interpolated\n");
253 int main(int argc, const char **argv) {
254 if (argc!=1 || isatty(0) || isatty(1))
255 { fputs("usage: interpolate <INPUT.cfm >OUTPUT.cfm\n",stderr); exit(8); }
257 characterise_input();
261 if (fclose(stdout)) diee("fclose stdout");