home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / gle / util / fitz / fitmain.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-29  |  6.3 KB  |  283 lines
  1. /*
  2.     Surface fitting program,  creates .z file from data points
  3. */
  4.  
  5. #include <ctype.h>
  6. #include <stdio.h>
  7. #include <stdlib.h>
  8. #include <string.h>
  9. #include <math.h>
  10. #define false (0)
  11. #define true (!false)
  12. #define gprint printf
  13. int load_data(char *s);
  14. int getstr(char *s);
  15. #ifdef unix
  16. #define VAXC 1
  17. #endif
  18. #ifdef VAXC
  19. #define farfree free
  20. #define farmalloc malloc
  21. #else
  22. #include <alloc.h>
  23. #endif
  24. typedef float real;
  25. int idbvip_(long *md, long *ncp, long *ndp, real *xd, real *yd, real *zd
  26.     ,long *nip, real *xi, real *yi, real *zi, long *iwk, real *wk);
  27. int setminmax(double v, double *min, double *max);
  28. int sort_data(int npnts,float *xd,float *yd,float *zd);
  29. int getrange(char *s, double  *min, double  *max, double  *step);
  30. int xyz_alloc;
  31. double getf(void);
  32. float *pntxyz;
  33. FILE *fp;
  34. long npnts;
  35. main()
  36. {
  37.     long md,ncp,ndp;
  38.     /* real xd[8],yd[8],zd[8]; */
  39.     char bfile[80],dfile[80],xrange[80];
  40.     real *xd, *yd, *zd;
  41.     double xmin,xmax,xstep,ymin,ymax,ystep,x,y;
  42.     real *rx, *ry, *rz;
  43.     long nx,ny;
  44.     long nrp;
  45.     /* real xi[888],yi[888],zi[888]; */
  46.     long *iwk;
  47.     real *wk;
  48.     long i,j,k;
  49.     
  50.     ymin = xmin = 10e10;
  51.     xmax = ymax = -xmin;
  52.  
  53.     printf("Data file conaining x y z data ? "); getstr(dfile);
  54.     load_data(dfile);
  55.     printf("Read %ld points from file {%s} \n",npnts/3,dfile);
  56.  
  57.     strcpy(bfile,dfile);
  58.     if (strchr(bfile,'.')!=NULL) *strchr(bfile,'.') = 0;
  59.     strcat(bfile,".z");
  60.     printf("Will write data to {%s} \n",bfile);
  61.     xd = malloc(sizeof(float) * (npnts/3+1));    
  62.     yd = malloc(sizeof(float) * (npnts/3+1));    
  63.     zd = malloc(sizeof(float) * (npnts/3+1));    
  64.     printf("Number of points to use for contouring each point ? [3] "); ncp = getf();
  65.     if (ncp==0) ncp = 3;
  66.     
  67.     for (j=0,i=0; i<npnts; i+=3) {
  68.         xd[j] = pntxyz[i];
  69.         yd[j] = pntxyz[i+1];
  70.         zd[j] = pntxyz[i+2];
  71.         setminmax(xd[j],&xmin,&xmax);
  72.         setminmax(yd[j],&ymin,&ymax);
  73.         j++;
  74.     }
  75.     farfree(pntxyz);
  76.     npnts = npnts/3;
  77.     sort_data(npnts,xd,yd,zd);
  78.  
  79.     for (i=0; i<npnts; i++) {
  80. /*        printf("%g %g %g \n",xd[i],yd[i],zd[i]); */
  81.         if (xd[i]==xd[i+1] && yd[i]==yd[i+1]) {
  82.             printf("Duplicate points %g %g %g \n"
  83.                 ,xd[i],yd[i],zd[i]);
  84.         }
  85.     }
  86.  
  87.     xstep = (xmax-xmin)/15.0;
  88.     ystep = (ymax-ymin)/15.0;
  89.     printf("Range of output x values [%g,%g,%g] ? ",xmin,xmax,xstep); getstr(xrange);
  90.     getrange(xrange,&xmin,&xmax,&xstep);
  91.     printf("Range of output y values [%g,%g,%g] ? ",ymin,ymax,ystep); getstr(xrange);
  92.     getrange(xrange,&ymin,&ymax,&ystep);
  93.  
  94.     xyz_alloc = 0;
  95.     nx = (xmax-xmin)/xstep + 1;
  96.     ny = (ymax-ymin)/ystep + 1;
  97.     rx = malloc(ny * nx * sizeof(float));
  98.     ry = malloc(nx * ny * sizeof(float));
  99.     rz = malloc(nx * ny * sizeof(float));
  100.     
  101.     k = 0;
  102.     for (j=0, y=ymin; j<ny; y+= ystep, j++) {
  103.         for (i=0, x=xmin; i<nx; x+= xstep, i++) {
  104.             rx[k] = x;
  105.             ry[k++] = y;
  106.         }
  107.     }    
  108.     ndp = npnts;
  109.     md = 1;
  110.     nrp = nx*ny;
  111.     k = 0;
  112.  
  113.     i = 27+ncp;
  114.     if (i<31) i = 31;
  115.     i = (i * ndp + nrp) * sizeof(*iwk) ;
  116.     j = 8*ndp*sizeof(*wk);
  117.     iwk = malloc(i);
  118.     wk = malloc(j);
  119.     printf("nx %ld, ny %ld, Work space iwk=%ld bytes wk=%ld bytes \n",nx,ny,i,j);
  120.     if (iwk==NULL || wk == NULL) {
  121.         printf("Unable to allocate enough workspace iwk=%ld wk=%ld  \n",i,j);
  122.         exit(1);
  123.     }
  124.  
  125.     idbvip_(&md,&ncp,&ndp,xd,yd,zd,&nrp,rx,ry,rz,iwk,wk);
  126.  
  127.     fp = fopen(bfile,"w");
  128.     if (fp==NULL) {
  129.         printf("Unable to Z file {%s} \n",bfile);
  130.         exit(1);
  131.     }
  132.  
  133.     fprintf(fp,"! nx %ld ny %ld xmin %g xmax %g ymin %g ymax %g\n",nx,ny,xmin,xmax,ymin,ymax);
  134.     printf("Y = ");
  135.     k = 0;
  136.     for (j=0, y=ymin; j<ny; y+= ystep, j++) {
  137.           printf("%g ",y);
  138.         for (i=0, x=xmin; i<nx; x+= xstep, i++) {
  139.             fprintf(fp,"%g ",rz[k++]);
  140.         }
  141.         fprintf(fp,"\n");
  142.     }
  143.     fclose(fp);
  144.     printf("\n");
  145. }
  146. pnt_alloc(int size)
  147. {
  148.     void *d;
  149.     if (size+10<xyz_alloc) return;
  150.     size = size*2;
  151.     d = farmalloc(size*sizeof(float));
  152.     if (d==NULL) {
  153.         gprint("Unable to allocate storage for data\n");
  154.         exit(1);
  155.     }
  156.     if (xyz_alloc>0) memcpy(d,pntxyz,xyz_alloc*sizeof(float));
  157.     xyz_alloc = size;
  158.     pntxyz = d;
  159. }
  160. static char buff[2001];
  161. load_data(char *fname)
  162. {
  163.     double v;
  164.     char *s;
  165.     FILE *df;
  166.     int i,nd,nc;
  167.     int done_err;
  168.     done_err = false;
  169.  
  170.     pnt_alloc(30);
  171.  
  172.     df = fopen(fname,"r");
  173.     if (df==NULL) {
  174.         printf("Unable to open data file %s \n",fname);
  175.         exit(1);
  176.     }
  177.     nd = 0;
  178.     for (;!feof(df);) {
  179.       if (fgets(buff,2000,df)!=NULL) {
  180.         if (strchr(buff,'!')!=NULL) *strchr(buff,'!') = 0;
  181.         nc = 0;
  182.         s = strtok(buff," \t\n,");
  183.         for (;s!=NULL;) {
  184.             v = atof(s);
  185.             pnt_alloc(nd);
  186.             if (isdigit(*s) || *s=='-' || *s=='+' || *s=='.') {
  187.                 pntxyz[nd++] = v; nc++;
  188.             } else gprint("Not a number {%s} \n",s);
  189.             s = strtok(NULL," \t\n,");
  190.         }
  191.         if (nc>0 && nc!=3) {
  192.             if (!done_err) gprint("Expecting 3 columns in data file, found %d \n {%s}\n (FATAL ERROR)\n",nc,buff);
  193.             exit(1);
  194.         }
  195.       }
  196.     }
  197.     fclose(df);
  198.     npnts = nd;
  199. }
  200. double getf()
  201. {
  202.     char buff[80];
  203.     getstr(buff);
  204.     return atof(buff);
  205. }
  206. getstr(char *s)
  207. {
  208.     char *ss;
  209.     ss = gets(s);
  210.     if (ss==NULL) {
  211.         printf("\nERROR, End of input file \n");
  212.         exit(1);
  213.     }
  214.     ss = strchr(s,'\n');
  215.     if (ss!=NULL) *ss = 0;
  216.     ss = strchr(s,'!');
  217.     if (ss!=NULL) *ss = 0;
  218. }
  219. getrange(char *s, double  *min, double  *max, double  *step)
  220. {
  221.     s = strtok(s,", :\n\t");
  222.     if (s!=NULL) *min = atof(s);
  223.     s = strtok(NULL,", :\n\t");
  224.     if (s!=NULL) *max = atof(s);
  225.     *step = (*max - *min) / 15.0;
  226.     s = strtok(NULL,", :\n\t");
  227.     if (s!=NULL) *step = atof(s);
  228. }
  229. setminmax(double v, double *min, double *max)
  230. {
  231.     if (v< *min) *min = v;
  232.     if (v> *max) *max = v;
  233. }
  234. float *xxx,*yyy,*zzz;
  235. myswap(int i, int j)
  236. {
  237.     static float a;
  238.     a = xxx[i]; xxx[i] = xxx[j]; xxx[j] = a;
  239.     a = yyy[i]; yyy[i] = yyy[j]; yyy[j] = a;
  240.     a = zzz[i]; zzz[i] = zzz[j]; zzz[j] = a;
  241. }
  242. mycmp(int i, float x1, float y1)
  243. {
  244.     static float a;
  245.     if (xxx[i]<x1) return -1;
  246.     if (xxx[i]>x1) return 1;
  247.     if (yyy[i]<y1) return -1;
  248.     if (yyy[i]>y1) return 1;
  249.     return 0;
  250. }
  251. int quick_sort(int nitems, int (*fswap) (int i, int j), int (*fcmp) (int i, float x1, float x2));
  252. sort_data(int npnts,float *xd,float *yd,float *zd)
  253. {
  254.     xxx = xd; yyy = yd; zzz = zd;
  255.     quick_sort(npnts,myswap,mycmp);        
  256. }
  257.  
  258. int (*ffswap) (int i, int j);
  259. int (*ffcmp) (int i, float x1, float x2);
  260. int qquick_sort(int left, int right);
  261. quick_sort(int nitems, int (*fswap) (int i, int j), int (*fcmp) (int i,
  262.     float x1, float x2))
  263. {
  264.     ffswap = fswap;
  265.     ffcmp = fcmp;
  266.     qquick_sort(0,nitems-1);
  267. }
  268. qquick_sort(int left, int right)
  269. {
  270.     int i,j,xx,yy;
  271.     float x1,x2;
  272.     i = left; j = right;
  273.     xx = (left+right)/2;
  274.     x1 = xxx[xx]; x2 = yyy[xx];
  275.     do {
  276.         while((*ffcmp)(i,x1,x2)<0 && i<right) i++;
  277.         while((*ffcmp)(j,x1,x2)>0 && j>left) j--;
  278.         if (i<=j) { (*ffswap)(i,j); i++; j--; }
  279.     } while (i<=j);
  280.     if (left<j) qquick_sort(left,j);
  281.     if (i<right) qquick_sort(i,right);
  282. }
  283.