|
|
1.1 ! root 1: #include <stdio.h> ! 2: #include <math.h> ! 3: #include "view2d.h" ! 4: #define MAXPTS 4096 ! 5: ! 6: int NX = 16, NY = 16; ! 7: char *progname; ! 8: ! 9: typedef struct Vector ! 10: { ! 11: float x, y, z; ! 12: } Vector; ! 13: ! 14: typedef struct Frame ! 15: { ! 16: Vector data[MAXPTS]; ! 17: double time; ! 18: int npts; ! 19: } Frame; ! 20: Frame frame[40]; ! 21: int interp = 1; ! 22: ! 23: float grid[200][200]; ! 24: float xmin, xmax, ymin, ymax; ! 25: float xstep, ystep; ! 26: short *data; ! 27: int verbose; ! 28: ! 29: main(argc, argv) ! 30: char **argv; ! 31: { ! 32: register i, j; ! 33: int nframes; ! 34: float x; ! 35: ! 36: verbose = 0; ! 37: progname = *argv; ! 38: for(argc--, argv++; **argv == '-'; argv++) ! 39: switch(argv[0][1]) ! 40: { ! 41: case 'n': ! 42: NY = -1; ! 43: sscanf(&argv[0][2], "%d,%d", &NX, &NY); ! 44: if(NY < 0) NY = NX; ! 45: break; ! 46: case 'c': /* piecewise constant */ ! 47: interp = 2; ! 48: break; ! 49: case 's': /* scatterplot on black background */ ! 50: interp = 3; ! 51: break; ! 52: case 'i': ! 53: interp = 0; ! 54: break; ! 55: case 'v': ! 56: verbose = 1; ! 57: break; ! 58: default: ! 59: error("usage: vdata -c -s -i -v -nNX,NY file"); ! 60: } ! 61: if(*argv) ! 62: if(freopen(*argv, "r", stdin)==NULL) ! 63: error("freopen(%s) failed",*argv); ! 64: if((NX>200)||(NY>200)) error("NX, NY must be <= 200"); ! 65: data = (short *)Malloc(NX*NY*sizeof(short)); ! 66: ! 67: if(interp) ! 68: { ! 69: xmax = ymax = -HUGE; ! 70: xmin = ymin = HUGE; ! 71: } ! 72: for(nframes = 0; readf(&frame[nframes]); nframes++){ ! 73: if(nframes>38) error("too many frames"); ! 74: } ! 75: if(nframes<=0) error("empty input"); ! 76: if(verbose)fprintf(stderr, "nx=%d ny=%d nframes=%d\n", NX, NY, nframes); ! 77: ! 78: if(interp) ! 79: { ! 80: xstep = (xmax-xmin)/(NX-1); ! 81: ystep = (ymax-ymin)/(NY-1); ! 82: if(verbose)fprintf(stderr, "%.5g<=x<=%.5g %.5g<=y<=%.5g\n", ! 83: xmin, xmax, ymin, ymax); ! 84: } ! 85: else ! 86: { ! 87: xmin = 0; xstep = 1; xmax = xmin + (NX-1); ! 88: ymin = 0; ystep = 1; ymax = ymin + (NY-1); ! 89: } ! 90: for(i = 0; i < nframes; i++){ ! 91: fit(&frame[i]); ! 92: out(frame[i].time,frame[i].npts); ! 93: } ! 94: exit(0); ! 95: } ! 96: ! 97: readf(s) ! 98: Frame *s; ! 99: { ! 100: register Vector *v; ! 101: register i; ! 102: float f; ! 103: ! 104: if(scanf("%d %E", &s->npts, &s->time) != 2) return(0); ! 105: if(verbose) fprintf(stderr,"%d %g\n",s->npts,s->time); ! 106: if(s->npts>MAXPTS) error("too many data points"); ! 107: if(interp) ! 108: for(v = s->data, i = 0; i < s->npts; i++, v++) ! 109: { ! 110: scanf("%f %f %f", &v->x, &v->y, &v->z); ! 111: if(v->x < xmin) xmin = v->x; ! 112: if(v->x > xmax) xmax = v->x; ! 113: if(v->y < ymin) ymin = v->y; ! 114: if(v->y > ymax) ymax = v->y; ! 115: } ! 116: else ! 117: for(v = s->data, i = 0; i < s->npts; i++, v++) ! 118: { ! 119: scanf("%f", &v->z); ! 120: } ! 121: return(1); ! 122: } ! 123: ! 124: fit(s) ! 125: Frame *s; ! 126: { ! 127: register Vector *v; ! 128: Vector h; ! 129: float g, x, y; ! 130: register i, j, k; ! 131: Vector *c1, *c2, *c3; ! 132: float d, d1, d2, d3; ! 133: ! 134: if( (interp==1) || (interp==2) ) { ! 135: for(j = 0; j < NY; j++) ! 136: for(i = 0; i < NX; i++) { ! 137: x = xmin + i*xstep; ! 138: y = ymin + j*ystep; ! 139: g = grid[i][j]; ! 140: d1 = d2 = d3 = HUGE; ! 141: for(v = s->data, k = 0; k < s->npts; k++, v++) { ! 142: d = (v->x - x)*(v->x - x) + ! 143: (v->y - y)*(v->y - y); ! 144: if(d < d1) { ! 145: c3 = c2; d3 = d2; ! 146: c2 = c1; d2 = d1; ! 147: c1 = v; d1 = d; ! 148: } else if(d < d2) { ! 149: c3 = c2; d3 = d2; ! 150: c2 = v; d2 = d; ! 151: } else if(d < d3) { ! 152: c3 = v; d3 = d; ! 153: } ! 154: } ! 155: /* c1 c2 c3 three closest points */ ! 156: if(interp==1) { ! 157: h.x = x; h.y = y; h.z = g; ! 158: interp_(&h, c1, c2, c3); ! 159: grid[i][j] = h.z; ! 160: } else if(interp==2) ! 161: grid[i][j] = c1->z; ! 162: } ! 163: } else if( interp==3 ) { ! 164: for(j = 0; j < NY; j++) ! 165: for(i = 0; i < NX; i++) ! 166: grid[i][j] = -1e26; ! 167: for(v = s->data, k = 0; k < s->npts; k++, v++) { ! 168: i = (NX-1)*(v->x-xmin)/(xmax-xmin); ! 169: j = (NY-1)*(v->y-ymin)/(ymax-ymin); ! 170: if((i<0)||(i>=NX)||(j<0)||(j>=NY)) ! 171: error("can't: i,j=%d,%d\n",i,j); ! 172: grid[i][j] = v->z; ! 173: } ! 174: } else { /* option -i */ ! 175: for(v = s->data, j = 0; j < NY; j++) ! 176: for(i = 0; i < NX; i++, v++) ! 177: grid[i][j] = v->z; ! 178: } ! 179: } ! 180: ! 181: out(tim,npt) ! 182: double tim; ! 183: int npt; ! 184: { ! 185: register short *d, i, j; ! 186: register float zmax, zmin, zrange, f, v2; ! 187: short u, v; ! 188: extern double pow2(); ! 189: ! 190: zmin = HUGE; ! 191: zmax = -HUGE; ! 192: for(j = 0; j < NY; j++) ! 193: for(i = 0; i < NX; i++){ ! 194: if(grid[i][j] < -1e25 ) continue; ! 195: if(grid[i][j] < zmin) zmin = grid[i][j]; ! 196: if(grid[i][j] > zmax) zmax = grid[i][j]; ! 197: } ! 198: ! 199: if(verbose) fprintf(stderr, "time=%g npts=%d %g<=z<=%g\n", ! 200: tim, npt, zmin, zmax); ! 201: zmax = (zmax<0)?0:zmax; ! 202: zmin = (zmin>0)?0:zmin; ! 203: zrange = zmax - zmin; ! 204: if(zrange<=0) ! 205: v=0, u=0, v2=1, f=0; ! 206: else{ ! 207: f = 2*BIG-2; ! 208: v = log(zrange)/log(2.)-14; ! 209: v2 = pow2(v); ! 210: while( 2*zrange < f*v2 ){ v--; v2 /= 2; } ! 211: f = -(zmin+zmax)/(2*v2); ! 212: u = f; ! 213: if( (f<-BIG)||(f>BIG) ){ ! 214: if(verbose)fprintf(stderr,"zmin=%g zmax=%g u=%d v=%d\n", ! 215: zmin,zmax,u,v); ! 216: error("u out of bounds f=%g",f); ! 217: } ! 218: } ! 219: d = data; ! 220: for(j = 0; j < NY; j++) ! 221: for(i = 0; i < NX; i++) ! 222: if(grid[i][j] < -1e25){ *d++ = -BIG-1; } ! 223: else{ *d++ = f + grid[i][j]/v2; } ! 224: view2d(1, NX, NY, tim, u, v, 0, 0, 0, data); ! 225: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.