VIEWDIMS=33 44 55
TARGETS= minimise primer lumpy.cfm sgtatham.cfm ring.cfm \
- $(addprefix view+, $(VIEWDIMS))
+ interpolate-44 \
+ $(addprefix view-, $(VIEWDIMS))
SGTATHAM=sgtatham
CWARNS= -Wall -Wwrite-strings -Wpointer-arith -Werror -Wshadow
ring.cfm: oldmoebius-converter prime.data /dev/null ../moebius/a.out
./$^ -o$@
+view-%: view+%.o mgraph+%.o common.o
+ $(CC) $(CFLAGS) -o $@ $^ $(LIBGSL) -L/usr/X11R6/lib -lX11
-# this ridiculous repetition is due to make being too lame
+interpolate-%: interpolate+%.o mgraph+%.o common.o
+ $(CC) $(CFLAGS) -o $@ $^ $(LIBGSL)
-view+%: view+%.o mgraph+%.o common.o
- $(CC) $(CFLAGS) -o $@ $^ $(LIBGSL) -L/usr/X11R6/lib -lX11
+# this ridiculous repetition is due to make being too lame
-view+%.o: view.c
+view+%.o: view.c
$(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFBITS=$* $< -o $@
-mgraph+%.o: mgraph.c
+mgraph+%.o: mgraph.c
$(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFBITS=$* $< -o $@
-.PRECIOUS: view+%.o mgraph+%.o
+interpolate+%.o: interpolate.c
+ $(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFBITS=$* $< -o $@
+.PRECIOUS: view+%.o mgraph+%.o interpolate+%.o
clean:
rm -f prime.data $(TARGETS)
realclean: clean
rm -f best
+%.d:
+
-include *.d
--- /dev/null
+/*
+ * 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;
+}