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