home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / RADIANCE / SRC / GEN / GENSURF.C < prev    next >
C/C++ Source or Header  |  1993-10-07  |  13KB  |  577 lines
  1. #ifndef lint
  2. static char SCCSid[] = "@(#)gensurf.c 2.3 2/5/92 LBL";
  3. #endif
  4.  
  5. /* Copyright (c) 1989 Regents of the University of California */
  6.  
  7. /*
  8.  *  gensurf.c - program to generate functional surfaces
  9.  *
  10.  *    Parametric functions x(s,t), y(s,t) and z(s,t)
  11.  *  specify the surface, which is tesselated into an m by n
  12.  *  array of paired triangles.
  13.  *    The surface normal is defined by the right hand
  14.  *  rule applied to (s,t).
  15.  *
  16.  *    4/3/87
  17.  */
  18.  
  19. #include  "standard.h"
  20.  
  21. char  XNAME[] =        "X`SYS`";        /* x function name */
  22. char  YNAME[] =        "Y`SYS`";        /* y function name */
  23. char  ZNAME[] =        "Z`SYS`";        /* z function name */
  24.  
  25. #define  ABS(x)        ((x)>=0 ? (x) : -(x))
  26.  
  27. #define  pvect(p)    printf(vformat, (p)[0], (p)[1], (p)[2])
  28.  
  29. char  vformat[] = "%15.9g %15.9g %15.9g\n";
  30. char  tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
  31. char  texname[] = "Phong";
  32.  
  33. int  smooth = 0;        /* apply smoothing? */
  34.  
  35. char  *modname, *surfname;
  36.  
  37.                 /* recorded data flags */
  38. #define  HASBORDER    01
  39. #define  TRIPLETS    02
  40.                 /* a data structure */
  41. struct {
  42.     int    flags;            /* data type */
  43.     short    m, n;            /* number of s and t values */
  44.     FLOAT    *data;            /* the data itself, s major sort */
  45. } datarec;            /* our recorded data */
  46.  
  47. double  l_hermite(), l_bezier(), l_bspline(), l_dataval();
  48. extern double  funvalue(), argument();
  49.  
  50. typedef struct {
  51.     FVECT  p;    /* vertex position */
  52.     FVECT  n;    /* average normal */
  53. } POINT;
  54.  
  55.  
  56. main(argc, argv)
  57. int  argc;
  58. char  *argv[];
  59. {
  60.     extern long    eclock;
  61.     POINT  *row0, *row1, *row2, *rp;
  62.     int  i, j, m, n;
  63.     char  stmp[256];
  64.  
  65.     varset("PI", ':', PI);
  66.     funset("hermite", 5, ':', l_hermite);
  67.     funset("bezier", 5, ':', l_bezier);
  68.     funset("bspline", 5, ':', l_bspline);
  69.  
  70.     if (argc < 8)
  71.         goto userror;
  72.  
  73.     for (i = 8; i < argc; i++)
  74.         if (!strcmp(argv[i], "-e"))
  75.             scompile(argv[++i], NULL, 0);
  76.         else if (!strcmp(argv[i], "-f"))
  77.             fcompile(argv[++i]);
  78.         else if (!strcmp(argv[i], "-s"))
  79.             smooth++;
  80.         else
  81.             goto userror;
  82.  
  83.     modname = argv[1];
  84.     surfname = argv[2];
  85.     m = atoi(argv[6]);
  86.     n = atoi(argv[7]);
  87.     if (m <= 0 || n <= 0)
  88.         goto userror;
  89.     if (!strcmp(argv[5], "-") || access(argv[5], 4) == 0) {    /* file? */
  90.         funset(ZNAME, 2, ':', l_dataval);
  91.         if (!strcmp(argv[5],argv[3]) && !strcmp(argv[5],argv[4])) {
  92.             loaddata(argv[5], m, n, 3);
  93.             funset(XNAME, 2, ':', l_dataval);
  94.             funset(YNAME, 2, ':', l_dataval);
  95.         } else {
  96.             loaddata(argv[5], m, n, 1);
  97.             sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
  98.             scompile(stmp, NULL, 0);
  99.             sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
  100.             scompile(stmp, NULL, 0);
  101.         }
  102.     } else {
  103.         sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
  104.         scompile(stmp, NULL, 0);
  105.         sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
  106.         scompile(stmp, NULL, 0);
  107.         sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
  108.         scompile(stmp, NULL, 0);
  109.     }
  110.     row0 = (POINT *)malloc((n+3)*sizeof(POINT));
  111.     row1 = (POINT *)malloc((n+3)*sizeof(POINT));
  112.     row2 = (POINT *)malloc((n+3)*sizeof(POINT));
  113.     if (row0 == NULL || row1 == NULL || row2 == NULL) {
  114.         fprintf(stderr, "%s: out of memory\n", argv[0]);
  115.         quit(1);
  116.     }
  117.     row0++; row1++; row2++;
  118.                         /* print header */
  119.     printhead(argc, argv);
  120.     eclock = 0;
  121.                         /* initialize */
  122.     comprow(-1.0/m, row0, n);
  123.     comprow(0.0, row1, n);
  124.     comprow(1.0/m, row2, n);
  125.     compnorms(row0, row1, row2, n);
  126.                         /* for each row */
  127.     for (i = 0; i < m; i++) {
  128.                         /* compute next row */
  129.         rp = row0;
  130.         row0 = row1;
  131.         row1 = row2;
  132.         row2 = rp;
  133.         comprow((double)(i+2)/m, row2, n);
  134.         compnorms(row0, row1, row2, n);
  135.  
  136.         for (j = 0; j < n; j++) {
  137.                             /* put polygons */
  138.             if ((i+j) & 1)
  139.                 putsquare(&row0[j], &row1[j],
  140.                         &row0[j+1], &row1[j+1]);
  141.             else
  142.                 putsquare(&row1[j], &row1[j+1],
  143.                         &row0[j], &row0[j+1]);
  144.         }
  145.     }
  146.  
  147.     quit(0);
  148.  
  149. userror:
  150.     fprintf(stderr, "Usage: %s material name ", argv[0]);
  151.     fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-e expr][-f file]\n");
  152.     quit(1);
  153. }
  154.  
  155.  
  156. loaddata(file, m, n, pointsize)        /* load point data from file */
  157. char  *file;
  158. int  m, n;
  159. int  pointsize;
  160. {
  161.     extern char  *fgetword();
  162.     FILE  *fp;
  163.     char  word[64];
  164.     register int  size;
  165.     register FLOAT  *dp;
  166.  
  167.     datarec.flags = HASBORDER;        /* assume border values */
  168.     datarec.m = m+1;
  169.     datarec.n = n+1;
  170.     size = datarec.m*datarec.n*pointsize;
  171.     if (pointsize == 3)
  172.         datarec.flags |= TRIPLETS;
  173.     dp = (FLOAT *)malloc(size*sizeof(FLOAT));
  174.     if ((datarec.data = dp) == NULL) {
  175.         fputs("Out of memory\n", stderr);
  176.         exit(1);
  177.     }
  178.     if (!strcmp(file, "-")) {
  179.         file = "<stdin>";
  180.         fp = stdin;
  181.     } else if ((fp = fopen(file, "r")) == NULL) {
  182.         fputs(file, stderr);
  183.         fputs(": cannot open\n", stderr);
  184.         exit(1);
  185.     }
  186.     while (size > 0 && fgetword(word, sizeof(word), fp) != NULL) {
  187.         if (!isflt(word)) {
  188.             fprintf(stderr, "%s: garbled data value: %s\n",
  189.                     file, word);
  190.             exit(1);
  191.         }
  192.         *dp++ = atof(word);
  193.         size--;
  194.     }
  195.     if (size == (m+n+1)*pointsize) {    /* no border after all */
  196.         dp = (FLOAT *)realloc((char *)datarec.data,
  197.                 m*n*pointsize*sizeof(FLOAT));
  198.         if (dp != NULL)
  199.             datarec.data = dp;
  200.         datarec.flags &= ~HASBORDER;
  201.         datarec.m = m;
  202.         datarec.n = n;
  203.         size = 0;
  204.     }
  205.     if (datarec.m < 2 || datarec.n < 2 || size != 0 ||
  206.             fgetword(word, sizeof(word), fp) != NULL) {
  207.         fputs(file, stderr);
  208.         fputs(": bad number of data points\n", stderr);
  209.         exit(1);
  210.     }
  211.     fclose(fp);
  212. }
  213.  
  214.  
  215. double
  216. l_dataval(nam)                /* return recorded data value */
  217. char  *nam;
  218. {
  219.     double  u, v;
  220.     register int  i, j;
  221.     register FLOAT  *dp;
  222.     double  d00, d01, d10, d11;
  223.                         /* compute coordinates */
  224.     u = argument(1); v = argument(2);
  225.     if (datarec.flags & HASBORDER) {
  226.         i = u *= datarec.m-1;
  227.         j = v *= datarec.n-1;
  228.     } else {
  229.         i = u = u*datarec.m - .5;
  230.         j = v = v*datarec.n - .5;
  231.     }
  232.     if (i < 0) i = 0;
  233.     else if (i > datarec.m-2) i = datarec.m-2;
  234.     if (j < 0) j = 0;
  235.     else if (j > datarec.n-2) j = datarec.n-2;
  236.                         /* compute value */
  237.     if (datarec.flags & TRIPLETS) {
  238.         dp = datarec.data + 3*(j*datarec.m + i);
  239.         if (nam == ZNAME)
  240.             dp += 2;
  241.         else if (nam == YNAME)
  242.             dp++;
  243.         d00 = dp[0]; d01 = dp[3];
  244.         dp += 3*datarec.m;
  245.         d10 = dp[0]; d11 = dp[3];
  246.     } else {
  247.         dp = datarec.data + j*datarec.m + i;
  248.         d00 = dp[0]; d01 = dp[1];
  249.         dp += datarec.m;
  250.         d10 = dp[0]; d11 = dp[1];
  251.     }
  252.                         /* bilinear interpolation */
  253.     return((j+1-v)*((i+1-u)*d00+(u-i)*d01)+(v-j)*((i+1-u)*d10+(u-i)*d11));
  254. }
  255.  
  256.  
  257. putsquare(p0, p1, p2, p3)        /* put out a square */
  258. POINT  *p0, *p1, *p2, *p3;
  259. {
  260.     static int  nout = 0;
  261.     FVECT  norm[4];
  262.     int  axis;
  263.     FVECT  v1, v2, vc1, vc2;
  264.     int  ok1, ok2;
  265.                     /* compute exact normals */
  266.     fvsum(v1, p1->p, p0->p, -1.0);
  267.     fvsum(v2, p2->p, p0->p, -1.0);
  268.     fcross(vc1, v1, v2);
  269.     ok1 = normalize(vc1) != 0.0;
  270.     fvsum(v1, p2->p, p3->p, -1.0);
  271.     fvsum(v2, p1->p, p3->p, -1.0);
  272.     fcross(vc2, v1, v2);
  273.     ok2 = normalize(vc2) != 0.0;
  274.     if (!(ok1 | ok2))
  275.         return;
  276.                     /* compute normal interpolation */
  277.     axis = norminterp(norm, p0, p1, p2, p3);
  278.  
  279.                     /* put out quadrilateral? */
  280.     if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
  281.         printf("\n%s ", modname);
  282.         if (axis != -1) {
  283.             printf("texfunc %s\n", texname);
  284.             printf(tsargs);
  285.             printf("0\n13\t%d\n", axis);
  286.             pvect(norm[0]);
  287.             pvect(norm[1]);
  288.             pvect(norm[2]);
  289.             fvsum(v1, norm[3], vc1, -0.5);
  290.             fvsum(v1, v1, vc2, -0.5);
  291.             pvect(v1);
  292.             printf("\n%s ", texname);
  293.         }
  294.         printf("polygon %s.%d\n", surfname, ++nout);
  295.         printf("0\n0\n12\n");
  296.         pvect(p0->p);
  297.         pvect(p1->p);
  298.         pvect(p3->p);
  299.         pvect(p2->p);
  300.         return;
  301.     }
  302.                     /* put out triangles? */
  303.     if (ok1) {
  304.         printf("\n%s ", modname);
  305.         if (axis != -1) {
  306.             printf("texfunc %s\n", texname);
  307.             printf(tsargs);
  308.             printf("0\n13\t%d\n", axis);
  309.             pvect(norm[0]);
  310.             pvect(norm[1]);
  311.             pvect(norm[2]);
  312.             fvsum(v1, norm[3], vc1, -1.0);
  313.             pvect(v1);
  314.             printf("\n%s ", texname);
  315.         }
  316.         printf("polygon %s.%d\n", surfname, ++nout);
  317.         printf("0\n0\n9\n");
  318.         pvect(p0->p);
  319.         pvect(p1->p);
  320.         pvect(p2->p);
  321.     }
  322.     if (ok2) {
  323.         printf("\n%s ", modname);
  324.         if (axis != -1) {
  325.             printf("texfunc %s\n", texname);
  326.             printf(tsargs);
  327.             printf("0\n13\t%d\n", axis);
  328.             pvect(norm[0]);
  329.             pvect(norm[1]);
  330.             pvect(norm[2]);
  331.             fvsum(v2, norm[3], vc2, -1.0);
  332.             pvect(v2);
  333.             printf("\n%s ", texname);
  334.         }
  335.         printf("polygon %s.%d\n", surfname, ++nout);
  336.         printf("0\n0\n9\n");
  337.         pvect(p2->p);
  338.         pvect(p1->p);
  339.         pvect(p3->p);
  340.     }
  341. }
  342.  
  343.  
  344. comprow(s, row, siz)            /* compute row of values */
  345. double  s;
  346. register POINT  *row;
  347. int  siz;
  348. {
  349.     double  st[2];
  350.     int  end;
  351.     register int  i;
  352.     
  353.     if (smooth) {
  354.         i = -1;            /* compute one past each end */
  355.         end = siz+1;
  356.     } else {
  357.         if (s < -FTINY || s > 1.0+FTINY)
  358.             return;
  359.         i = 0;
  360.         end = siz;
  361.     }
  362.     st[0] = s;
  363.     while (i <= end) {
  364.         st[1] = (double)i/siz;
  365.         row[i].p[0] = funvalue(XNAME, 2, st);
  366.         row[i].p[1] = funvalue(YNAME, 2, st);
  367.         row[i].p[2] = funvalue(ZNAME, 2, st);
  368.         i++;
  369.     }
  370. }
  371.  
  372.  
  373. compnorms(r0, r1, r2, siz)        /* compute row of averaged normals */
  374. register POINT  *r0, *r1, *r2;
  375. int  siz;
  376. {
  377.     FVECT  v1, v2;
  378.     register int  i;
  379.  
  380.     if (!smooth)            /* not needed if no smoothing */
  381.         return;
  382.                     /* compute middle points */
  383.     while (siz-- >= 0) {
  384.         fvsum(v1, r2[0].p, r0[0].p, -1.0);
  385.         fvsum(v2, r1[1].p, r1[-1].p, -1.0);
  386.         fcross(r1[0].n, v1, v2);
  387.         normalize(r1[0].n);
  388.         r0++; r1++; r2++;
  389.     }
  390. }
  391.  
  392.  
  393. int
  394. norminterp(resmat, p0, p1, p2, p3)    /* compute normal interpolation */
  395. register FVECT  resmat[4];
  396. POINT  *p0, *p1, *p2, *p3;
  397. {
  398. #define u  ((ax+1)%3)
  399. #define v  ((ax+2)%3)
  400.  
  401.     register int  ax;
  402.     MAT4  eqnmat;
  403.     FVECT  v1;
  404.     register int  i, j;
  405.  
  406.     if (!smooth)            /* no interpolation if no smoothing */
  407.         return(-1);
  408.                     /* find dominant axis */
  409.     VCOPY(v1, p0->n);
  410.     fvsum(v1, v1, p1->n, 1.0);
  411.     fvsum(v1, v1, p2->n, 1.0);
  412.     fvsum(v1, v1, p3->n, 1.0);
  413.     ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
  414.     ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
  415.                     /* assign equation matrix */
  416.     eqnmat[0][0] = p0->p[u]*p0->p[v];
  417.     eqnmat[0][1] = p0->p[u];
  418.     eqnmat[0][2] = p0->p[v];
  419.     eqnmat[0][3] = 1.0;
  420.     eqnmat[1][0] = p1->p[u]*p1->p[v];
  421.     eqnmat[1][1] = p1->p[u];
  422.     eqnmat[1][2] = p1->p[v];
  423.     eqnmat[1][3] = 1.0;
  424.     eqnmat[2][0] = p2->p[u]*p2->p[v];
  425.     eqnmat[2][1] = p2->p[u];
  426.     eqnmat[2][2] = p2->p[v];
  427.     eqnmat[2][3] = 1.0;
  428.     eqnmat[3][0] = p3->p[u]*p3->p[v];
  429.     eqnmat[3][1] = p3->p[u];
  430.     eqnmat[3][2] = p3->p[v];
  431.     eqnmat[3][3] = 1.0;
  432.                     /* invert matrix (solve system) */
  433.     if (!invmat(eqnmat, eqnmat))
  434.         return(-1);            /* no solution */
  435.                     /* compute result matrix */
  436.     for (j = 0; j < 4; j++)
  437.         for (i = 0; i < 3; i++)
  438.             resmat[j][i] =    eqnmat[j][0]*p0->n[i] +
  439.                     eqnmat[j][1]*p1->n[i] +
  440.                     eqnmat[j][2]*p2->n[i] +
  441.                     eqnmat[j][3]*p3->n[i];
  442.     return(ax);
  443.  
  444. #undef u
  445. #undef v
  446. }
  447.  
  448.  
  449. /*
  450.  * invmat - computes the inverse of mat into inverse.  Returns 1
  451.  * if there exists an inverse, 0 otherwise.  It uses Gaussian Elimination
  452.  * method.
  453.  */
  454.  
  455. invmat(inverse,mat)
  456. MAT4  inverse, mat;
  457. {
  458. #define SWAP(a,b,t) (t=a,a=b,b=t)
  459.  
  460.     MAT4  m4tmp;
  461.     register int i,j,k;
  462.     register double temp;
  463.  
  464.     copymat4(m4tmp, mat);
  465.                     /* set inverse to identity */
  466.     for (i = 0; i < 4; i++)
  467.         for (j = 0; j < 4; j++)
  468.             inverse[i][j] = i==j ? 1.0 : 0.0;
  469.  
  470.     for(i = 0; i < 4; i++) {
  471.         /* Look for row with largest pivot and swap rows */
  472.         temp = FTINY; j = -1;
  473.         for(k = i; k < 4; k++)
  474.             if(ABS(m4tmp[k][i]) > temp) {
  475.                 temp = ABS(m4tmp[k][i]);
  476.                 j = k;
  477.                 }
  478.         if(j == -1)    /* No replacing row -> no inverse */
  479.             return(0);
  480.         if (j != i)
  481.             for(k = 0; k < 4; k++) {
  482.                 SWAP(m4tmp[i][k],m4tmp[j][k],temp);
  483.                 SWAP(inverse[i][k],inverse[j][k],temp);
  484.                 }
  485.  
  486.         temp = m4tmp[i][i];
  487.         for(k = 0; k < 4; k++) {
  488.             m4tmp[i][k] /= temp;
  489.             inverse[i][k] /= temp;
  490.             }
  491.         for(j = 0; j < 4; j++) {
  492.             if(j != i) {
  493.                 temp = m4tmp[j][i];
  494.                 for(k = 0; k < 4; k++) {
  495.                     m4tmp[j][k] -= m4tmp[i][k]*temp;
  496.                     inverse[j][k] -= inverse[i][k]*temp;
  497.                     }
  498.                 }
  499.             }
  500.         }
  501.     return(1);
  502.  
  503. #undef SWAP
  504. }
  505.  
  506.  
  507. eputs(msg)
  508. char  *msg;
  509. {
  510.     fputs(msg, stderr);
  511. }
  512.  
  513.  
  514. wputs(msg)
  515. char  *msg;
  516. {
  517.     eputs(msg);
  518. }
  519.  
  520.  
  521. quit(code)
  522. {
  523.     exit(code);
  524. }
  525.  
  526.  
  527. printhead(ac, av)        /* print command header */
  528. register int  ac;
  529. register char  **av;
  530. {
  531.     putchar('#');
  532.     while (ac--) {
  533.         putchar(' ');
  534.         fputs(*av++, stdout);
  535.     }
  536.     putchar('\n');
  537. }
  538.  
  539.  
  540. double
  541. l_hermite()            
  542. {
  543.     double  t;
  544.     
  545.     t = argument(5);
  546.     return(    argument(1)*((2.0*t-3.0)*t*t+1.0) +
  547.         argument(2)*(-2.0*t+3.0)*t*t +
  548.         argument(3)*((t-2.0)*t+1.0)*t +
  549.         argument(4)*(t-1.0)*t*t );
  550. }
  551.  
  552.  
  553. double
  554. l_bezier()
  555. {
  556.     double  t;
  557.  
  558.     t = argument(5);
  559.     return(    argument(1) * (1.+t*(-3.+t*(3.-t))) +
  560.         argument(2) * 3.*t*(1.+t*(-2.+t)) +
  561.         argument(3) * 3.*t*t*(1.-t) +
  562.         argument(4) * t*t*t );
  563. }
  564.  
  565.  
  566. double
  567. l_bspline()
  568. {
  569.     double  t;
  570.  
  571.     t = argument(5);
  572.     return(    argument(1) * (1./6.+t*(-1./2.+t*(1./2.-1./6.*t))) +
  573.         argument(2) * (2./3.+t*t*(-1.+1./2.*t)) +
  574.         argument(3) * (1./6.+t*(1./2.+t*(1./2.-1./2.*t))) +
  575.         argument(4) * (1./6.*t*t*t) );
  576. }
  577.