Return to regrid.c CVS log

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

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

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