home *** CD-ROM | disk | FTP | other *** search
/ gondwana.ecr.mu.oz.au/pub/ / Graphics.tar / Graphics / VOGLE.ZIP / SRC / PATCHES.C < prev    next >
C/C++ Source or Header  |  2000-02-11  |  10KB  |  484 lines
  1. #include "vogle.h"
  2.  
  3. static int    ntsegs, tsegs = 10, nusegs, usegs = 10, 
  4.         ntcurves = 10, nucurves = 10, 
  5.         ntiter = 1, nuiter = 1,
  6.         makeprec_called = 0;
  7.  
  8. static Matrix    tbasis, ubasis, et, eu,
  9.         ones =    {
  10.                 {1.0, 1.0, 1.0, 1.0},
  11.                 {1.0, 1.0, 1.0, 1.0},
  12.                 {1.0, 1.0, 1.0, 1.0},
  13.                 {1.0, 1.0, 1.0, 1.0}
  14.             };
  15.             
  16.  
  17. void        premulttensor(), multtensor(), 
  18.         copytensor(), copytensortrans(),
  19.         drpatch(), drcurve(), 
  20.         rpatch(), transformtensor();
  21.  
  22. static    void    makeprec(), replace(), iterate(), extract();
  23.  
  24. static    float    addemup();
  25.  
  26. /*
  27.  * patchbasis
  28.  *
  29.  *    Specify the two basis matrices for a patch
  30.  */
  31. void
  32. patchbasis(tb, ub)
  33.         Matrix     tb, ub;
  34. {
  35.         if(!vdevice.initialised)
  36.                 verror("patchbasis: vogle not initialised");
  37.  
  38.         copymatrix(tbasis, tb);
  39.     copytranspose(ubasis, ub);
  40. }
  41.  
  42. /*
  43.  * patchprecision
  44.  *
  45.  *    Specify the lower limit on the number of line segments used
  46.  * to draw a curve as part of a patch. The actual number used varies
  47.  * with the number of curve segments in the "other" direction.
  48.  */
  49. void
  50. patchprecision(tseg, useg)
  51.         int     tseg, useg;
  52. {
  53.         if (!vdevice.initialised)
  54.                 verror("patchprecision: vogle not initialised");
  55.  
  56.         if(tseg > 0 && useg > 0) {
  57.                 tsegs = tseg;
  58.                 usegs = useg;
  59.         } else
  60.                 verror("patchprecision: number of segments <= 0");
  61.     /*
  62.      * Set up the difference matrices
  63.      */
  64.     makeprec();
  65. }
  66.  
  67. /*
  68.  * patchcurves
  69.  *
  70.  *    Specify the number of curves to be drawn in each direction 
  71.  *    on a patch.
  72.  */
  73. void
  74. patchcurves(nt, nu)
  75.         int     nt, nu;
  76. {
  77.  
  78.         if (!vdevice.initialised)
  79.                 verror("patchcurves: vogle not initialised");
  80.  
  81.         if(nt > 0 && nu > 0) {
  82.                 ntcurves = nt;
  83.                 nucurves = nu;
  84.         } else
  85.                 verror("patchcurves: number of patch curves <= 0");
  86.  
  87.     /*
  88.      * Set up the difference matrices
  89.      */
  90.     makeprec();
  91. }
  92.  
  93. /*
  94.  * makeprec
  95.  *
  96.  *    Makes up the two precision matrices for a patch
  97.  */
  98. static    void
  99. makeprec()
  100. {
  101.     float    n2, n3;
  102.  
  103.     /*
  104.      * Find ntsegs, nusegs, ntiter, nuiter....
  105.      * ie. the actual number of curve segments of a tcurve,
  106.      *     the actual number of curve segments of a ucurve,
  107.      *     and the number of times to iterate in each direction.
  108.      */
  109.  
  110.     ntsegs = tsegs;
  111.     ntiter = ntsegs / (nucurves - 1);
  112.     if (ntsegs > ntiter * (nucurves - 1)) 
  113.         ntsegs = (++ntiter) * (nucurves - 1);
  114.  
  115.     nusegs = usegs;
  116.     nuiter = nusegs / (ntcurves - 1);
  117.     if (nusegs > nuiter * (ntcurves - 1)) 
  118.         nusegs = (++nuiter) * (ntcurves - 1);
  119.  
  120.     /*
  121.      * Doing the t precision matrix.....
  122.      */
  123.     identmatrix(et);
  124.     n2 = (float)(ntsegs * ntsegs);
  125.     n3 = (float)(ntsegs * n2);
  126.  
  127.     et[0][0] = et[2][2] = et[3][3] = 0.0;
  128.     et[1][0] = 1.0 / n3;
  129.     et[1][1] = 1.0 / n2;
  130.     et[2][0] = et[3][0] = 6.0 / n3;
  131.     et[2][1] = 2.0 / n2;
  132.     et[1][2] = 1.0 / (float)ntsegs;
  133.     et[0][3] = 1.0;
  134.  
  135.     /*
  136.      * Make the Transpose of eu
  137.      */
  138.     identmatrix(eu);
  139.     n2 = (float)(nusegs * nusegs);
  140.     n3 = (float)(nusegs * n2);
  141.  
  142.     eu[0][0] = eu[2][2] = eu[3][3] = 0.0;
  143.     eu[0][1] = 1.0 / n3;
  144.     eu[1][1] = 1.0 / n2;
  145.     eu[0][2] = eu[0][3] = 6.0 / n3;
  146.     eu[1][2] = 2.0 / n2;
  147.     eu[2][1] = 1.0 / (float)nusegs;
  148.     eu[3][0] = 1.0;
  149.  
  150.     makeprec_called = 1;
  151. }
  152.  
  153.  
  154. /*
  155.  * patch
  156.  *
  157.  *    Draws a bicubic patch. (ie. all the w coords a 1 and the
  158.  *    basis matrices don't change that)
  159.  */
  160. void
  161. patch(geomx, geomy, geomz)
  162.     Matrix    geomx, geomy, geomz;
  163. {
  164.     rpatch(geomx, geomy, geomz, ones);
  165. }
  166.  
  167. /*
  168.  * rpatch
  169.  *
  170.  *    Draws rational bicubic patches.
  171.  *
  172.  *    Reference: J. H. Clark, Parametric Curves, Surfaces and volumes in 
  173.  *    computer graphics and computer aided Geometric Design.
  174.  *    Technical report No. 221, Nov 1981.
  175.  *    Computer Systems Lab. Dept's of Elecrical Eng. and Computer Science,
  176.  *    Standford University, Standford, California 94305.
  177.  */
  178. void
  179. rpatch(geomx, geomy, geomz, geomw)
  180.         Matrix  geomx, geomy, geomz, geomw;
  181. {
  182.  
  183.     Tensor    S, R;
  184.         Matrix    tmp, tmp2;
  185.     float    xlast, ylast, zlast;
  186.         int    i, j, sync;
  187.     Token    *tok;
  188.  
  189.         if (!vdevice.initialised)
  190.                 verror("patch: vogle not initialised");
  191.  
  192.     /* 
  193.      *  Form S = et . tbasis . Gtensor . ubasisT . euT
  194.      */
  195.  
  196.     if (!makeprec_called)
  197.         makeprec();
  198.  
  199.     mult4x4(tmp, et, tbasis);
  200.     mult4x4(tmp2, ubasis, eu);
  201.  
  202.     /*
  203.      * Load the geometry matrices into S.
  204.      */
  205.     for (i = 0; i < 4; i++)
  206.         for (j = 0; j < 4; j++) {
  207.             S[0][i][j] = geomx[i][j];
  208.             S[1][i][j] = geomy[i][j];
  209.             S[2][i][j] = geomz[i][j];
  210.             S[3][i][j] = geomw[i][j];
  211.         }
  212.  
  213.     premulttensor(R, tbasis, S);
  214.     multtensor(S, ubasis, R);
  215.  
  216.     /*
  217.      * Find the last point on the curve.
  218.      */
  219.     xlast = addemup(S[0]);
  220.     ylast = addemup(S[1]);
  221.     zlast = addemup(S[2]);
  222.  
  223.     /*
  224.       * Multiply the precision matrices in.
  225.      */
  226.     premulttensor(R, et, S);
  227.     multtensor(S, eu, R);
  228.  
  229.     if (vdevice.inobject) {
  230.         tok = newtokens(74);
  231.         tok[0].i = RPATCH;
  232.         tok[1].f = xlast;
  233.         tok[2].f = ylast;
  234.         tok[3].f = zlast;
  235.         tok[4].i = ntcurves;
  236.         tok[5].i = nucurves;
  237.         tok[6].i = ntsegs;
  238.         tok[7].i = nusegs;
  239.         tok[8].i = ntiter;
  240.         tok[9].i = nuiter;
  241.  
  242.         tok += 10;
  243.         for (i = 0; i < 4; i++)
  244.             for (j = 0; j < 4; j++) {
  245.                 (tok++)->f = S[0][i][j];
  246.                 (tok++)->f = S[1][i][j];
  247.                 (tok++)->f = S[2][i][j];
  248.                 (tok++)->f = S[3][i][j];
  249.             }
  250.  
  251.         return;
  252.     }
  253.  
  254.     /*
  255.      * Multiply by the current transformation....
  256.      */
  257.     transformtensor(S, vdevice.transmat->m);
  258.  
  259.     /*
  260.      * Draw the patch....
  261.      */
  262.     if ((sync = vdevice.sync))
  263.         vdevice.sync = 0;
  264.  
  265.     drpatch(S, ntcurves, nucurves, ntsegs, nusegs, ntiter, nuiter);
  266.  
  267.     if (sync) {
  268.         vdevice.sync = 1;
  269.         (*vdevice.dev.Vsync)();
  270.     }
  271.  
  272.     /*
  273.      * Set the current (untransformed) world spot....
  274.      */
  275.     vdevice.cpW[V_X] = xlast;
  276.     vdevice.cpW[V_Y] = ylast;
  277.     vdevice.cpW[V_Z] = zlast;
  278. }
  279.  
  280. /*
  281.  * transformtensor
  282.  *
  283.  *    Transform the tensor S by the matrix m
  284.  */
  285. void
  286. transformtensor(S, m)
  287.     Tensor    S;
  288.     Matrix    m;
  289. {
  290.     Matrix    tmp, tmp2;
  291.     register    int    i;
  292.  
  293.     for (i = 0; i < 4; i++) {
  294.         extract(tmp, S, i);
  295.         mult4x4(tmp2, tmp, m);
  296.         replace(S, tmp2, i);
  297.     }
  298. }
  299.  
  300. /*
  301.  * replace
  302.  *
  303.  *    Does the reverse of extract.
  304.  */
  305. static    void
  306. replace(a, b, k)
  307.     Tensor    a;
  308.     Matrix    b;
  309.     int    k;
  310. {
  311.     int    i, j;
  312.  
  313.     /*
  314.      * Not unwound because it only gets called once per patch.
  315.      */
  316.     for (i = 0; i < 4; i++)
  317.         for (j = 0; j < 4; j++)
  318.             a[j][i][k] = b[i][j];
  319. }
  320.  
  321. /*
  322.  * drpatch
  323.  *
  324.  *    Actually does the work of drawing a patch.
  325.  */
  326. void
  327. drpatch(R, ntcurves, nucurves, ntsegs, nusegs, ntiter, nuiter)
  328.     Tensor    R;
  329.     int    ntcurves, nucurves, ntsegs, nusegs, ntiter, nuiter;
  330. {
  331.     Tensor    S;
  332.     Matrix    tmp;
  333.     int    i;
  334.  
  335.     /*
  336.          *  Copy R transposed into S
  337.          */
  338.     copytensortrans(S, R);
  339.  
  340.     for (i = 0; i < ntcurves; i++) {
  341.         extract(tmp, R, 0);
  342.         drcurve(ntsegs, tmp);
  343.         iterate(R, nuiter);
  344.     }
  345.  
  346.     /*
  347.      * Now using S...
  348.      */
  349.     for (i = 0; i < nucurves; i++) {
  350.         extract(tmp, S, 0);
  351.         drcurve(nusegs, tmp);
  352.         iterate(S, ntiter);
  353.     }
  354. }
  355.  
  356. /*
  357.  * iterate
  358.  *
  359.  *    Iterates the forward difference tensor R
  360.  */
  361. static    void
  362. iterate(R, n)
  363.     register Tensor    R;
  364.     int    n;
  365. {
  366.     register int    it;
  367.  
  368. /*
  369.  * Anyone for an unwound loop or two???
  370.  *
  371.  *    for (it = 0; it < n; it++) {
  372.  *        for (i = 0; i < 4; i++) 
  373.  *            for (j = 0; j < 4; j++)
  374.  *                for (k = 0; k < 3; k++)
  375.  *                    R[i][j][k] += R[i][j][k+1];
  376.  *    }
  377.  */
  378.     for (it = 0; it < n; it++) {
  379.         R[0][0][0] += R[0][0][1];
  380.         R[0][0][1] += R[0][0][2];
  381.         R[0][0][2] += R[0][0][3];
  382.  
  383.         R[0][1][0] += R[0][1][1];
  384.         R[0][1][1] += R[0][1][2];
  385.         R[0][1][2] += R[0][1][3];
  386.  
  387.         R[0][2][0] += R[0][2][1];
  388.         R[0][2][1] += R[0][2][2];
  389.         R[0][2][2] += R[0][2][3];
  390.  
  391.         R[0][3][0] += R[0][3][1];
  392.         R[0][3][1] += R[0][3][2];
  393.         R[0][3][2] += R[0][3][3];
  394.  
  395.         R[1][0][0] += R[1][0][1];
  396.         R[1][0][1] += R[1][0][2];
  397.         R[1][0][2] += R[1][0][3];
  398.  
  399.         R[1][1][0] += R[1][1][1];
  400.         R[1][1][1] += R[1][1][2];
  401.         R[1][1][2] += R[1][1][3];
  402.  
  403.         R[1][2][0] += R[1][2][1];
  404.         R[1][2][1] += R[1][2][2];
  405.         R[1][2][2] += R[1][2][3];
  406.  
  407.         R[1][3][0] += R[1][3][1];
  408.         R[1][3][1] += R[1][3][2];
  409.         R[1][3][2] += R[1][3][3];
  410.  
  411.         R[2][0][0] += R[2][0][1];
  412.         R[2][0][1] += R[2][0][2];
  413.         R[2][0][2] += R[2][0][3];
  414.  
  415.         R[2][1][0] += R[2][1][1];
  416.         R[2][1][1] += R[2][1][2];
  417.         R[2][1][2] += R[2][1][3];
  418.  
  419.         R[2][2][0] += R[2][2][1];
  420.         R[2][2][1] += R[2][2][2];
  421.         R[2][2][2] += R[2][2][3];
  422.  
  423.         R[2][3][0] += R[2][3][1];
  424.         R[2][3][1] += R[2][3][2];
  425.         R[2][3][2] += R[2][3][3];
  426.  
  427.         R[3][0][0] += R[3][0][1];
  428.         R[3][0][1] += R[3][0][2];
  429.         R[3][0][2] += R[3][0][3];
  430.  
  431.         R[3][1][0] += R[3][1][1];
  432.         R[3][1][1] += R[3][1][2];
  433.         R[3][1][2] += R[3][1][3];
  434.  
  435.         R[3][2][0] += R[3][2][1];
  436.         R[3][2][1] += R[3][2][2];
  437.         R[3][2][2] += R[3][2][3];
  438.  
  439.         R[3][3][0] += R[3][3][1];
  440.         R[3][3][1] += R[3][3][2];
  441.         R[3][3][2] += R[3][3][3];
  442.     }
  443. }
  444.  
  445. /*
  446.  * Extract the k'th column of the tensor a into the matrix b.
  447.  */
  448. static    void
  449. extract(b, a, k)
  450.     register Tensor a;
  451.     register Matrix    b;
  452.     register int    k;
  453. {
  454.     b[0][0] = a[0][0][k];
  455.     b[0][1] = a[1][0][k];
  456.     b[0][2] = a[2][0][k];
  457.     b[0][3] = a[3][0][k];
  458.  
  459.     b[1][0] = a[0][1][k];
  460.     b[1][1] = a[1][1][k];
  461.     b[1][2] = a[2][1][k];
  462.     b[1][3] = a[3][1][k];
  463.  
  464.     b[2][0] = a[0][2][k];
  465.     b[2][1] = a[1][2][k];
  466.     b[2][2] = a[2][2][k];
  467.     b[2][3] = a[3][2][k];
  468.  
  469.     b[3][0] = a[0][3][k];
  470.     b[3][1] = a[1][3][k];
  471.     b[3][2] = a[2][3][k];
  472.     b[3][3] = a[3][3][k];
  473. }
  474.  
  475. static    float
  476. addemup(m)
  477.     Matrix    m;
  478. {
  479.     return (m[0][0] + m[0][1] + m[0][2] + m[0][3] +
  480.         m[1][0] + m[1][1] + m[1][2] + m[1][3] +
  481.         m[2][0] + m[2][1] + m[2][2] + m[2][3] +
  482.         m[3][0] + m[3][1] + m[3][2] + m[3][3]);
  483. }
  484.