Return to regrid.c CVS log

Up to [Research Unix] / researchv10no / cmd / view2d

Annotation of researchv10no/cmd/view2d/regrid.c, revision 1.1

1.1     ! root        1: #include <stdio.h>
        !             2: #include "view2d.h"
        !             3: #define NMAX 1030
        !             4: char *progname;
        !             5: 
        !             6: short timewarp;
        !             7: double ts, te;
        !             8: extern Rd2d rd;
        !             9: int verbose;
        !            10: 
        !            11: typedef struct Frame {
        !            12:   double time;
        !            13:   short *p;
        !            14: } Frame;
        !            15: Frame frame[3];
        !            16: Frame *here, *tween, *there, *tframe;
        !            17: 
        !            18: short nx, ny;  /* input grid */
        !            19: short mx, my;  /* output grid */
        !            20: float x[NMAX], y[NMAX];
        !            21: int detrend;
        !            22: 
        !            23: main(argc, argv)
        !            24:   int argc;
        !            25:   char **argv;
        !            26: {
        !            27:   int i, j;
        !            28:   short *in;  /* input values */
        !            29:   short *mid; /* intermediate values */
        !            30:   int nfr;    /* number of frames to be drawn */
        !            31:   int ifr;
        !            32:   double tsi, tei;  /* starting and ending time to be imposed */
        !            33:   int warpi = 0;  /* internal timewarp */
        !            34:   double twtime, timestep, t;
        !            35:   short *a, *b, *c;
        !            36:   float *d;   /* trend plane */
        !            37:   int fd;
        !            38:   int tinterp = 0;
        !            39:   int ulim = 0;
        !            40:   int timebar = 0;
        !            41:   int firstfr = 1;  /* reading first frame? */
        !            42:   float ufmin, ufmax;
        !            43:   char *Malloc();
        !            44: 
        !            45:   timewarp = 0;
        !            46:   detrend = 0;
        !            47:   verbose = 0;
        !            48:   progname = argv[0];
        !            49:   here  = frame+0;
        !            50:   tween = frame+1;
        !            51:   there = frame+2;
        !            52:   nfr = -1; mx = -1; my = -1;   /* flag uninitialized values */
        !            53: 
        !            54:   for(argc--, argv++; *argv; argv++){
        !            55:     if(**argv == '-' ){
        !            56:       switch(argv[0][1]) {
        !            57:       case 'b':
        !            58:        i = sscanf(&argv[0][2], "%d", &timebar);
        !            59:        if(i!=1) timebar = -1;
        !            60:        break;
        !            61:       case 'f':
        !            62:         i = sscanf(&argv[0][2], "%d", &nfr);
        !            63:         tinterp = 1;
        !            64:         break;
        !            65:       case 'm':
        !            66:         i = sscanf(&argv[0][2], "%e, %e", &ufmin, &ufmax);
        !            67:         if(i!=2) error("bad fmin,fmax");
        !            68:         ulim++;
        !            69:         break;
        !            70:       case 'n':
        !            71:         i = sscanf(&argv[0][2], "%hd, %hd", &mx, &my);
        !            72:         if(i == 1) { my = mx; i = 2; }
        !            73:         if((i!=2)||(mx<0)||(mx>NMAX)||(my<0)||(my>NMAX)) error("bad NX, NY");
        !            74:         break;
        !            75:       case 'r':
        !            76:         detrend++;
        !            77:         break;
        !            78:       case 't':
        !            79:         warpi = sscanf(&argv[0][2], "%E, %E", &tsi, &tei);
        !            80:         if(warpi==0) error("bad TS,TE");
        !            81:         break;
        !            82:       case 'v':
        !            83:         verbose++;
        !            84:         break;
        !            85:       }
        !            86:     }else{
        !            87:       if(fd) error("can only read one file");
        !            88:       fd = Open(*argv,0);
        !            89:     }
        !            90:   }
        !            91: 
        !            92:   if(!tinterp){  /* if not interpolating, can use rd.c's -t processing */
        !            93:     timewarp=warpi;
        !            94:     ts = tsi;
        !            95:     if(timewarp==2) te = tei;
        !            96:   }
        !            97:   rd2dh(fd,&nx,&ny);
        !            98:   if(warpi<=0) tsi = ts;
        !            99:   if(warpi<=1) tei = te;
        !           100:   if(ulim){
        !           101:     rd.fmin = ufmin;
        !           102:     rd.fmax = ufmax;
        !           103:     g_rang2();
        !           104:   }
        !           105:   if(nfr==-1) nfr = rd.nfr;
        !           106:   if((nx > NMAX) || (ny > NMAX)) error("n too large");
        !           107:   if(mx==-1){ mx = nx; my = ny; }
        !           108:   if(verbose){
        !           109:     fprintf(stderr,"ts=%g te=%g\n",tsi,tei);
        !           110:     fprintf(stderr,"nx=%d ny=%d nframes=%d\n",nx,ny,rd.nfr);
        !           111:     fprintf(stderr,"mx=%d my=%d\n",mx,my);
        !           112:     fprintf(stderr,"global starting_time=%g ending_time=%g\n",rd.ts,rd.te);
        !           113:     fprintf(stderr,"fmin=%g fmax=%g\n",rd.fmin,rd.fmax);
        !           114:     fprintf(stderr,"pmin=%d pmax=%d\n",rd.pmin,rd.pmax);
        !           115:     fprintf(stderr,"u=%d v=%d\n",rd.u,rd.v);
        !           116:   }
        !           117:   if(tinterp&&(rd.nfr==1)) error("can't interpolate %d frames!",rd.nfr);
        !           118:   if(timebar&&(rd.ts==rd.te)){
        !           119:      fprintf(stderr,"couldn't determine boundary times for timebar\n");
        !           120:      timebar = 0;
        !           121:   }
        !           122:   for(i=0; i<nx; i++) x[i] = i/(nx-1.);
        !           123:   for(j=0; j<ny; j++) y[j] = j/(ny-1.);
        !           124: 
        !           125:   in = (short *)Malloc(nx*ny*sizeof(short));
        !           126:   mid = (short *)Malloc(mx*ny*sizeof(short));
        !           127:   if(timebar){
        !           128:     if(timebar==-1) timebar = my/30;
        !           129:     if(timebar<1) timebar = 1;
        !           130:     i = mx*(my+3*timebar)*sizeof(short);
        !           131:   }else{
        !           132:     i = mx*my*sizeof(short);
        !           133:   }
        !           134:   j = (nx*ny*sizeof(short)+mx*ny*sizeof(short)+3*i)/1000;
        !           135:   if(detrend) j += (8*mx*my)*sizeof(float)/1000;
        !           136:   if(j>900) fprintf(stderr,"using %dkB workspace\n",j);
        !           137:   here->p  = 3*timebar*mx+(short *)Malloc(i);
        !           138:   tween->p = 3*timebar*mx+(short *)Malloc(i);
        !           139:   there->p = 3*timebar*mx+(short *)Malloc(i);
        !           140:   if(detrend) d = (float *)Malloc(mx*my*sizeof(float));
        !           141: 
        !           142:   if(tinterp){
        !           143:     ifr = 1;
        !           144:     twtime = tsi;
        !           145:     timestep = (tei-tsi)/((nfr-1.)+1e-20);
        !           146:     rdframe(here,in,mid);
        !           147:     if(detrend){ fit(here->p,d); rmfit(here->p,d); }
        !           148:     /* should fit first interpolated frame, but that's a nuisance */
        !           149:     rdframe(there,in,mid);
        !           150:     if(detrend) rmfit(there->p,d);
        !           151:     while( ifr<=nfr ){
        !           152:       while( twtime>there->time ){
        !           153:         tframe = there; there = here; here = tframe; /* swap */
        !           154:         if( rdframe(there,in,mid)==0 ) error("unexpected EOF!");
        !           155:         if(detrend) rmfit(there->p,d);
        !           156:       }
        !           157:       t = (twtime-here->time)/((there->time-here->time)+1e-30);
        !           158:       if(verbose)
        !           159:         fprintf(stderr,"time=%g %g(%g,%g)\n",twtime,t,here->time,there->time);
        !           160:       if( (t<0) || (t>1) ){
        !           161:         fprintf(stderr,"here=%g there=%g\n",here->time,there->time);
        !           162:         error("can't happen: extrap t=%g twtime=%g",t,twtime);
        !           163:         }
        !           164:       for(i=mx*my, a=here->p, b=there->p, c=tween->p; i>0; i--){
        !           165:         if( (*a < -BIG)||(*b < -BIG) ){
        !           166:           *c++ = -BIG-1;  a++; b++;
        !           167:         }else{
        !           168:           *c++ = (1-t) * *a++ + t * *b++;
        !           169:         }
        !           170:       }
        !           171:       wrframe(twtime,tween->p,timebar);
        !           172:       ifr++;
        !           173:       twtime = tsi+(ifr-1)*timestep;
        !           174:       if(twtime>te) twtime=te;
        !           175:       if(twtime>rd.te) break;
        !           176:     }
        !           177:   }else{  /* no time interpolation */
        !           178:     while( rdframe(here,in,mid) ){
        !           179:       if(verbose) fprintf(stderr,"time=%g\n",here->time);
        !           180:       if(detrend){
        !           181:         if(firstfr){ fit(here->p,d); firstfr=0; }
        !           182:         rmfit(here->p,d);
        !           183:       }
        !           184:       wrframe(here->time,here->p,timebar);
        !           185:     }
        !           186:   }
        !           187:   exit(0);
        !           188: }
        !           189: 
        !           190: 
        !           191: wrframe(t,p,timebar)
        !           192:   double t;
        !           193:   short p[];
        !           194:   int timebar;
        !           195: {
        !           196:   if(detrend){
        !           197:     view2d(1,mx,my,t,rd.u,rd.v,0,0,0,p);
        !           198:   }else{
        !           199:     if(timebar){
        !           200:       int i,j;
        !           201:       int k = mx*(t-rd.ts)/(rd.te-rd.ts);
        !           202:       short *q;
        !           203:       short space = -BIG-1;  /* black */
        !           204:       short bar = rd.pmin + .8*((int)rd.pmax-rd.pmin);   /* bar color */
        !           205:       short base = rd.pmin;   /* base-bar color */
        !           206:       q = &p[-3*timebar*mx];
        !           207:       for(i=0;i<timebar;i++){
        !           208:         for(j=0;j<k;j++){ *q++ = bar; }
        !           209:         for(j=k;j<mx;j++){ *q++ = base; }
        !           210:       }
        !           211:       for(i=0;i<2*timebar*mx;i++){ *q++ = space; }
        !           212:     }
        !           213:     view2d(1,mx,my+3*timebar,t,rd.u,rd.v,1,rd.pmin,rd.pmax,&p[-3*timebar*mx]);
        !           214:   }
        !           215: }
        !           216: 
        !           217: 
        !           218: int                           /* returns 0 on EOF */
        !           219: rdframe ( h, u, v )
        !           220:   Frame *h;
        !           221:   short *u;
        !           222:   short *v;
        !           223: {
        !           224:   int i, ii, j, k, m;
        !           225:   float s, t;
        !           226:   short *w;
        !           227:   int umin, umax;
        !           228: 
        !           229:  if( rd2di( &h->time, u ) == 0 ) return(0);
        !           230: 
        !           231:  if((mx==nx)&&(my==ny)){
        !           232:   w = h->p;
        !           233:   m=nx*ny;
        !           234:   for(i = 0; i < m; i++){ *w++ = *u++; }
        !           235:  }else{
        !           236:   /* interpolate in x */
        !           237:   ii = 0;
        !           238:   for(i = 0; i < mx; i++){
        !           239:     t = i/(mx-1.);
        !           240:     while( (ii<nx-1) && (t>=x[ii+1]) ) ii++;
        !           241:     if(t == x[ii]) {
        !           242:       for(j = 0; j < ny; j++) {
        !           243:         v[i+j*mx] = u[ii+j*nx];
        !           244:       }
        !           245:     }else{
        !           246:       s = (t-x[ii])/(x[ii+1]-x[ii]);
        !           247:       for(j = 0; j < ny; j++){
        !           248:         k = ii+j*nx;
        !           249:         if( (u[k] < -BIG)||(u[k+1] < -BIG) ){
        !           250:           v[i+j*mx] = -BIG-1;
        !           251:          }else{
        !           252:           v[i+j*mx] = (u[k] + s*(u[k+1]-u[k]));
        !           253:          }
        !           254:       }
        !           255:     }
        !           256:   }
        !           257: 
        !           258:   /* interpolate in y */
        !           259:   w = h->p;
        !           260:   ii = 0;
        !           261:   for(i = 0; i < my; i++){
        !           262:     t = i/(my-1.);
        !           263:     while( (ii<ny-1) && (t>=y[ii+1]) ) ii++;
        !           264:     if(t == y[ii]){
        !           265:       for(j = 0; j < mx; j++){
        !           266:         w[j+i*mx] = v[j+ii*mx];
        !           267:       }
        !           268:     }else{
        !           269:       s = (t-y[ii])/(y[ii+1]-y[ii]);
        !           270:       for(j = 0; j < mx; j++){
        !           271:         k = j+ii*mx;
        !           272:         if( (v[k] < -BIG)||(v[k+mx] < -BIG) ){
        !           273:           w[j+i*mx] = -BIG-1;
        !           274:          }else{
        !           275:           w[j+i*mx] = (v[k] + s*(v[k+mx]-v[k]));
        !           276:          }
        !           277:       }
        !           278:     }
        !           279:   }
        !           280:  }
        !           281:   return(1);
        !           282: }
        !           283: 
        !           284: 
        !           285: fit(p,d)
        !           286:   short *p;
        !           287:   float *d;  /* l2 fit to p */
        !           288: {
        !           289:   float *z;  /* floating copy of first frame data */
        !           290:   float *zz; /* pointer into z */
        !           291:   short *l;  /* shares storage with z */
        !           292:   float *work, plane[3];  /* workspace for tl2fit */
        !           293:   float x, y;
        !           294:   float u2, v2, pow2();
        !           295:   int i, j;
        !           296:   int imx, imy, imxmy;
        !           297:   char *Malloc();
        !           298: 
        !           299:   i = mx*my*sizeof(float);
        !           300:   z        = (float *)Malloc(i);
        !           301:   work     = (float *)Malloc(6*i);
        !           302:   imx = mx; imy = my;  imxmy = mx*my;
        !           303:   u2 = rd.u;  v2 = pow2(rd.v);
        !           304:   for( i=mx*my, zz=z; i>0; i--){
        !           305:     *zz++ = ((*p++)-u2)*v2;
        !           306:   }
        !           307:   tl2fit_(&imx,&imy,&imxmy,z,work,plane);
        !           308:   free(z); free(work);
        !           309:   if(verbose) fprintf(stderr,"l2 plane: %g + %g x + %g y for 0<=x,y<=1\n",
        !           310:     plane[0], plane[1], plane[2]);
        !           311: 
        !           312:   for(j=1; j<=my; j++){
        !           313:     y = (j-1)/(my-1.);
        !           314:     for(i=1; i<=mx; i++){
        !           315:       x = (i-1)/(mx-1.);
        !           316:       *d++ = (plane[0] + x*plane[1] + y*plane[2])/v2 + u2;
        !           317:     }
        !           318:   }
        !           319: }
        !           320: 
        !           321: rmfit(p,d)            /*  p -= d */
        !           322:   short *p;
        !           323:   float *d;
        !           324: {
        !           325:   int i;
        !           326:   float f;
        !           327:   short *p0, *p1;
        !           328:   p0 = p;
        !           329:   p1 = p+mx*my-1;
        !           330:   for(i=mx*my; i>0; i--, p++, d++){
        !           331:     f = *p - *d;
        !           332:     *p = (f>BIG)? BIG: ( (f<-BIG)? -BIG: f );
        !           333:   }
        !           334: }
        !           335: 

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.