a4ae0395cad71146ced1f5bbd9ec9faa119e7481
1 /* $Id: m_eval.c,v 1.1 2000/12/26 05:09:31 keithw Exp $ */
4 * Mesa 3-D graphics library
7 * Copyright (C) 1999-2000 Brian Paul All Rights Reserved.
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:
16 * The above copyright notice and this permission notice shall be included
17 * in all copies or substantial portions of the Software.
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.
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).
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.
45 static GLfloat inv_tab
[MAX_EVAL_ORDER
];
50 * Horner scheme for Bezier curves
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.
59 * A cubic Bezier curve with control points b0, b1, b2, b3 can be
62 * (([3] [3] ) [3] ) [3]
63 * c(t) = (([0]*s*b0 + [1]*t*b1)*s + [2]*t^2*b2)*s + [3]*t^2*b3
66 * where s=1-t and the binomial coefficients [i]. These can
67 * be computed iteratively using the identity:
70 * [i] = (n-i+1)/i * [i-1] and [0] = 1
75 _math_horner_bezier_curve(const GLfloat
*cp
, GLfloat
*out
, GLfloat t
,
76 GLuint dim
, GLuint order
)
79 GLuint i
, k
, bincoeff
;
87 out
[k
] = s
*cp
[k
] + bincoeff
*t
*cp
[dim
+k
];
89 for(i
=2, cp
+=2*dim
, powert
=t
*t
; i
<order
; i
++, powert
*=t
, cp
+=dim
)
92 bincoeff
*= inv_tab
[i
];
95 out
[k
] = s
*out
[k
] + bincoeff
*powert
*cp
[k
];
98 else /* order=1 -> constant curve */
106 * Tensor product Bezier surfaces
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.
114 * To store the curve control polygon additional storage
115 * for max(uorder,vorder) points is needed in the
120 _math_horner_bezier_surf(GLfloat
*cn
, GLfloat
*out
, GLfloat u
, GLfloat v
,
121 GLuint dim
, GLuint uorder
, GLuint vorder
)
123 GLfloat
*cp
= cn
+ uorder
*vorder
*dim
;
124 GLuint i
, uinc
= vorder
*dim
;
131 GLuint j
, k
, bincoeff
;
133 /* Compute the control polygon for the surface-curve in u-direction */
134 for(j
=0; j
<vorder
; j
++)
136 GLfloat
*ucp
= &cn
[j
*dim
];
138 /* Each control point is the point for parameter u on a */
139 /* curve defined by the control polygons in u-direction */
144 cp
[j
*dim
+k
] = s
*ucp
[k
] + bincoeff
*u
*ucp
[uinc
+k
];
146 for(i
=2, ucp
+=2*uinc
, poweru
=u
*u
; i
<uorder
;
147 i
++, poweru
*=u
, ucp
+=uinc
)
149 bincoeff
*= uorder
-i
;
150 bincoeff
*= inv_tab
[i
];
153 cp
[j
*dim
+k
] = s
*cp
[j
*dim
+k
] + bincoeff
*poweru
*ucp
[k
];
157 /* Evaluate curve point in v */
158 _math_horner_bezier_curve(cp
, out
, v
, dim
, vorder
);
160 else /* uorder=1 -> cn defines a curve in v */
161 _math_horner_bezier_curve(cn
, out
, v
, dim
, vorder
);
163 else /* vorder <= uorder */
169 /* Compute the control polygon for the surface-curve in u-direction */
170 for(i
=0; i
<uorder
; i
++, cn
+= uinc
)
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 */
176 _math_horner_bezier_curve(cn
, &cp
[i
*dim
], v
, dim
, vorder
);
179 /* Evaluate curve point in u */
180 _math_horner_bezier_curve(cp
, out
, u
, dim
, uorder
);
182 else /* vorder=1 -> cn defines a curve in u */
183 _math_horner_bezier_curve(cn
, out
, u
, dim
, uorder
);
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.
197 * De Casteljau needs additional storage for uorder*vorder
198 * values in the control net cn.
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
)
206 GLfloat
*dcn
= cn
+ uorder
*vorder
*dim
;
207 GLfloat us
= 1.0-u
, vs
= 1.0-v
;
209 GLuint minorder
= uorder
< vorder
? uorder
: vorder
;
210 GLuint uinc
= vorder
*dim
;
211 GLuint dcuinc
= vorder
;
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. */
219 #define CN(I,J,K) cn[(I)*uinc+(J)*dim+(K)]
220 #define DCN(I, J) dcn[(I)*dcuinc+(J)]
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
));
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
));
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
));
240 else if(minorder
== uorder
)
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
);
248 for(j
=0; j
<vorder
-1; j
++)
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);
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);
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
++)
263 /* for the derivative in u */
264 DCN(1,j
) = vs
*DCN(1,j
) + v
*DCN(1,j
+1);
266 /* for the `point' */
267 DCN(0,j
) = vs
*DCN(0,j
) + v
*DCN(0,j
+1);
270 /* derivative direction in v */
271 dv
[k
] = DCN(0,1) - DCN(0,0);
273 /* derivative direction in u */
274 du
[k
] = vs
*DCN(1,0) + v
*DCN(1,1);
276 /* last linear de Casteljau step */
277 out
[k
] = vs
*DCN(0,0) + v
*DCN(0,1);
280 else /* minorder == vorder */
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
++)
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);
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);
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
++)
302 /* for the derivative in v */
303 DCN(i
,1) = us
*DCN(i
,1) + u
*DCN(i
+1,1);
305 /* for the `point' */
306 DCN(i
,0) = us
*DCN(i
,0) + u
*DCN(i
+1,0);
309 /* derivative direction in u */
310 du
[k
] = DCN(1,0) - DCN(0,0);
312 /* derivative direction in v */
313 dv
[k
] = us
*DCN(0,1) + u
*DCN(1,1);
315 /* last linear de Casteljau step */
316 out
[k
] = us
*DCN(0,0) + u
*DCN(1,0);
320 else if(uorder
== vorder
)
324 /* first bilinear de Casteljau step */
325 for(i
=0; i
<uorder
-1; i
++)
327 DCN(i
,0) = us
*CN(i
,0,k
) + u
*CN(i
+1,0,k
);
328 for(j
=0; j
<vorder
-1; j
++)
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);
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
++)
339 DCN(i
,0) = us
*DCN(i
,0) + u
*DCN(i
+1,0);
340 for(j
=0; j
<vorder
-h
; j
++)
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);
347 /* derivative direction in u */
348 du
[k
] = vs
*(DCN(1,0) - DCN(0,0)) +
349 v
*(DCN(1,1) - DCN(0,1));
351 /* derivative direction in v */
352 dv
[k
] = us
*(DCN(0,1) - DCN(0,0)) +
353 u
*(DCN(1,1) - DCN(1,0));
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));
360 else if(minorder
== uorder
)
364 /* first bilinear de Casteljau step */
365 for(i
=0; i
<uorder
-1; i
++)
367 DCN(i
,0) = us
*CN(i
,0,k
) + u
*CN(i
+1,0,k
);
368 for(j
=0; j
<vorder
-1; j
++)
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);
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
++)
379 DCN(i
,0) = us
*DCN(i
,0) + u
*DCN(i
+1,0);
380 for(j
=0; j
<vorder
-h
; j
++)
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);
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
++)
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);
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);
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
++)
405 /* for the derivative in u */
406 DCN(2,j
) = vs
*DCN(2,j
) + v
*DCN(2,j
+1);
408 /* for the `point' */
409 DCN(0,j
) = vs
*DCN(0,j
) + v
*DCN(0,j
+1);
412 /* derivative direction in v */
413 dv
[k
] = DCN(0,1) - DCN(0,0);
415 /* derivative direction in u */
416 du
[k
] = vs
*DCN(2,0) + v
*DCN(2,1);
418 /* last linear de Casteljau step */
419 out
[k
] = vs
*DCN(0,0) + v
*DCN(0,1);
422 else /* minorder == vorder */
426 /* first bilinear de Casteljau step */
427 for(i
=0; i
<uorder
-1; i
++)
429 DCN(i
,0) = us
*CN(i
,0,k
) + u
*CN(i
+1,0,k
);
430 for(j
=0; j
<vorder
-1; j
++)
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);
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
++)
441 DCN(i
,0) = us
*DCN(i
,0) + u
*DCN(i
+1,0);
442 for(j
=0; j
<vorder
-h
; j
++)
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);
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
++)
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);
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);
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
++)
467 /* for the derivative in v */
468 DCN(i
,2) = us
*DCN(i
,2) + u
*DCN(i
+1,2);
470 /* for the `point' */
471 DCN(i
,0) = us
*DCN(i
,0) + u
*DCN(i
+1,0);
474 /* derivative direction in u */
475 du
[k
] = DCN(1,0) - DCN(0,0);
477 /* derivative direction in v */
478 dv
[k
] = us
*DCN(0,2) + u
*DCN(1,2);
480 /* last linear de Casteljau step */
481 out
[k
] = us
*DCN(0,0) + u
*DCN(1,0);
490 * Do one-time initialization for evaluators.
492 void _math_init_eval( void )
496 /* KW: precompute 1/x for useful x.
498 for (i
= 1 ; i
< MAX_EVAL_ORDER
; i
++)
499 inv_tab
[i
] = 1.0 / i
;