From 802d7ffd6900d7c17335caa957f1a3752d3f5ccb Mon Sep 17 00:00:00 2001 From: Ian Jackson Date: Sat, 19 Jan 2008 11:10:15 +0000 Subject: [PATCH] interpolation more flexible --- Makefile | 4 +- interpolate.c | 166 ++++++++++++++++++++++++++++++-------------------- 2 files changed, 102 insertions(+), 68 deletions(-) diff --git a/Makefile b/Makefile index 79e8edf..bb869a1 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,8 @@ -VIEWDIMS=33 64 125 +VIEWDIMS=33 64 125 246 TARGETS= minimise primer lumpy.cfm sgtatham.cfm ring.cfm \ - interpolate-64 \ + interpolate-64 interpolate-125 interpolate-246 \ $(addprefix view-, $(VIEWDIMS)) SGTATHAM=sgtatham diff --git a/interpolate.c b/interpolate.c index 778e37c..9225b8f 100644 --- a/interpolate.c +++ b/interpolate.c @@ -4,11 +4,8 @@ #include "mgraph.h" -#define OXBITS (XBITS-1) -#define OX (1< INT_MAX) fail("input file is not reasonable whole number of vertices\n"); - fprintf(stderr,"XBITS=%d XY=%d*%d=%d\n", XBITS, X,Y,N); - fprintf(stderr,"OXBITS=%d OXY=%d*%d=%d\n", OXBITS, OX,OY,ON); - fprintf(stderr,"sizeof(double)=%d st_size=%lu\n", - (int)sizeof(double), (unsigned long)stab.st_size); + printf("XBITS=%d XY=%d*%d=%d", XBITS, X,Y,N); - if (stab.st_size != (ON*sizeof(double)*D3)) - fail("input file wrong size\n"); + oN= stab.st_size / (sizeof(double)*D3); + + for (shift=1, inc=2; + shift%2d", oy,y); - for (ox=0; ox%2d", oy,y); + for (ox=0; ox=X) { y= (Y-1)-y; x-=X; } errno= 0; - ov= (oy << OXBITS) | ox; + ov= (oy << oXBITS) | ox; v= (y << YSHIFT) | x; - fprintf(stderr, " 0%02o->0x%02x", ov, v); + printf(" 0%02o->0x%02x", ov, v); if (fread(all.a[v], sizeof(double), D3, stdin) != D3) diee("\nread input"); computed_count[v]= -1; } - fputc('\n',stderr); + putchar('\n'); } } @@ -94,19 +108,19 @@ static void traverse_next(Traverse *t) { t->v= v2; } -#if 0 - -static void interpolate(void) { +static void interpolate_lin4point(void) { /* four points P Q R S, although P and S may be missing * interpolate in QR finding M. */ int xq,yq,eqr, vp,vq,vm,vr,vs, k; double pqtarg[D3], srtarg[D3]; + if (inc != 2) fail("cannot do >2x lin4point interpolation\n"); + for (eqr=1; eqr!=4; eqr+=5, eqr%=V6) { /* each old edge exactly once */ - fprintf(stderr,"eqr=%d\n",eqr); - for (yq=0; yq>1)&1); xq>1)&1); xq=0) @@ -141,7 +155,7 @@ static void interpolate(void) { K srtarg[k]= all.a[vq][k]; K { - const double alpha= 3.0; + const double alpha= 6; all.a[vm][k]= 0.5 * (all.a[vq][k] + all.a[vr][k]); pqtarg[k]= 0.5 * (pqtarg[k] + all.a[vm][k]); srtarg[k]= 0.5 * (srtarg[k] + all.a[vm][k]); @@ -150,18 +164,14 @@ static void interpolate(void) { note_computed(vm,k); } } - fputc('\n',stderr); + putchar('\n'); } } } -#endif - -#if 1 - -static void interpolate_line(int startvertex, - int direction /* edge number */, - int nmax) { +static void interpolate_gsl_line + (int startvertex, int direction /* edge number */, int nmax, + const gsl_interp_type *aperiodic, const gsl_interp_type *periodic) { double xa[nmax], ya[D3][nmax]; const gsl_interp_type *it; gsl_interp *interp; @@ -170,9 +180,9 @@ static void interpolate_line(int startvertex, Traverse traverse; #define STARTV (traverse.v=startvertex, traverse.e=direction) -#define NEXTV (traverse_next(&traverse), fprintf(stderr,"-%02x",traverse.v)) +#define NEXTV (traverse_next(&traverse), printf("-%02x",traverse.v)) - fprintf(stderr,"interpolate_line 0x%02x,%d nmax=%d ", + printf("interpolate_line 0x%02x,%d nmax=%d ", startvertex,direction,nmax); for (i=0, STARTV; @@ -182,27 +192,27 @@ static void interpolate_line(int startvertex, assert(i < nmax); xa[i]= i; K ya[k][i]= all.a[traverse.v][k]; - fputc('*',stderr); - if (i && traverse.v==startvertex) { fputc('@',stderr); break; } + putchar('*'); + if (i && traverse.v==startvertex) { putchar('@'); break; } i++; NEXTV; if (traverse.v<0) break; NEXTV; } if (traverse.v>=0) { - it= gsl_interp_akima_periodic; + it= periodic; nold= i+1; nnew= i; } else { - it= gsl_interp_akima; + it= aperiodic; nold= i; nnew= i-1; } - fprintf(stderr," n=%d->%d loop=0x%02x", nold,nnew,traverse.v); + printf(" n=%d->%d loop=0x%02x", nold,nnew,traverse.v); K { - fprintf(stderr,"\n k=%d ", k); + printf("\n k=%d ", k); interp= gsl_interp_alloc(it,nold); if (!interp) diee("gsl_interp_alloc"); accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc"); @@ -212,7 +222,7 @@ static void interpolate_line(int startvertex, i=0); NEXTV; assert(traverse.v>=0); - fputc('#',stderr); + putchar('#'); GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel, &all.a[traverse.v][k]) ); note_computed(traverse.v, k); @@ -221,49 +231,73 @@ static void interpolate_line(int startvertex, gsl_interp_accel_free(accel); gsl_interp_free(interp); } - fprintf(stderr," done\n"); + printf(" done\n"); } -static void interpolate(void) { +static void interpolate_gsl(const gsl_interp_type *aperiodic, + const gsl_interp_type *periodic) { int x,y; - for (y=0; y<(Y+1)/2; y+=INC) { - interpolate_line(y<>1)&1), 0, OX*2+1); + for (y=0; y<(Y+1)/2; y+=inc) { + interpolate_gsl_line(y<>1)&1), 0, oX*2+1, + aperiodic, periodic); } - for (x=0; xOUTPUT.cfm\n",stderr); exit(8); } + if (argc!=4 || + strncmp(argv[1],"-a",2) || strlen(argv[1])!=3 || + argv[2][0]=='-' || + strncmp(argv[3],"-o",2)) { + fputs("usage: interpolate -aK INPUT.cfm -oOUTPUT.cfm\n",stderr); + exit(8); + } + if (!freopen(argv[2],"rb",stdin)) diee("open input"); characterise_input(); read_input(); - interpolate(); - write_output(); - if (fclose(stdout)) diee("fclose stdout"); + + switch (argv[1][2]) { + case '4': interpolate_lin4point(); break; + case 'a': interpolate_gsl(gsl_interp_akima, + gsl_interp_akima_periodic); break; + case 'c': interpolate_gsl(gsl_interp_cspline, + gsl_interp_cspline_periodic); break; + default: + fail("unknown interpolation algorithm\n"); + } + + write_output(argv[3]+2); return 0; } -- 2.30.2