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