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