Back to index

radiance  4R0+20100331
gensurf.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: gensurf.c,v 2.18 2007/10/18 03:59:54 greg Exp $";
00003 #endif
00004 /*
00005  *  gensurf.c - program to generate functional surfaces
00006  *
00007  *     Parametric functions x(s,t), y(s,t) and z(s,t)
00008  *  specify the surface, which is tesselated into an m by n
00009  *  array of paired triangles.
00010  *     The surface normal is defined by the right hand
00011  *  rule applied to (s,t).
00012  *
00013  *     4/3/87
00014  *
00015  *     4/16/02       Added conditional vertex output
00016  */
00017 
00018 #include  "standard.h"
00019 
00020 #include  "paths.h"
00021 #include  "resolu.h"
00022 #include  "rterror.h"
00023 #include  "calcomp.h"
00024 
00025 char  XNAME[] =             "X`SYS";             /* x function name */
00026 char  YNAME[] =             "Y`SYS";             /* y function name */
00027 char  ZNAME[] =             "Z`SYS";             /* z function name */
00028 
00029 char  VNAME[] =      "valid";             /* valid vertex name */
00030 
00031 #define  ABS(x)             ((x)>=0 ? (x) : -(x))
00032 
00033 #define  ZEROVECT(v) (DOT(v,v) <= FTINY*FTINY)
00034 
00035 #define  pvect(p)    printf(vformat, (p)[0], (p)[1], (p)[2])
00036 
00037 char  vformat[] = "%18.12g %18.12g %18.12g\n";
00038 char  tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
00039 char  texname[] = "Phong";
00040 
00041 int  smooth = 0;            /* apply smoothing? */
00042 int  objout = 0;            /* output .OBJ format? */
00043 
00044 char  *modname, *surfname;
00045 
00046                             /* recorded data flags */
00047 #define  HASBORDER   01
00048 #define  TRIPLETS    02
00049                             /* a data structure */
00050 struct {
00051        int    flags;               /* data type */
00052        short  m, n;                /* number of s and t values */
00053        RREAL  *data;               /* the data itself, s major sort */
00054 } datarec;                  /* our recorded data */
00055 
00056 /* XXX this is redundant with rt/noise3.c, should go to a library */
00057 double  l_hermite(), l_bezier(), l_bspline(), l_dataval();
00058 
00059 typedef struct {
00060        int  valid;   /* point is valid (vertex number) */
00061        FVECT  p;     /* vertex position */
00062        FVECT  n;     /* average normal */
00063        RREAL  uv[2]; /* (u,v) position */
00064 } POINT;
00065 
00066 
00067 void loaddata(char *file, int m, int n, int pointsize);
00068 double l_dataval(char *nam);
00069 void putobjrow(POINT *rp, int n);
00070 void putsquare(POINT *p0, POINT *p1, POINT *p2, POINT *p3);
00071 void comprow(double s, POINT *row, int siz);
00072 void compnorms(POINT *r0, POINT *r1, POINT *r2, int siz);
00073 int norminterp(FVECT resmat[4], POINT *p0, POINT *p1, POINT *p2, POINT *p3);
00074 
00075 
00076 int
00077 main(argc, argv)
00078 int  argc;
00079 char  *argv[];
00080 {
00081        POINT  *row0, *row1, *row2, *rp;
00082        int  i, j, m, n;
00083        char  stmp[256];
00084 
00085        varset("PI", ':', PI);
00086        funset("hermite", 5, ':', l_hermite);
00087        funset("bezier", 5, ':', l_bezier);
00088        funset("bspline", 5, ':', l_bspline);
00089 
00090        if (argc < 8)
00091               goto userror;
00092 
00093        for (i = 8; i < argc; i++)
00094               if (!strcmp(argv[i], "-e"))
00095                      scompile(argv[++i], NULL, 0);
00096               else if (!strcmp(argv[i], "-f"))
00097                      fcompile(argv[++i]);
00098               else if (!strcmp(argv[i], "-s"))
00099                      smooth++;
00100               else if (!strcmp(argv[i], "-o"))
00101                      objout++;
00102               else
00103                      goto userror;
00104 
00105        modname = argv[1];
00106        surfname = argv[2];
00107        m = atoi(argv[6]);
00108        n = atoi(argv[7]);
00109        if (m <= 0 || n <= 0)
00110               goto userror;
00111        if (!strcmp(argv[5], "-") || access(argv[5], 4) == 0) { /* file? */
00112               funset(ZNAME, 2, ':', l_dataval);
00113               if (!strcmp(argv[5],argv[3]) && !strcmp(argv[5],argv[4])) {
00114                      loaddata(argv[5], m, n, 3);
00115                      funset(XNAME, 2, ':', l_dataval);
00116                      funset(YNAME, 2, ':', l_dataval);
00117               } else {
00118                      loaddata(argv[5], m, n, 1);
00119                      sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
00120                      scompile(stmp, NULL, 0);
00121                      sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
00122                      scompile(stmp, NULL, 0);
00123               }
00124        } else {
00125               sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
00126               scompile(stmp, NULL, 0);
00127               sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
00128               scompile(stmp, NULL, 0);
00129               sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
00130               scompile(stmp, NULL, 0);
00131        }
00132        row0 = (POINT *)malloc((n+3)*sizeof(POINT));
00133        row1 = (POINT *)malloc((n+3)*sizeof(POINT));
00134        row2 = (POINT *)malloc((n+3)*sizeof(POINT));
00135        if (row0 == NULL || row1 == NULL || row2 == NULL) {
00136               fprintf(stderr, "%s: out of memory\n", argv[0]);
00137               quit(1);
00138        }
00139        row0++; row1++; row2++;
00140                                           /* print header */
00141        fputs("# ", stdout);
00142        printargs(argc, argv, stdout);
00143        eclock = 0;
00144                                           /* initialize */
00145        comprow(-1.0/m, row0, n);
00146        comprow(0.0, row1, n);
00147        comprow(1.0/m, row2, n);
00148        compnorms(row0, row1, row2, n);
00149        if (objout) {
00150               printf("\nusemtl %s\n\n", modname);
00151               putobjrow(row1, n);
00152        }
00153                                           /* for each row */
00154        for (i = 0; i < m; i++) {
00155                                           /* compute next row */
00156               rp = row0;
00157               row0 = row1;
00158               row1 = row2;
00159               row2 = rp;
00160               comprow((double)(i+2)/m, row2, n);
00161               compnorms(row0, row1, row2, n);
00162               if (objout)
00163                      putobjrow(row1, n);
00164 
00165               for (j = 0; j < n; j++) {
00166                      int  orient = (j & 1);
00167                                                  /* put polygons */
00168                      if (!(row0[j].valid && row1[j+1].valid))
00169                             orient = 1;
00170                      else if (!(row1[j].valid && row0[j+1].valid))
00171                             orient = 0;
00172                      if (orient)
00173                             putsquare(&row0[j], &row1[j],
00174                                           &row0[j+1], &row1[j+1]);
00175                      else
00176                             putsquare(&row1[j], &row1[j+1],
00177                                           &row0[j], &row0[j+1]);
00178               }
00179        }
00180 
00181        return 0;
00182 
00183 userror:
00184        fprintf(stderr, "Usage: %s material name ", argv[0]);
00185        fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-o][-e expr][-f file]\n");
00186        return 1;
00187 }
00188 
00189 
00190 void
00191 loaddata(            /* load point data from file */
00192        char  *file,
00193        int  m,
00194        int  n,
00195        int  pointsize
00196 )
00197 {
00198        FILE  *fp;
00199        char  word[64];
00200        register int  size;
00201        register RREAL  *dp;
00202 
00203        datarec.flags = HASBORDER;         /* assume border values */
00204        datarec.m = m+1;
00205        datarec.n = n+1;
00206        size = datarec.m*datarec.n*pointsize;
00207        if (pointsize == 3)
00208               datarec.flags |= TRIPLETS;
00209        dp = (RREAL *)malloc(size*sizeof(RREAL));
00210        if ((datarec.data = dp) == NULL) {
00211               fputs("Out of memory\n", stderr);
00212               exit(1);
00213        }
00214        if (!strcmp(file, "-")) {
00215               file = "<stdin>";
00216               fp = stdin;
00217        } else if ((fp = fopen(file, "r")) == NULL) {
00218               fputs(file, stderr);
00219               fputs(": cannot open\n", stderr);
00220               exit(1);
00221        }
00222        while (size > 0 && fgetword(word, sizeof(word), fp) != NULL) {
00223               if (!isflt(word)) {
00224                      fprintf(stderr, "%s: garbled data value: %s\n",
00225                                    file, word);
00226                      exit(1);
00227               }
00228               *dp++ = atof(word);
00229               size--;
00230        }
00231        if (size == (m+n+1)*pointsize) {   /* no border after all */
00232               dp = (RREAL *)realloc((void *)datarec.data,
00233                             m*n*pointsize*sizeof(RREAL));
00234               if (dp != NULL)
00235                      datarec.data = dp;
00236               datarec.flags &= ~HASBORDER;
00237               datarec.m = m;
00238               datarec.n = n;
00239               size = 0;
00240        }
00241        if (datarec.m < 2 || datarec.n < 2 || size != 0 ||
00242                      fgetword(word, sizeof(word), fp) != NULL) {
00243               fputs(file, stderr);
00244               fputs(": bad number of data points\n", stderr);
00245               exit(1);
00246        }
00247        fclose(fp);
00248 }
00249 
00250 
00251 double
00252 l_dataval(                         /* return recorded data value */
00253        char  *nam
00254 )
00255 {
00256        double  u, v;
00257        register int  i, j;
00258        register RREAL  *dp;
00259        double  d00, d01, d10, d11;
00260                                           /* compute coordinates */
00261        u = argument(1); v = argument(2);
00262        if (datarec.flags & HASBORDER) {
00263               i = u *= datarec.m-1;
00264               j = v *= datarec.n-1;
00265        } else {
00266               i = u = u*datarec.m - .5;
00267               j = v = v*datarec.n - .5;
00268        }
00269        if (i < 0) i = 0;
00270        else if (i > datarec.m-2) i = datarec.m-2;
00271        if (j < 0) j = 0;
00272        else if (j > datarec.n-2) j = datarec.n-2;
00273                                           /* compute value */
00274        if (datarec.flags & TRIPLETS) {
00275               dp = datarec.data + 3*(j*datarec.m + i);
00276               if (nam == ZNAME)
00277                      dp += 2;
00278               else if (nam == YNAME)
00279                      dp++;
00280               d00 = dp[0]; d01 = dp[3];
00281               dp += 3*datarec.m;
00282               d10 = dp[0]; d11 = dp[3];
00283        } else {
00284               dp = datarec.data + j*datarec.m + i;
00285               d00 = dp[0]; d01 = dp[1];
00286               dp += datarec.m;
00287               d10 = dp[0]; d11 = dp[1];
00288        }
00289                                           /* bilinear interpolation */
00290        return((j+1-v)*((i+1-u)*d00+(u-i)*d01)+(v-j)*((i+1-u)*d10+(u-i)*d11));
00291 }
00292 
00293 
00294 void
00295 putobjrow(                  /* output vertex row to .OBJ */
00296        register POINT  *rp,
00297        int  n
00298 )
00299 {
00300        static int    nverts = 0;
00301 
00302        for ( ; n-- >= 0; rp++) {
00303               if (!rp->valid)
00304                      continue;
00305               fputs("v ", stdout);
00306               pvect(rp->p);
00307               if (smooth && !ZEROVECT(rp->n))
00308                      printf("\tvn %.9g %.9g %.9g\n",
00309                                    rp->n[0], rp->n[1], rp->n[2]);
00310               printf("\tvt %.9g %.9g\n", rp->uv[0], rp->uv[1]);
00311               rp->valid = ++nverts;
00312        }
00313 }
00314 
00315 
00316 void
00317 putsquare(           /* put out a square */
00318        POINT *p0,
00319        POINT *p1,
00320        POINT *p2,
00321        POINT *p3
00322 )
00323 {
00324        static int  nout = 0;
00325        FVECT  norm[4];
00326        int  axis;
00327        FVECT  v1, v2, vc1, vc2;
00328        int  ok1, ok2;
00329                                    /* compute exact normals */
00330        ok1 = (p0->valid && p1->valid && p2->valid);
00331        if (ok1) {
00332               VSUB(v1, p1->p, p0->p);
00333               VSUB(v2, p2->p, p0->p);
00334               fcross(vc1, v1, v2);
00335               ok1 = (normalize(vc1) != 0.0);
00336        }
00337        ok2 = (p1->valid && p2->valid && p3->valid);
00338        if (ok2) {
00339               VSUB(v1, p2->p, p3->p);
00340               VSUB(v2, p1->p, p3->p);
00341               fcross(vc2, v1, v2);
00342               ok2 = (normalize(vc2) != 0.0);
00343        }
00344        if (!(ok1 | ok2))
00345               return;
00346        if (objout) {               /* output .OBJ faces */
00347               int    p0n=0, p1n=0, p2n=0, p3n=0;
00348               if (smooth) {
00349                      if (!ZEROVECT(p0->n))
00350                             p0n = p0->valid;
00351                      if (!ZEROVECT(p1->n))
00352                             p1n = p1->valid;
00353                      if (!ZEROVECT(p2->n))
00354                             p2n = p2->valid;
00355                      if (!ZEROVECT(p3->n))
00356                             p3n = p3->valid;
00357               }
00358               if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
00359                      printf("f %d/%d/%d %d/%d/%d %d/%d/%d %d/%d/%d\n",
00360                                    p0->valid, p0->valid, p0n,
00361                                    p1->valid, p1->valid, p1n,
00362                                    p3->valid, p3->valid, p3n,
00363                                    p2->valid, p2->valid, p2n);
00364                      return;
00365               }
00366               if (ok1)
00367                      printf("f %d/%d/%d %d/%d/%d %d/%d/%d\n",
00368                                    p0->valid, p0->valid, p0n,
00369                                    p1->valid, p1->valid, p1n,
00370                                    p2->valid, p2->valid, p2n);
00371               if (ok2)
00372                      printf("f %d/%d/%d %d/%d/%d %d/%d/%d\n",
00373                                    p2->valid, p2->valid, p2n,
00374                                    p1->valid, p1->valid, p1n,
00375                                    p3->valid, p3->valid, p3n);
00376               return;
00377        }
00378                                    /* compute normal interpolation */
00379        axis = norminterp(norm, p0, p1, p2, p3);
00380 
00381                                    /* put out quadrilateral? */
00382        if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
00383               printf("\n%s ", modname);
00384               if (axis != -1) {
00385                      printf("texfunc %s\n", texname);
00386                      printf(tsargs);
00387                      printf("0\n13\t%d\n", axis);
00388                      pvect(norm[0]);
00389                      pvect(norm[1]);
00390                      pvect(norm[2]);
00391                      fvsum(v1, norm[3], vc1, -0.5);
00392                      fvsum(v1, v1, vc2, -0.5);
00393                      pvect(v1);
00394                      printf("\n%s ", texname);
00395               }
00396               printf("polygon %s.%d\n", surfname, ++nout);
00397               printf("0\n0\n12\n");
00398               pvect(p0->p);
00399               pvect(p1->p);
00400               pvect(p3->p);
00401               pvect(p2->p);
00402               return;
00403        }
00404                                    /* put out triangles? */
00405        if (ok1) {
00406               printf("\n%s ", modname);
00407               if (axis != -1) {
00408                      printf("texfunc %s\n", texname);
00409                      printf(tsargs);
00410                      printf("0\n13\t%d\n", axis);
00411                      pvect(norm[0]);
00412                      pvect(norm[1]);
00413                      pvect(norm[2]);
00414                      fvsum(v1, norm[3], vc1, -1.0);
00415                      pvect(v1);
00416                      printf("\n%s ", texname);
00417               }
00418               printf("polygon %s.%d\n", surfname, ++nout);
00419               printf("0\n0\n9\n");
00420               pvect(p0->p);
00421               pvect(p1->p);
00422               pvect(p2->p);
00423        }
00424        if (ok2) {
00425               printf("\n%s ", modname);
00426               if (axis != -1) {
00427                      printf("texfunc %s\n", texname);
00428                      printf(tsargs);
00429                      printf("0\n13\t%d\n", axis);
00430                      pvect(norm[0]);
00431                      pvect(norm[1]);
00432                      pvect(norm[2]);
00433                      fvsum(v2, norm[3], vc2, -1.0);
00434                      pvect(v2);
00435                      printf("\n%s ", texname);
00436               }
00437               printf("polygon %s.%d\n", surfname, ++nout);
00438               printf("0\n0\n9\n");
00439               pvect(p2->p);
00440               pvect(p1->p);
00441               pvect(p3->p);
00442        }
00443 }
00444 
00445 
00446 void
00447 comprow(                    /* compute row of values */
00448        double  s,
00449        register POINT  *row,
00450        int  siz
00451 )
00452 {
00453        double  st[2];
00454        int  end;
00455        int  checkvalid;
00456        register int  i;
00457        
00458        if (smooth) {
00459               i = -1;                     /* compute one past each end */
00460               end = siz+1;
00461        } else {
00462               if (s < -FTINY || s > 1.0+FTINY)
00463                      return;
00464               i = 0;
00465               end = siz;
00466        }
00467        st[0] = s;
00468        checkvalid = (fundefined(VNAME) == 2);
00469        while (i <= end) {
00470               st[1] = (double)i/siz;
00471               if (checkvalid && funvalue(VNAME, 2, st) <= 0.0) {
00472                      row[i].valid = 0;
00473                      row[i].p[0] = row[i].p[1] = row[i].p[2] = 0.0;
00474                      row[i].uv[0] = row[i].uv[1] = 0.0;
00475               } else {
00476                      row[i].valid = 1;
00477                      row[i].p[0] = funvalue(XNAME, 2, st);
00478                      row[i].p[1] = funvalue(YNAME, 2, st);
00479                      row[i].p[2] = funvalue(ZNAME, 2, st);
00480                      row[i].uv[0] = st[0];
00481                      row[i].uv[1] = st[1];
00482               }
00483               i++;
00484        }
00485 }
00486 
00487 
00488 void
00489 compnorms(           /* compute row of averaged normals */
00490        register POINT  *r0,
00491        register POINT  *r1,
00492        register POINT  *r2,
00493        int  siz
00494 )
00495 {
00496        FVECT  v1, v2;
00497 
00498        if (!smooth)                /* not needed if no smoothing */
00499               return;
00500                                    /* compute row 1 normals */
00501        while (siz-- >= 0) {
00502               if (!r1[0].valid)
00503                      continue;
00504               if (!r0[0].valid) {
00505                      if (!r2[0].valid) {
00506                             r1[0].n[0] = r1[0].n[1] = r1[0].n[2] = 0.0;
00507                             continue;
00508                      }
00509                      fvsum(v1, r2[0].p, r1[0].p, -1.0);
00510               } else if (!r2[0].valid)
00511                      fvsum(v1, r1[0].p, r0[0].p, -1.0);
00512               else
00513                      fvsum(v1, r2[0].p, r0[0].p, -1.0);
00514               if (!r1[-1].valid) {
00515                      if (!r1[1].valid) {
00516                             r1[0].n[0] = r1[0].n[1] = r1[0].n[2] = 0.0;
00517                             continue;
00518                      }
00519                      fvsum(v2, r1[1].p, r1[0].p, -1.0);
00520               } else if (!r1[1].valid)
00521                      fvsum(v2, r1[0].p, r1[-1].p, -1.0);
00522               else
00523                      fvsum(v2, r1[1].p, r1[-1].p, -1.0);
00524               fcross(r1[0].n, v1, v2);
00525               normalize(r1[0].n);
00526               r0++; r1++; r2++;
00527        }
00528 }
00529 
00530 
00531 int
00532 norminterp(   /* compute normal interpolation */
00533        register FVECT  resmat[4],
00534        POINT  *p0,
00535        POINT  *p1,
00536        POINT  *p2,
00537        POINT  *p3
00538 )
00539 {
00540 #define u  ((ax+1)%3)
00541 #define v  ((ax+2)%3)
00542 
00543        register int  ax;
00544        MAT4  eqnmat;
00545        FVECT  v1;
00546        register int  i, j;
00547 
00548        if (!smooth)                /* no interpolation if no smoothing */
00549               return(-1);
00550                                    /* find dominant axis */
00551        VCOPY(v1, p0->n);
00552        fvsum(v1, v1, p1->n, 1.0);
00553        fvsum(v1, v1, p2->n, 1.0);
00554        fvsum(v1, v1, p3->n, 1.0);
00555        ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
00556        ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
00557                                    /* assign equation matrix */
00558        eqnmat[0][0] = p0->p[u]*p0->p[v];
00559        eqnmat[0][1] = p0->p[u];
00560        eqnmat[0][2] = p0->p[v];
00561        eqnmat[0][3] = 1.0;
00562        eqnmat[1][0] = p1->p[u]*p1->p[v];
00563        eqnmat[1][1] = p1->p[u];
00564        eqnmat[1][2] = p1->p[v];
00565        eqnmat[1][3] = 1.0;
00566        eqnmat[2][0] = p2->p[u]*p2->p[v];
00567        eqnmat[2][1] = p2->p[u];
00568        eqnmat[2][2] = p2->p[v];
00569        eqnmat[2][3] = 1.0;
00570        eqnmat[3][0] = p3->p[u]*p3->p[v];
00571        eqnmat[3][1] = p3->p[u];
00572        eqnmat[3][2] = p3->p[v];
00573        eqnmat[3][3] = 1.0;
00574                                    /* invert matrix (solve system) */
00575        if (!invmat4(eqnmat, eqnmat))
00576               return(-1);                 /* no solution */
00577                                    /* compute result matrix */
00578        for (j = 0; j < 4; j++)
00579               for (i = 0; i < 3; i++)
00580                      resmat[j][i] =       eqnmat[j][0]*p0->n[i] +
00581                                    eqnmat[j][1]*p1->n[i] +
00582                                    eqnmat[j][2]*p2->n[i] +
00583                                    eqnmat[j][3]*p3->n[i];
00584        return(ax);
00585 
00586 #undef u
00587 #undef v
00588 }
00589 
00590 
00591 double
00592 l_hermite(char *nm)
00593 {
00594        double  t;
00595        
00596        t = argument(5);
00597        return(       argument(1)*((2.0*t-3.0)*t*t+1.0) +
00598               argument(2)*(-2.0*t+3.0)*t*t +
00599               argument(3)*((t-2.0)*t+1.0)*t +
00600               argument(4)*(t-1.0)*t*t );
00601 }
00602 
00603 
00604 double
00605 l_bezier(char *nm)
00606 {
00607        double  t;
00608 
00609        t = argument(5);
00610        return(       argument(1) * (1.+t*(-3.+t*(3.-t))) +
00611               argument(2) * 3.*t*(1.+t*(-2.+t)) +
00612               argument(3) * 3.*t*t*(1.-t) +
00613               argument(4) * t*t*t );
00614 }
00615 
00616 
00617 double
00618 l_bspline(char *nm)
00619 {
00620        double  t;
00621 
00622        t = argument(5);
00623        return(       argument(1) * (1./6.+t*(-1./2.+t*(1./2.-1./6.*t))) +
00624               argument(2) * (2./3.+t*t*(-1.+1./2.*t)) +
00625               argument(3) * (1./6.+t*(1./2.+t*(1./2.-1./2.*t))) +
00626               argument(4) * (1./6.*t*t*t) );
00627 }