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 (at) cip.informatik.uni-erlangen.de) and 30 * Volker Weiss (vrweiss (at) 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
