Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | RSS feed

  1.  
  2. /*
  3.  * Mesa 3-D graphics library
  4.  *
  5.  * Copyright (C) 1999-2001  Brian Paul   All Rights Reserved.
  6.  *
  7.  * Permission is hereby granted, free of charge, to any person obtaining a
  8.  * copy of this software and associated documentation files (the "Software"),
  9.  * to deal in the Software without restriction, including without limitation
  10.  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  11.  * and/or sell copies of the Software, and to permit persons to whom the
  12.  * Software is furnished to do so, subject to the following conditions:
  13.  *
  14.  * The above copyright notice and this permission notice shall be included
  15.  * in all copies or substantial portions of the Software.
  16.  *
  17.  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  18.  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19.  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  20.  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  21.  * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  22.  * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  23.  * OTHER DEALINGS IN THE SOFTWARE.
  24.  */
  25.  
  26.  
  27. /*
  28.  * eval.c was written by
  29.  * Bernd Barsuhn (bdbarsuh@cip.informatik.uni-erlangen.de) and
  30.  * Volker Weiss (vrweiss@cip.informatik.uni-erlangen.de).
  31.  *
  32.  * My original implementation of evaluators was simplistic and didn't
  33.  * compute surface normal vectors properly.  Bernd and Volker applied
  34.  * used more sophisticated methods to get better results.
  35.  *
  36.  * Thanks guys!
  37.  */
  38.  
  39.  
  40. #include "main/glheader.h"
  41. #include "main/config.h"
  42. #include "m_eval.h"
  43.  
  44. static GLfloat inv_tab[MAX_EVAL_ORDER];
  45.  
  46.  
  47.  
  48. /*
  49.  * Horner scheme for Bezier curves
  50.  *
  51.  * Bezier curves can be computed via a Horner scheme.
  52.  * Horner is numerically less stable than the de Casteljau
  53.  * algorithm, but it is faster. For curves of degree n
  54.  * the complexity of Horner is O(n) and de Casteljau is O(n^2).
  55.  * Since stability is not important for displaying curve
  56.  * points I decided to use the Horner scheme.
  57.  *
  58.  * A cubic Bezier curve with control points b0, b1, b2, b3 can be
  59.  * written as
  60.  *
  61.  *        (([3]        [3]     )     [3]       )     [3]
  62.  * c(t) = (([0]*s*b0 + [1]*t*b1)*s + [2]*t^2*b2)*s + [3]*t^2*b3
  63.  *
  64.  *                                           [n]
  65.  * where s=1-t and the binomial coefficients [i]. These can
  66.  * be computed iteratively using the identity:
  67.  *
  68.  * [n]               [n  ]             [n]
  69.  * [i] = (n-i+1)/i * [i-1]     and     [0] = 1
  70.  */
  71.  
  72.  
  73. void
  74. _math_horner_bezier_curve(const GLfloat * cp, GLfloat * out, GLfloat t,
  75.                           GLuint dim, GLuint order)
  76. {
  77.    GLfloat s, powert, bincoeff;
  78.    GLuint i, k;
  79.  
  80.    if (order >= 2) {
  81.       bincoeff = (GLfloat) (order - 1);
  82.       s = 1.0F - t;
  83.  
  84.       for (k = 0; k < dim; k++)
  85.          out[k] = s * cp[k] + bincoeff * t * cp[dim + k];
  86.  
  87.       for (i = 2, cp += 2 * dim, powert = t * t; i < order;
  88.            i++, powert *= t, cp += dim) {
  89.          bincoeff *= (GLfloat) (order - i);
  90.          bincoeff *= inv_tab[i];
  91.  
  92.          for (k = 0; k < dim; k++)
  93.             out[k] = s * out[k] + bincoeff * powert * cp[k];
  94.       }
  95.    }
  96.    else {                       /* order=1 -> constant curve */
  97.  
  98.       for (k = 0; k < dim; k++)
  99.          out[k] = cp[k];
  100.    }
  101. }
  102.  
  103. /*
  104.  * Tensor product Bezier surfaces
  105.  *
  106.  * Again the Horner scheme is used to compute a point on a
  107.  * TP Bezier surface. First a control polygon for a curve
  108.  * on the surface in one parameter direction is computed,
  109.  * then the point on the curve for the other parameter
  110.  * direction is evaluated.
  111.  *
  112.  * To store the curve control polygon additional storage
  113.  * for max(uorder,vorder) points is needed in the
  114.  * control net cn.
  115.  */
  116.  
  117. void
  118. _math_horner_bezier_surf(GLfloat * cn, GLfloat * out, GLfloat u, GLfloat v,
  119.                          GLuint dim, GLuint uorder, GLuint vorder)
  120. {
  121.    GLfloat *cp = cn + uorder * vorder * dim;
  122.    GLuint i, uinc = vorder * dim;
  123.  
  124.    if (vorder > uorder) {
  125.       if (uorder >= 2) {
  126.          GLfloat s, poweru, bincoeff;
  127.          GLuint j, k;
  128.  
  129.          /* Compute the control polygon for the surface-curve in u-direction */
  130.          for (j = 0; j < vorder; j++) {
  131.             GLfloat *ucp = &cn[j * dim];
  132.  
  133.             /* Each control point is the point for parameter u on a */
  134.             /* curve defined by the control polygons in u-direction */
  135.             bincoeff = (GLfloat) (uorder - 1);
  136.             s = 1.0F - u;
  137.  
  138.             for (k = 0; k < dim; k++)
  139.                cp[j * dim + k] = s * ucp[k] + bincoeff * u * ucp[uinc + k];
  140.  
  141.             for (i = 2, ucp += 2 * uinc, poweru = u * u; i < uorder;
  142.                  i++, poweru *= u, ucp += uinc) {
  143.                bincoeff *= (GLfloat) (uorder - i);
  144.                bincoeff *= inv_tab[i];
  145.  
  146.                for (k = 0; k < dim; k++)
  147.                   cp[j * dim + k] =
  148.                      s * cp[j * dim + k] + bincoeff * poweru * ucp[k];
  149.             }
  150.          }
  151.  
  152.          /* Evaluate curve point in v */
  153.          _math_horner_bezier_curve(cp, out, v, dim, vorder);
  154.       }
  155.       else                      /* uorder=1 -> cn defines a curve in v */
  156.          _math_horner_bezier_curve(cn, out, v, dim, vorder);
  157.    }
  158.    else {                       /* vorder <= uorder */
  159.  
  160.       if (vorder > 1) {
  161.          GLuint i;
  162.  
  163.          /* Compute the control polygon for the surface-curve in u-direction */
  164.          for (i = 0; i < uorder; i++, cn += uinc) {
  165.             /* For constant i all cn[i][j] (j=0..vorder) are located */
  166.             /* on consecutive memory locations, so we can use        */
  167.             /* horner_bezier_curve to compute the control points     */
  168.  
  169.             _math_horner_bezier_curve(cn, &cp[i * dim], v, dim, vorder);
  170.          }
  171.  
  172.          /* Evaluate curve point in u */
  173.          _math_horner_bezier_curve(cp, out, u, dim, uorder);
  174.       }
  175.       else                      /* vorder=1 -> cn defines a curve in u */
  176.          _math_horner_bezier_curve(cn, out, u, dim, uorder);
  177.    }
  178. }
  179.  
  180. /*
  181.  * The direct de Casteljau algorithm is used when a point on the
  182.  * surface and the tangent directions spanning the tangent plane
  183.  * should be computed (this is needed to compute normals to the
  184.  * surface). In this case the de Casteljau algorithm approach is
  185.  * nicer because a point and the partial derivatives can be computed
  186.  * at the same time. To get the correct tangent length du and dv
  187.  * must be multiplied with the (u2-u1)/uorder-1 and (v2-v1)/vorder-1.
  188.  * Since only the directions are needed, this scaling step is omitted.
  189.  *
  190.  * De Casteljau needs additional storage for uorder*vorder
  191.  * values in the control net cn.
  192.  */
  193.  
  194. void
  195. _math_de_casteljau_surf(GLfloat * cn, GLfloat * out, GLfloat * du,
  196.                         GLfloat * dv, GLfloat u, GLfloat v, GLuint dim,
  197.                         GLuint uorder, GLuint vorder)
  198. {
  199.    GLfloat *dcn = cn + uorder * vorder * dim;
  200.    GLfloat us = 1.0F - u, vs = 1.0F - v;
  201.    GLuint h, i, j, k;
  202.    GLuint minorder = uorder < vorder ? uorder : vorder;
  203.    GLuint uinc = vorder * dim;
  204.    GLuint dcuinc = vorder;
  205.  
  206.    /* Each component is evaluated separately to save buffer space  */
  207.    /* This does not drasticaly decrease the performance of the     */
  208.    /* algorithm. If additional storage for (uorder-1)*(vorder-1)   */
  209.    /* points would be available, the components could be accessed  */
  210.    /* in the innermost loop which could lead to less cache misses. */
  211.  
  212. #define CN(I,J,K) cn[(I)*uinc+(J)*dim+(K)]
  213. #define DCN(I, J) dcn[(I)*dcuinc+(J)]
  214.    if (minorder < 3) {
  215.       if (uorder == vorder) {
  216.          for (k = 0; k < dim; k++) {
  217.             /* Derivative direction in u */
  218.             du[k] = vs * (CN(1, 0, k) - CN(0, 0, k)) +
  219.                v * (CN(1, 1, k) - CN(0, 1, k));
  220.  
  221.             /* Derivative direction in v */
  222.             dv[k] = us * (CN(0, 1, k) - CN(0, 0, k)) +
  223.                u * (CN(1, 1, k) - CN(1, 0, k));
  224.  
  225.             /* bilinear de Casteljau step */
  226.             out[k] = us * (vs * CN(0, 0, k) + v * CN(0, 1, k)) +
  227.                u * (vs * CN(1, 0, k) + v * CN(1, 1, k));
  228.          }
  229.       }
  230.       else if (minorder == uorder) {
  231.          for (k = 0; k < dim; k++) {
  232.             /* bilinear de Casteljau step */
  233.             DCN(1, 0) = CN(1, 0, k) - CN(0, 0, k);
  234.             DCN(0, 0) = us * CN(0, 0, k) + u * CN(1, 0, k);
  235.  
  236.             for (j = 0; j < vorder - 1; j++) {
  237.                /* for the derivative in u */
  238.                DCN(1, j + 1) = CN(1, j + 1, k) - CN(0, j + 1, k);
  239.                DCN(1, j) = vs * DCN(1, j) + v * DCN(1, j + 1);
  240.  
  241.                /* for the `point' */
  242.                DCN(0, j + 1) = us * CN(0, j + 1, k) + u * CN(1, j + 1, k);
  243.                DCN(0, j) = vs * DCN(0, j) + v * DCN(0, j + 1);
  244.             }
  245.  
  246.             /* remaining linear de Casteljau steps until the second last step */
  247.             for (h = minorder; h < vorder - 1; h++)
  248.                for (j = 0; j < vorder - h; j++) {
  249.                   /* for the derivative in u */
  250.                   DCN(1, j) = vs * DCN(1, j) + v * DCN(1, j + 1);
  251.  
  252.                   /* for the `point' */
  253.                   DCN(0, j) = vs * DCN(0, j) + v * DCN(0, j + 1);
  254.                }
  255.  
  256.             /* derivative direction in v */
  257.             dv[k] = DCN(0, 1) - DCN(0, 0);
  258.  
  259.             /* derivative direction in u */
  260.             du[k] = vs * DCN(1, 0) + v * DCN(1, 1);
  261.  
  262.             /* last linear de Casteljau step */
  263.             out[k] = vs * DCN(0, 0) + v * DCN(0, 1);
  264.          }
  265.       }
  266.       else {                    /* minorder == vorder */
  267.  
  268.          for (k = 0; k < dim; k++) {
  269.             /* bilinear de Casteljau step */
  270.             DCN(0, 1) = CN(0, 1, k) - CN(0, 0, k);
  271.             DCN(0, 0) = vs * CN(0, 0, k) + v * CN(0, 1, k);
  272.             for (i = 0; i < uorder - 1; i++) {
  273.                /* for the derivative in v */
  274.                DCN(i + 1, 1) = CN(i + 1, 1, k) - CN(i + 1, 0, k);
  275.                DCN(i, 1) = us * DCN(i, 1) + u * DCN(i + 1, 1);
  276.  
  277.                /* for the `point' */
  278.                DCN(i + 1, 0) = vs * CN(i + 1, 0, k) + v * CN(i + 1, 1, k);
  279.                DCN(i, 0) = us * DCN(i, 0) + u * DCN(i + 1, 0);
  280.             }
  281.  
  282.             /* remaining linear de Casteljau steps until the second last step */
  283.             for (h = minorder; h < uorder - 1; h++)
  284.                for (i = 0; i < uorder - h; i++) {
  285.                   /* for the derivative in v */
  286.                   DCN(i, 1) = us * DCN(i, 1) + u * DCN(i + 1, 1);
  287.  
  288.                   /* for the `point' */
  289.                   DCN(i, 0) = us * DCN(i, 0) + u * DCN(i + 1, 0);
  290.                }
  291.  
  292.             /* derivative direction in u */
  293.             du[k] = DCN(1, 0) - DCN(0, 0);
  294.  
  295.             /* derivative direction in v */
  296.             dv[k] = us * DCN(0, 1) + u * DCN(1, 1);
  297.  
  298.             /* last linear de Casteljau step */
  299.             out[k] = us * DCN(0, 0) + u * DCN(1, 0);
  300.          }
  301.       }
  302.    }
  303.    else if (uorder == vorder) {
  304.       for (k = 0; k < dim; k++) {
  305.          /* first bilinear de Casteljau step */
  306.          for (i = 0; i < uorder - 1; i++) {
  307.             DCN(i, 0) = us * CN(i, 0, k) + u * CN(i + 1, 0, k);
  308.             for (j = 0; j < vorder - 1; j++) {
  309.                DCN(i, j + 1) = us * CN(i, j + 1, k) + u * CN(i + 1, j + 1, k);
  310.                DCN(i, j) = vs * DCN(i, j) + v * DCN(i, j + 1);
  311.             }
  312.          }
  313.  
  314.          /* remaining bilinear de Casteljau steps until the second last step */
  315.          for (h = 2; h < minorder - 1; h++)
  316.             for (i = 0; i < uorder - h; i++) {
  317.                DCN(i, 0) = us * DCN(i, 0) + u * DCN(i + 1, 0);
  318.                for (j = 0; j < vorder - h; j++) {
  319.                   DCN(i, j + 1) = us * DCN(i, j + 1) + u * DCN(i + 1, j + 1);
  320.                   DCN(i, j) = vs * DCN(i, j) + v * DCN(i, j + 1);
  321.                }
  322.             }
  323.  
  324.          /* derivative direction in u */
  325.          du[k] = vs * (DCN(1, 0) - DCN(0, 0)) + v * (DCN(1, 1) - DCN(0, 1));
  326.  
  327.          /* derivative direction in v */
  328.          dv[k] = us * (DCN(0, 1) - DCN(0, 0)) + u * (DCN(1, 1) - DCN(1, 0));
  329.  
  330.          /* last bilinear de Casteljau step */
  331.          out[k] = us * (vs * DCN(0, 0) + v * DCN(0, 1)) +
  332.             u * (vs * DCN(1, 0) + v * DCN(1, 1));
  333.       }
  334.    }
  335.    else if (minorder == uorder) {
  336.       for (k = 0; k < dim; k++) {
  337.          /* first bilinear de Casteljau step */
  338.          for (i = 0; i < uorder - 1; i++) {
  339.             DCN(i, 0) = us * CN(i, 0, k) + u * CN(i + 1, 0, k);
  340.             for (j = 0; j < vorder - 1; j++) {
  341.                DCN(i, j + 1) = us * CN(i, j + 1, k) + u * CN(i + 1, j + 1, k);
  342.                DCN(i, j) = vs * DCN(i, j) + v * DCN(i, j + 1);
  343.             }
  344.          }
  345.  
  346.          /* remaining bilinear de Casteljau steps until the second last step */
  347.          for (h = 2; h < minorder - 1; h++)
  348.             for (i = 0; i < uorder - h; i++) {
  349.                DCN(i, 0) = us * DCN(i, 0) + u * DCN(i + 1, 0);
  350.                for (j = 0; j < vorder - h; j++) {
  351.                   DCN(i, j + 1) = us * DCN(i, j + 1) + u * DCN(i + 1, j + 1);
  352.                   DCN(i, j) = vs * DCN(i, j) + v * DCN(i, j + 1);
  353.                }
  354.             }
  355.  
  356.          /* last bilinear de Casteljau step */
  357.          DCN(2, 0) = DCN(1, 0) - DCN(0, 0);
  358.          DCN(0, 0) = us * DCN(0, 0) + u * DCN(1, 0);
  359.          for (j = 0; j < vorder - 1; j++) {
  360.             /* for the derivative in u */
  361.             DCN(2, j + 1) = DCN(1, j + 1) - DCN(0, j + 1);
  362.             DCN(2, j) = vs * DCN(2, j) + v * DCN(2, j + 1);
  363.  
  364.             /* for the `point' */
  365.             DCN(0, j + 1) = us * DCN(0, j + 1) + u * DCN(1, j + 1);
  366.             DCN(0, j) = vs * DCN(0, j) + v * DCN(0, j + 1);
  367.          }
  368.  
  369.          /* remaining linear de Casteljau steps until the second last step */
  370.          for (h = minorder; h < vorder - 1; h++)
  371.             for (j = 0; j < vorder - h; j++) {
  372.                /* for the derivative in u */
  373.                DCN(2, j) = vs * DCN(2, j) + v * DCN(2, j + 1);
  374.  
  375.                /* for the `point' */
  376.                DCN(0, j) = vs * DCN(0, j) + v * DCN(0, j + 1);
  377.             }
  378.  
  379.          /* derivative direction in v */
  380.          dv[k] = DCN(0, 1) - DCN(0, 0);
  381.  
  382.          /* derivative direction in u */
  383.          du[k] = vs * DCN(2, 0) + v * DCN(2, 1);
  384.  
  385.          /* last linear de Casteljau step */
  386.          out[k] = vs * DCN(0, 0) + v * DCN(0, 1);
  387.       }
  388.    }
  389.    else {                       /* minorder == vorder */
  390.  
  391.       for (k = 0; k < dim; k++) {
  392.          /* first bilinear de Casteljau step */
  393.          for (i = 0; i < uorder - 1; i++) {
  394.             DCN(i, 0) = us * CN(i, 0, k) + u * CN(i + 1, 0, k);
  395.             for (j = 0; j < vorder - 1; j++) {
  396.                DCN(i, j + 1) = us * CN(i, j + 1, k) + u * CN(i + 1, j + 1, k);
  397.                DCN(i, j) = vs * DCN(i, j) + v * DCN(i, j + 1);
  398.             }
  399.          }
  400.  
  401.          /* remaining bilinear de Casteljau steps until the second last step */
  402.          for (h = 2; h < minorder - 1; h++)
  403.             for (i = 0; i < uorder - h; i++) {
  404.                DCN(i, 0) = us * DCN(i, 0) + u * DCN(i + 1, 0);
  405.                for (j = 0; j < vorder - h; j++) {
  406.                   DCN(i, j + 1) = us * DCN(i, j + 1) + u * DCN(i + 1, j + 1);
  407.                   DCN(i, j) = vs * DCN(i, j) + v * DCN(i, j + 1);
  408.                }
  409.             }
  410.  
  411.          /* last bilinear de Casteljau step */
  412.          DCN(0, 2) = DCN(0, 1) - DCN(0, 0);
  413.          DCN(0, 0) = vs * DCN(0, 0) + v * DCN(0, 1);
  414.          for (i = 0; i < uorder - 1; i++) {
  415.             /* for the derivative in v */
  416.             DCN(i + 1, 2) = DCN(i + 1, 1) - DCN(i + 1, 0);
  417.             DCN(i, 2) = us * DCN(i, 2) + u * DCN(i + 1, 2);
  418.  
  419.             /* for the `point' */
  420.             DCN(i + 1, 0) = vs * DCN(i + 1, 0) + v * DCN(i + 1, 1);
  421.             DCN(i, 0) = us * DCN(i, 0) + u * DCN(i + 1, 0);
  422.          }
  423.  
  424.          /* remaining linear de Casteljau steps until the second last step */
  425.          for (h = minorder; h < uorder - 1; h++)
  426.             for (i = 0; i < uorder - h; i++) {
  427.                /* for the derivative in v */
  428.                DCN(i, 2) = us * DCN(i, 2) + u * DCN(i + 1, 2);
  429.  
  430.                /* for the `point' */
  431.                DCN(i, 0) = us * DCN(i, 0) + u * DCN(i + 1, 0);
  432.             }
  433.  
  434.          /* derivative direction in u */
  435.          du[k] = DCN(1, 0) - DCN(0, 0);
  436.  
  437.          /* derivative direction in v */
  438.          dv[k] = us * DCN(0, 2) + u * DCN(1, 2);
  439.  
  440.          /* last linear de Casteljau step */
  441.          out[k] = us * DCN(0, 0) + u * DCN(1, 0);
  442.       }
  443.    }
  444. #undef DCN
  445. #undef CN
  446. }
  447.  
  448.  
  449. /*
  450.  * Do one-time initialization for evaluators.
  451.  */
  452. void
  453. _math_init_eval(void)
  454. {
  455.    GLuint i;
  456.  
  457.    /* KW: precompute 1/x for useful x.
  458.     */
  459.    for (i = 1; i < MAX_EVAL_ORDER; i++)
  460.       inv_tab[i] = 1.0F / i;
  461. }
  462.