1 //
2 // Copyright 2013 Pixar
3 //
4 // Licensed under the Apache License, Version 2.0 (the "Apache License")
5 // with the following modification; you may not use this file except in
6 // compliance with the Apache License and the following modification to it:
7 // Section 6. Trademarks. is deleted and replaced with:
8 //
9 // 6. Trademarks. This License does not grant permission to use the trade
10 // names, trademarks, service marks, or product names of the Licensor
11 // and its affiliates, except as required to comply with Section 4(c) of
12 // the License and to reproduce the content of the NOTICE file.
13 //
14 // You may obtain a copy of the Apache License at
15 //
16 // http://www.apache.org/licenses/LICENSE-2.0
17 //
18 // Unless required by applicable law or agreed to in writing, software
19 // distributed under the Apache License with the above modification is
20 // distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
21 // KIND, either express or implied. See the Apache License for the specific
22 // language governing permissions and limitations under the Apache License.
23 //
24
25 #include "../far/patchBasis.h"
26 #include "../far/patchDescriptor.h"
27
28 #include <cassert>
29 #include <cstring>
30 #include <cmath>
31 #include <cstdio>
32
33 namespace OpenSubdiv {
34 namespace OPENSUBDIV_VERSION {
35
36 namespace Far {
37 namespace internal {
38
39
40 //
41 // Basis support for quadrilateral patches:
42 //
43 // Quadrilateral patches are parameterized in terms of (s,t) as follows:
44 //
45 // (1,0) *---------* (1,1)
46 // | 3 2 |
47 // t | |
48 // | |
49 // | 0 1 |
50 // (0,0) *---------* (1,0)
51 // s
52 //
53
54
55 //
56 // Simple bilinear quad:
57 //
58 template <typename REAL>
59 int
EvalBasisLinear(REAL s,REAL t,REAL wP[4],REAL wDs[4],REAL wDt[4],REAL wDss[4],REAL wDst[4],REAL wDtt[4])60 EvalBasisLinear(REAL s, REAL t,
61 REAL wP[4], REAL wDs[4], REAL wDt[4],
62 REAL wDss[4], REAL wDst[4], REAL wDtt[4]) {
63
64 REAL sC = 1.0f - s;
65 REAL tC = 1.0f - t;
66
67 if (wP) {
68 wP[0] = sC * tC;
69 wP[1] = s * tC;
70 wP[2] = s * t;
71 wP[3] = sC * t;
72 }
73 if (wDs && wDt) {
74 wDs[0] = -tC;
75 wDs[1] = tC;
76 wDs[2] = t;
77 wDs[3] = -t;
78
79 wDt[0] = -sC;
80 wDt[1] = -s;
81 wDt[2] = s;
82 wDt[3] = sC;
83
84 if (wDss && wDst && wDtt) {
85 for(int i = 0; i < 4; ++i) {
86 wDss[i] = 0.0f;
87 wDtt[i] = 0.0f;
88 }
89
90 wDst[0] = 1.0f;
91 wDst[1] = -1.0f;
92 wDst[2] = -1.0f;
93 wDst[3] = 1.0f;
94 }
95 }
96 return 4;
97 }
98
99 //
100 // Bicubic BSpline patch:
101 //
102 // 12-----13------14-----15
103 // | | | |
104 // | | | |
105 // 8------9------10-----11
106 // | | t | |
107 // | | s | |
108 // 4------5------6------7
109 // | | | |
110 // | | | |
111 // O------1------2------3
112 //
113 // The basis if a bicubic BSpline patch is a tensor product, which we make
114 // use of here by evaluating, differentiating and combining basis functions
115 // in each of the two parametric directions.
116 //
117 // Not all 16 points will be present. The boundary mask indicates boundary
118 // edges beyond which phantom points are implicitly extrapolated. Weights
119 // for missing points are set to zero while those contributing to their
120 // implicit extrapolation will be adjusted.
121 //
122 namespace {
123 //
124 // Cubic BSpline curve basis evaluation:
125 //
126 template <typename REAL>
127 void
evalBSplineCurve(REAL t,REAL wP[4],REAL wDP[4],REAL wDP2[4])128 evalBSplineCurve(REAL t, REAL wP[4], REAL wDP[4], REAL wDP2[4]) {
129
130 REAL const one6th = (REAL)(1.0 / 6.0);
131
132 REAL t2 = t * t;
133 REAL t3 = t * t2;
134
135 wP[0] = one6th * (1.0f - 3.0f*(t - t2) - t3);
136 wP[1] = one6th * (4.0f - 6.0f*t2 + 3.0f*t3);
137 wP[2] = one6th * (1.0f + 3.0f*(t + t2 - t3));
138 wP[3] = one6th * ( t3);
139
140 if (wDP) {
141 wDP[0] = -0.5f*t2 + t - 0.5f;
142 wDP[1] = 1.5f*t2 - 2.0f*t;
143 wDP[2] = -1.5f*t2 + t + 0.5f;
144 wDP[3] = 0.5f*t2;
145 }
146 if (wDP2) {
147 wDP2[0] = - t + 1.0f;
148 wDP2[1] = 3.0f * t - 2.0f;
149 wDP2[2] = -3.0f * t + 1.0f;
150 wDP2[3] = t;
151 }
152 }
153
154 //
155 // Weight adjustments to account for phantom end points:
156 //
157 template <typename REAL>
158 void
adjustBSplineBoundaryWeights(int boundary,REAL w[16])159 adjustBSplineBoundaryWeights(int boundary, REAL w[16]) {
160
161 if ((boundary & 1) != 0) {
162 for (int i = 0; i < 4; ++i) {
163 w[i + 8] -= w[i + 0];
164 w[i + 4] += w[i + 0] * 2.0f;
165 w[i + 0] = 0.0f;
166 }
167 }
168 if ((boundary & 2) != 0) {
169 for (int i = 0; i < 16; i += 4) {
170 w[i + 1] -= w[i + 3];
171 w[i + 2] += w[i + 3] * 2.0f;
172 w[i + 3] = 0.0f;
173 }
174 }
175 if ((boundary & 4) != 0) {
176 for (int i = 0; i < 4; ++i) {
177 w[i + 4] -= w[i + 12];
178 w[i + 8] += w[i + 12] * 2.0f;
179 w[i + 12] = 0.0f;
180 }
181 }
182 if ((boundary & 8) != 0) {
183 for (int i = 0; i < 16; i += 4) {
184 w[i + 2] -= w[i + 0];
185 w[i + 1] += w[i + 0] * 2.0f;
186 w[i + 0] = 0.0f;
187 }
188 }
189 }
190
191 template <typename REAL>
192 void
boundBasisBSpline(int boundary,REAL wP[16],REAL wDs[16],REAL wDt[16],REAL wDss[16],REAL wDst[16],REAL wDtt[16])193 boundBasisBSpline(int boundary,
194 REAL wP[16], REAL wDs[16], REAL wDt[16],
195 REAL wDss[16], REAL wDst[16], REAL wDtt[16]) {
196
197 if (wP) {
198 adjustBSplineBoundaryWeights(boundary, wP);
199 }
200 if (wDs && wDt) {
201 adjustBSplineBoundaryWeights(boundary, wDs);
202 adjustBSplineBoundaryWeights(boundary, wDt);
203
204 if (wDss && wDst && wDtt) {
205 adjustBSplineBoundaryWeights(boundary, wDss);
206 adjustBSplineBoundaryWeights(boundary, wDst);
207 adjustBSplineBoundaryWeights(boundary, wDtt);
208 }
209 }
210 }
211
212 } // end namespace
213
214 template <typename REAL>
215 int
EvalBasisBSpline(REAL s,REAL t,REAL wP[16],REAL wDs[16],REAL wDt[16],REAL wDss[16],REAL wDst[16],REAL wDtt[16])216 EvalBasisBSpline(REAL s, REAL t,
217 REAL wP[16], REAL wDs[16], REAL wDt[16],
218 REAL wDss[16], REAL wDst[16], REAL wDtt[16]) {
219
220 REAL sWeights[4], tWeights[4], dsWeights[4], dtWeights[4], dssWeights[4], dttWeights[4];
221
222 evalBSplineCurve(s, wP ? sWeights : 0, wDs ? dsWeights : 0, wDss ? dssWeights : 0);
223 evalBSplineCurve(t, wP ? tWeights : 0, wDt ? dtWeights : 0, wDtt ? dttWeights : 0);
224
225 if (wP) {
226 for (int i = 0; i < 4; ++i) {
227 for (int j = 0; j < 4; ++j) {
228 wP[4*i+j] = sWeights[j] * tWeights[i];
229 }
230 }
231 }
232 if (wDs && wDt) {
233 for (int i = 0; i < 4; ++i) {
234 for (int j = 0; j < 4; ++j) {
235 wDs[4*i+j] = dsWeights[j] * tWeights[i];
236 wDt[4*i+j] = sWeights[j] * dtWeights[i];
237 }
238 }
239
240 if (wDss && wDst && wDtt) {
241 for (int i = 0; i < 4; ++i) {
242 for (int j = 0; j < 4; ++j) {
243 wDss[4*i+j] = dssWeights[j] * tWeights[i];
244 wDst[4*i+j] = dsWeights[j] * dtWeights[i];
245 wDtt[4*i+j] = sWeights[j] * dttWeights[i];
246 }
247 }
248 }
249 }
250 return 16;
251 }
252
253 //
254 // Bicubic Bezier patch:
255 //
256 // 12-----13------14-----15
257 // | | | |
258 // | | | |
259 // 8------9------10-----11
260 // | | | |
261 // | | | |
262 // 4------5------6------7
263 // | t | | |
264 // | s | | |
265 // O------1------2------3
266 //
267 // As was the case with the BSpline patch, a bicubic Bezier patch can also
268 // make use of its tensor product property by evaluating, differentiating
269 // and combining basis functions in each of the two parametric directions.
270 //
271 namespace {
272 //
273 // Cubic Bezier curve basis evaluation:
274 //
275 template <typename REAL>
276 void
evalBezierCurve(REAL t,REAL wP[4],REAL wDP[4],REAL wDP2[4])277 evalBezierCurve(REAL t, REAL wP[4], REAL wDP[4], REAL wDP2[4]) {
278
279 // The four uniform cubic Bezier basis functions (in terms of t and its
280 // complement tC) evaluated at t:
281 REAL t2 = t*t;
282 REAL tC = 1.0f - t;
283 REAL tC2 = tC * tC;
284
285 wP[0] = tC2 * tC;
286 wP[1] = tC2 * t * 3.0f;
287 wP[2] = t2 * tC * 3.0f;
288 wP[3] = t2 * t;
289
290 // Derivatives of the above four basis functions at t:
291 if (wDP) {
292 wDP[0] = -3.0f * tC2;
293 wDP[1] = 9.0f * t2 - 12.0f * t + 3.0f;
294 wDP[2] = -9.0f * t2 + 6.0f * t;
295 wDP[3] = 3.0f * t2;
296 }
297
298 // Second derivatives of the basis functions at t:
299 if (wDP2) {
300 wDP2[0] = 6.0f * tC;
301 wDP2[1] = 18.0f * t - 12.0f;
302 wDP2[2] = -18.0f * t + 6.0f;
303 wDP2[3] = 6.0f * t;
304 }
305 }
306 } // end namespace
307
308 template <typename REAL>
309 int
EvalBasisBezier(REAL s,REAL t,REAL wP[16],REAL wDs[16],REAL wDt[16],REAL wDss[16],REAL wDst[16],REAL wDtt[16])310 EvalBasisBezier(REAL s, REAL t,
311 REAL wP[16], REAL wDs[16], REAL wDt[16],
312 REAL wDss[16], REAL wDst[16], REAL wDtt[16]) {
313
314 REAL sWeights[4], tWeights[4], dsWeights[4], dtWeights[4], dssWeights[4], dttWeights[4];
315
316 evalBezierCurve(s, wP ? sWeights : 0, wDs ? dsWeights : 0, wDss ? dssWeights : 0);
317 evalBezierCurve(t, wP ? tWeights : 0, wDt ? dtWeights : 0, wDtt ? dttWeights : 0);
318
319 if (wP) {
320 for (int i = 0; i < 4; ++i) {
321 for (int j = 0; j < 4; ++j) {
322 wP[4*i+j] = sWeights[j] * tWeights[i];
323 }
324 }
325 }
326 if (wDs && wDt) {
327 for (int i = 0; i < 4; ++i) {
328 for (int j = 0; j < 4; ++j) {
329 wDs[4*i+j] = dsWeights[j] * tWeights[i];
330 wDt[4*i+j] = sWeights[j] * dtWeights[i];
331 }
332 }
333
334 if (wDss && wDst && wDtt) {
335 for (int i = 0; i < 4; ++i) {
336 for (int j = 0; j < 4; ++j) {
337 wDss[4*i+j] = dssWeights[j] * tWeights[i];
338 wDst[4*i+j] = dsWeights[j] * dtWeights[i];
339 wDtt[4*i+j] = sWeights[j] * dttWeights[i];
340 }
341 }
342 }
343 }
344 return 16;
345 }
346
347 //
348 // Cubic Gregory patch:
349 //
350 // P3 e3- e2+ P2
351 // 15------17-------11--------10
352 // | | | |
353 // | | | |
354 // | | f3- | f2+ |
355 // | 19 13 |
356 // e3+ 16-----18 14-----12 e2-
357 // | f3+ f2- |
358 // | |
359 // | |
360 // | f0- f1+ |
361 // e0- 2------4 8------6 e1+
362 // | 3 9 |
363 // | | f0+ | f1- |
364 // | t | | |
365 // | s | | |
366 // O--------1--------7--------5
367 // P0 e0+ e1- P1
368 //
369 // The 20-point cubic Gregory patch is an extension of the 16-point bicubic
370 // Bezier patch with the 4 interior points of the Bezier patch replaced with
371 // pairs of points (face points -- fi+ and fi-) that are rationally combined.
372 //
373 // The point ordering of the Gregory patch deviates considerably from the
374 // BSpline and Bezier patches by grouping the 5 points at each corner and
375 // ordering the groups by corner index.
376 //
377 template <typename REAL>
378 int
EvalBasisGregory(REAL s,REAL t,REAL point[20],REAL wDs[20],REAL wDt[20],REAL wDss[20],REAL wDst[20],REAL wDtt[20])379 EvalBasisGregory(REAL s, REAL t,
380 REAL point[20], REAL wDs[20], REAL wDt[20],
381 REAL wDss[20], REAL wDst[20], REAL wDtt[20]) {
382
383 // Indices of boundary and interior points and their corresponding Bezier points
384 // (this can be reduced with more direct indexing and unrolling of loops):
385 //
386 static int const boundaryGregory[12] = { 0, 1, 7, 5, 2, 6, 16, 12, 15, 17, 11, 10 };
387 static int const boundaryBezSCol[12] = { 0, 1, 2, 3, 0, 3, 0, 3, 0, 1, 2, 3 };
388 static int const boundaryBezTRow[12] = { 0, 0, 0, 0, 1, 1, 2, 2, 3, 3, 3, 3 };
389
390 static int const interiorGregory[8] = { 3, 4, 8, 9, 13, 14, 18, 19 };
391 static int const interiorBezSCol[8] = { 1, 1, 2, 2, 2, 2, 1, 1 };
392 static int const interiorBezTRow[8] = { 1, 1, 1, 1, 2, 2, 2, 2 };
393
394 //
395 // Bezier basis functions are denoted with B while the rational multipliers for the
396 // interior points will be denoted G -- so we have B(s), B(t) and G(s,t):
397 //
398 // Directional Bezier basis functions B at s and t:
399 REAL Bs[4], Bds[4], Bdss[4];
400 REAL Bt[4], Bdt[4], Bdtt[4];
401
402 evalBezierCurve(s, Bs, wDs ? Bds : 0, wDss ? Bdss : 0);
403 evalBezierCurve(t, Bt, wDt ? Bdt : 0, wDtt ? Bdtt : 0);
404
405 // Rational multipliers G at s and t:
406 REAL sC = 1.0f - s;
407 REAL tC = 1.0f - t;
408
409 // Use <= here to avoid compiler warnings -- the sums should always be non-negative:
410 REAL df0 = s + t; df0 = (df0 <= 0.0f) ? (REAL)1.0f : (1.0f / df0);
411 REAL df1 = sC + t; df1 = (df1 <= 0.0f) ? (REAL)1.0f : (1.0f / df1);
412 REAL df2 = sC + tC; df2 = (df2 <= 0.0f) ? (REAL)1.0f : (1.0f / df2);
413 REAL df3 = s + tC; df3 = (df3 <= 0.0f) ? (REAL)1.0f : (1.0f / df3);
414
415 REAL G[8] = { s*df0, t*df0, t*df1, sC*df1, sC*df2, tC*df2, tC*df3, s*df3 };
416
417 // Combined weights for boundary and interior points:
418 for (int i = 0; i < 12; ++i) {
419 point[boundaryGregory[i]] = Bs[boundaryBezSCol[i]] * Bt[boundaryBezTRow[i]];
420 }
421 for (int i = 0; i < 8; ++i) {
422 point[interiorGregory[i]] = Bs[interiorBezSCol[i]] * Bt[interiorBezTRow[i]] * G[i];
423 }
424
425 //
426 // For derivatives, the basis functions for the interior points are rational and ideally
427 // require appropriate differentiation, i.e. product rule for the combination of B and G
428 // and the quotient rule for the rational G itself. As initially proposed by Loop et al
429 // though, the approximation using the 16 Bezier points arising from the G(s,t) has
430 // proved adequate (and is what the GPU shaders use) so we continue to use that here.
431 //
432 // An implementation of the true derivatives is provided and conditionally compiled for
433 // those that require it, e.g.:
434 //
435 // dclyde's note: skipping half of the product rule like this does seem to change the
436 // result a lot in my tests. This is not a runtime bottleneck for cloth sims anyway
437 // so I'm just using the accurate version.
438 //
439 if (wDs && wDt) {
440 bool find_second_partials = wDs && wDst && wDtt;
441
442 // Combined weights for boundary points -- simple tensor products:
443 for (int i = 0; i < 12; ++i) {
444 int iDst = boundaryGregory[i];
445 int tRow = boundaryBezTRow[i];
446 int sCol = boundaryBezSCol[i];
447
448 wDs[iDst] = Bds[sCol] * Bt[tRow];
449 wDt[iDst] = Bdt[tRow] * Bs[sCol];
450
451 if (find_second_partials) {
452 wDss[iDst] = Bdss[sCol] * Bt[tRow];
453 wDst[iDst] = Bds[sCol] * Bdt[tRow];
454 wDtt[iDst] = Bs[sCol] * Bdtt[tRow];
455 }
456 }
457
458 #ifndef OPENSUBDIV_GREGORY_EVAL_TRUE_DERIVATIVES
459 // Approximation to the true Gregory derivatives by differentiating the Bezier patch
460 // unique to the given (s,t), i.e. having F = (g^+ * f^+) + (g^- * f^-) as its four
461 // interior points:
462 //
463 // Combined weights for interior points -- tensor products with G+ or G-:
464 for (int i = 0; i < 8; ++i) {
465 int iDst = interiorGregory[i];
466 int tRow = interiorBezTRow[i];
467 int sCol = interiorBezSCol[i];
468
469 wDs[iDst] = Bds[sCol] * Bt[tRow] * G[i];
470 wDt[iDst] = Bdt[tRow] * Bs[sCol] * G[i];
471
472 if (find_second_partials) {
473 wDss[iDst] = Bdss[sCol] * Bt[tRow] * G[i];
474 wDst[iDst] = Bds[sCol] * Bdt[tRow] * G[i];
475 wDtt[iDst] = Bs[sCol] * Bdtt[tRow] * G[i];
476 }
477 }
478 #else
479 // True Gregory derivatives using appropriate differentiation of composite functions:
480 //
481 // Note that for G(s,t) = N(s,t) / D(s,t), all N' and D' are trivial constants (which
482 // simplifies things for higher order derivatives). And while each pair of functions
483 // G (i.e. the G+ and G- corresponding to points f+ and f-) must sum to 1 to ensure
484 // Bezier equivalence (when f+ = f-), the pairs of G' must similarly sum to 0. So we
485 // can potentially compute only one of the pair and negate the result for the other
486 // (and with 4 or 8 computations involving these constants, this is all very SIMD
487 // friendly...) but for now we treat all 8 independently for simplicity.
488 //
489 //REAL N[8] = { s, t, t, sC, sC, tC, tC, s };
490 REAL D[8] = { df0, df0, df1, df1, df2, df2, df3, df3 };
491
492 static REAL const Nds[8] = { 1.0f, 0.0f, 0.0f, -1.0f, -1.0f, 0.0f, 0.0f, 1.0f };
493 static REAL const Ndt[8] = { 0.0f, 1.0f, 1.0f, 0.0f, 0.0f, -1.0f, -1.0f, 0.0f };
494
495 static REAL const Dds[8] = { 1.0f, 1.0f, -1.0f, -1.0f, -1.0f, -1.0f, 1.0f, 1.0f };
496 static REAL const Ddt[8] = { 1.0f, 1.0f, 1.0f, 1.0f, -1.0f, -1.0f, -1.0f, -1.0f };
497
498 // Combined weights for interior points -- combinations of B, B', G and G':
499 for (int i = 0; i < 8; ++i) {
500 int iDst = interiorGregory[i];
501 int tRow = interiorBezTRow[i];
502 int sCol = interiorBezSCol[i];
503
504 // Quotient rule for G' (re-expressed in terms of G to simplify (and D = 1/D)):
505 REAL Gds = (Nds[i] - Dds[i] * G[i]) * D[i];
506 REAL Gdt = (Ndt[i] - Ddt[i] * G[i]) * D[i];
507
508 // Product rule combining B and B' with G and G':
509 wDs[iDst] = (Bds[sCol] * G[i] + Bs[sCol] * Gds) * Bt[tRow];
510 wDt[iDst] = (Bdt[tRow] * G[i] + Bt[tRow] * Gdt) * Bs[sCol];
511
512 if (find_second_partials) {
513 REAL Dsqr_inv = D[i]*D[i];
514
515 REAL Gdss = 2.0f * Dds[i] * Dsqr_inv * (G[i] * Dds[i] - Nds[i]);
516 REAL Gdst = Dsqr_inv * (2.0f * G[i] * Dds[i] * Ddt[i] - Nds[i] * Ddt[i] - Ndt[i] * Dds[i]);
517 REAL Gdtt = 2.0f * Ddt[i] * Dsqr_inv * (G[i] * Ddt[i] - Ndt[i]);
518
519 wDss[iDst] = (Bdss[sCol] * G[i] + 2.0f * Bds[sCol] * Gds + Bs[sCol] * Gdss) * Bt[tRow];
520 wDst[iDst] = Bt[tRow] * (Bs[sCol] * Gdst + Bds[sCol] * Gdt) +
521 Bdt[tRow] * (Bds[sCol] * G[i] + Bs[sCol] * Gds);
522 wDtt[iDst] = (Bdtt[tRow] * G[i] + 2.0f * Bdt[tRow] * Gdt + Bt[tRow] * Gdtt) * Bs[sCol];
523 }
524 }
525 #endif
526 }
527 return 20;
528 }
529
530
531 //
532 // Basis support for triangular patches:
533 //
534 // Triangular patches may be evaluated in barycentric (trivariate) or
535 // bivariate form, depending on the complexity of their basis functions.
536 // The parametric orientation for a triangle is as follows:
537 //
538 // (1,0)
539 // *
540 // . .
541 // t . 2 .
542 // . .
543 // . 0 1 .
544 // (0,0) *---------* (1,0)
545 // s
546 //
547 // With the origin (0,0) -- barycentric (0,0,w = 1) -- oriented at the
548 // corner V0, the corners V0, V1, and V2 correspond to barycentric
549 // coordinates W, U and V. This is consistent with GPU tessellation
550 // shaders, but not with many publications where the corners correspond
551 // more intuitively to U, V and W.
552 //
553
554
555 //
556 // Simple linear triangle:
557 //
558 template <typename REAL>
559 int
EvalBasisLinearTri(REAL s,REAL t,REAL wP[3],REAL wDs[3],REAL wDt[3],REAL wDss[3],REAL wDst[3],REAL wDtt[3])560 EvalBasisLinearTri(REAL s, REAL t,
561 REAL wP[3], REAL wDs[3], REAL wDt[3],
562 REAL wDss[3], REAL wDst[3], REAL wDtt[3]) {
563
564 if (wP) {
565 wP[0] = 1.0f - s - t;
566 wP[1] = s;
567 wP[2] = t;
568 }
569 if (wDs && wDt) {
570 wDs[0] = -1.0f;
571 wDs[1] = 1.0f;
572 wDs[2] = 0.0f;
573
574 wDt[0] = -1.0f;
575 wDt[1] = 0.0f;
576 wDt[2] = 1.0f;
577
578 if (wDss && wDst && wDtt) {
579 wDss[0] = wDss[1] = wDss[2] = 0.0f;
580 wDst[0] = wDst[1] = wDst[2] = 0.0f;
581 wDtt[0] = wDtt[1] = wDtt[2] = 0.0f;
582 }
583 }
584 return 3;
585 }
586
587
588 //
589 // Quartic Box spline triangle:
590 //
591 // Points for the quartic triangular Box spline (representing regular
592 // patches for Loop subdivision) are as follows:
593 //
594 // 10-----11
595 // . . . .
596 // . . . .
597 // 7-----8-----9
598 // . . . . . .
599 // . . . . . .
600 // 3-----4-----5-----6
601 // . . . . . .
602 // . . . . . /
603 // 0-----1-----2
604 //
605 // Stam provided the basis functions for these patches in terms of barycentric
606 // coordinates (u,v,w) (see Stam's "Evaluation of Loop Subdivision Surfaces").
607 // Unfortunately, unlike the basis functions for a quartic Bezier triangle,
608 // they are not very compact -- 3 functions involving 9 quartic terms and 3
609 // others involving 15 quartic terms. (In contrast, the maximum number of
610 // terms in bivariate form is 15.)
611 //
612 // Since we also need to differentiate with respect to u and v, we eliminate w
613 // and use the coefficient matrix C multiplied by the set of monomials M
614 // evaluated at (u,v), i.e. the full set of basis functions is:
615 //
616 // B(u,v) = C * M(u,v)
617 //
618 // where
619 //
620 // M(u,v) = { 1, u,v, uu,uv,vv, uuu,uuv,uvv,vvv, uuuu,uuuv,uuvv,uvvv,vvvv }
621 //
622 // and the 12 x 15 matrix C is as follows, scaled by a common factor of 1/12:
623 //
624 // { 1, -2,-4, 0, 6, 6, 2, 0, -6, -4, -1, -2, 0, 2, 1 },
625 // { 1, 2,-2, 0, -6, 0, -4, 0, 6, 2, 2, 4, 0, -2, -1 },
626 // { 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, -1, -2, 0, 0, 0 },
627 // { 1, -4,-2, 6, 6, 0, -4, -6, 0, 2, 1, 2, 0, -2, -1 },
628 // { 6, 0, 0, -12,-12,-12, 8, 12, 12, 8, -1, -2, 0, -2, -1 },
629 // { 1, 4, 2, 6, 6, 0, -4, -6,-12, -4, -1, -2, 0, 4, 2 },
630 // { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0 },
631 // { 1, -2, 2, 0, -6, 0, 2, 6, 0, -4, -1, -2, 0, 4, 2 },
632 // { 1, 2, 4, 0, 6, 6, -4,-12, -6, -4, 2, 4, 0, -2, -1 },
633 // { 0, 0, 0, 0, 0, 0, 2, 6, 6, 2, -1, -2, 0, -2, -1 },
634 // { 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, -2, -1 },
635 // { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1 }
636 //
637 // Differentiating the monomials and refactoring yields a unique set of
638 // coefficients for each of the derivatives, which we multiply by M(u,v).
639 //
640 namespace {
641 template <typename REAL>
642 inline void
evalBivariateMonomialsQuartic(REAL s,REAL t,REAL M[])643 evalBivariateMonomialsQuartic(REAL s, REAL t, REAL M[]) {
644
645 M[0] = 1.0;
646
647 M[1] = s;
648 M[2] = t;
649
650 M[3] = s * s;
651 M[4] = s * t;
652 M[5] = t * t;
653
654 M[6] = M[3] * s;
655 M[7] = M[4] * s;
656 M[8] = M[4] * t;
657 M[9] = M[5] * t;
658
659 M[10] = M[6] * s;
660 M[11] = M[7] * s;
661 M[12] = M[3] * M[5];
662 M[13] = M[8] * t;
663 M[14] = M[9] * t;
664 }
665
666 template <typename REAL>
667 void
evalBoxSplineTriDerivWeights(REAL const stMonomials[],int ds,int dt,REAL w[])668 evalBoxSplineTriDerivWeights(REAL const stMonomials[], int ds, int dt, REAL w[]) {
669
670 REAL const * M = stMonomials;
671
672 REAL S = 1.0;
673
674 int totalOrder = ds + dt;
675 if (totalOrder == 0) {
676 S *= (REAL) (1.0 / 12.0);
677
678 w[0] = S * (1 - 2*M[1] - 4*M[2] + 6*M[4] + 6*M[5] + 2*M[6] - 6*M[8] - 4*M[9] - M[10] - 2*M[11] + 2*M[13] + M[14]);
679 w[1] = S * (1 + 2*M[1] - 2*M[2] - 6*M[4] - 4*M[6] + 6*M[8] + 2*M[9] + 2*M[10] + 4*M[11] - 2*M[13] - M[14]);
680 w[2] = S * ( 2*M[6] - M[10] - 2*M[11] );
681 w[3] = S * (1 - 4*M[1] - 2*M[2] + 6*M[3] + 6*M[4] - 4*M[6] - 6*M[7] + 2*M[9] + M[10] + 2*M[11] - 2*M[13] - M[14]);
682 w[4] = S * (6 -12*M[3] -12*M[4] -12*M[5] + 8*M[6] +12*M[7] +12*M[8] + 8*M[9] - M[10] - 2*M[11] - 2*M[13] - M[14]);
683 w[5] = S * (1 + 4*M[1] + 2*M[2] + 6*M[3] + 6*M[4] - 4*M[6] - 6*M[7] -12*M[8] - 4*M[9] - M[10] - 2*M[11] + 4*M[13] + 2*M[14]);
684 w[6] = S * ( M[10] + 2*M[11] );
685 w[7] = S * (1 - 2*M[1] + 2*M[2] - 6*M[4] + 2*M[6] + 6*M[7] - 4*M[9] - M[10] - 2*M[11] + 4*M[13] + 2*M[14]);
686 w[8] = S * (1 + 2*M[1] + 4*M[2] + 6*M[4] + 6*M[5] - 4*M[6] -12*M[7] - 6*M[8] - 4*M[9] + 2*M[10] + 4*M[11] - 2*M[13] - M[14]);
687 w[9] = S * ( 2*M[6] + 6*M[7] + 6*M[8] + 2*M[9] - M[10] - 2*M[11] - 2*M[13] - M[14]);
688 w[10] = S * ( 2*M[9] - 2*M[13] - M[14]);
689 w[11] = S * ( 2*M[13] + M[14]);
690 } else if (totalOrder == 1) {
691 S *= (REAL) (1.0 / 6.0);
692
693 if (ds) {
694 w[0] = S * (-1 + 3*M[2] + 3*M[3] - 3*M[5] - 2*M[6] - 3*M[7] + M[9]);
695 w[1] = S * ( 1 - 3*M[2] - 6*M[3] + 3*M[5] + 4*M[6] + 6*M[7] - M[9]);
696 w[2] = S * ( 3*M[3] - 2*M[6] - 3*M[7] );
697 w[3] = S * (-2 + 6*M[1] + 3*M[2] - 6*M[3] - 6*M[4] + 2*M[6] + 3*M[7] - M[9]);
698 w[4] = S * ( -12*M[1] - 6*M[2] +12*M[3] +12*M[4] + 6*M[5] - 2*M[6] - 3*M[7] - M[9]);
699 w[5] = S * ( 2 + 6*M[1] + 3*M[2] - 6*M[3] - 6*M[4] - 6*M[5] - 2*M[6] - 3*M[7] + 2*M[9]);
700 w[6] = S * ( 2*M[6] + 3*M[7] );
701 w[7] = S * (-1 - 3*M[2] + 3*M[3] + 6*M[4] - 2*M[6] - 3*M[7] + 2*M[9]);
702 w[8] = S * ( 1 + 3*M[2] - 6*M[3] -12*M[4] - 3*M[5] + 4*M[6] + 6*M[7] - M[9]);
703 w[9] = S * ( 3*M[3] + 6*M[4] + 3*M[5] - 2*M[6] - 3*M[7] - M[9]);
704 w[10] = S * ( - M[9]);
705 w[11] = S * ( M[9]);
706 } else {
707 w[0] = S * (-2 + 3*M[1] + 6*M[2] - 6*M[4] - 6*M[5] - M[6] + 3*M[8] + 2*M[9]);
708 w[1] = S * (-1 - 3*M[1] + 6*M[4] + 3*M[5] + 2*M[6] - 3*M[8] - 2*M[9]);
709 w[2] = S * ( - M[6] );
710 w[3] = S * (-1 + 3*M[1] - 3*M[3] + 3*M[5] + M[6] - 3*M[8] - 2*M[9]);
711 w[4] = S * ( - 6*M[1] -12*M[2] + 6*M[3] +12*M[4] +12*M[5] - M[6] - 3*M[8] - 2*M[9]);
712 w[5] = S * ( 1 + 3*M[1] - 3*M[3] -12*M[4] - 6*M[5] - M[6] + 6*M[8] + 4*M[9]);
713 w[6] = S * ( + M[6] );
714 w[7] = S * ( 1 - 3*M[1] + 3*M[3] - 6*M[5] - M[6] + 6*M[8] + 4*M[9]);
715 w[8] = S * ( 2 + 3*M[1] + 6*M[2] - 6*M[3] - 6*M[4] - 6*M[5] + 2*M[6] - 3*M[8] - 2*M[9]);
716 w[9] = S * ( + 3*M[3] + 6*M[4] + 3*M[5] - M[6] - 3*M[8] - 2*M[9]);
717 w[10] = S * ( 3*M[5] - 3*M[8] - 2*M[9]);
718 w[11] = S * ( 3*M[8] + 2*M[9]);
719 }
720 } else if (totalOrder == 2) {
721 if (ds == 2) {
722 w[0] = S * ( + M[1] - M[3] - M[4]);
723 w[1] = S * ( - 2*M[1] + 2*M[3] + 2*M[4]);
724 w[2] = S * ( M[1] - M[3] - M[4]);
725 w[3] = S * ( 1 - 2*M[1] - M[2] + M[3] + M[4]);
726 w[4] = S * (-2 + 4*M[1] + 2*M[2] - M[3] - M[4]);
727 w[5] = S * ( 1 - 2*M[1] - M[2] - M[3] - M[4]);
728 w[6] = S * ( M[3] + M[4]);
729 w[7] = S * ( + M[1] + M[2] - M[3] - M[4]);
730 w[8] = S * ( - 2*M[1] - 2*M[2] + 2*M[3] + 2*M[4]);
731 w[9] = S * ( M[1] + M[2] - M[3] - M[4]);
732 w[10] = 0;
733 w[11] = 0;
734 } else if (dt == 2) {
735 w[0] = S * ( 1 - M[1] - 2*M[2] + M[4] + M[5]);
736 w[1] = S * ( + M[1] + M[2] - M[4] - M[5]);
737 w[2] = 0;
738 w[3] = S * ( + M[2] - M[4] - M[5]);
739 w[4] = S * (-2 + 2*M[1] + 4*M[2] - M[4] - M[5]);
740 w[5] = S * ( - 2*M[1] - 2*M[2] + 2*M[4] + 2*M[5]);
741 w[6] = 0;
742 w[7] = S * ( - 2*M[2] + 2*M[4] + 2*M[5]);
743 w[8] = S * ( 1 - M[1] - 2*M[2] - M[4] - M[5]);
744 w[9] = S * ( + M[1] + M[2] - M[4] - M[5]);
745 w[10] = S * ( M[2] - M[4] - M[5]);
746 w[11] = S * ( M[4] + M[5]);
747 } else {
748 S *= (REAL) (1.0 / 2.0);
749
750 w[0] = S * ( 1 - 2*M[2] - M[3] + M[5]);
751 w[1] = S * (-1 + 2*M[2] + 2*M[3] - M[5]);
752 w[2] = S * ( - M[3] );
753 w[3] = S * ( 1 - 2*M[1] + M[3] - M[5]);
754 w[4] = S * (-2 + 4*M[1] + 4*M[2] - M[3] - M[5]);
755 w[5] = S * ( 1 - 2*M[1] - 4*M[2] - M[3] + 2*M[5]);
756 w[6] = S * ( + M[3] );
757 w[7] = S * (-1 + 2*M[1] - M[3] + 2*M[5]);
758 w[8] = S * ( 1 - 4*M[1] - 2*M[2] + 2*M[3] - M[5]);
759 w[9] = S * ( + 2*M[1] + 2*M[2] - M[3] - M[5]);
760 w[10] = S * ( - M[5]);
761 w[11] = S * ( M[5]);
762 }
763 } else {
764 assert(totalOrder <= 2);
765 }
766 }
767
768 template <typename REAL>
769 void
adjustBoxSplineTriBoundaryWeights(int boundaryMask,REAL weights[])770 adjustBoxSplineTriBoundaryWeights(int boundaryMask, REAL weights[]) {
771
772 if (boundaryMask == 0) return;
773
774 //
775 // Determine boundary edges and vertices from the lower 3 and upper
776 // 2 bits of the 5-bit mask:
777 //
778 int upperBits = (boundaryMask >> 3) & 0x3;
779 int lowerBits = boundaryMask & 7;
780
781 int eBits = lowerBits;
782 int vBits = 0;
783
784 if (upperBits == 1) {
785 // Boundary vertices only:
786 vBits = eBits;
787 eBits = 0;
788 } else if (upperBits == 2) {
789 // Opposite vertex bit is edge bit rotated one to the right:
790 vBits = ((eBits & 1) << 2) | (eBits >> 1);
791 }
792
793 bool edge0IsBoundary = (eBits & 1) != 0;
794 bool edge1IsBoundary = (eBits & 2) != 0;
795 bool edge2IsBoundary = (eBits & 4) != 0;
796
797 //
798 // Adjust weights for the 4 boundary points and 3 interior points
799 // to account for the 3 phantom points adjacent to each
800 // boundary edge:
801 //
802 if (edge0IsBoundary) {
803 REAL w0 = weights[0];
804 if (edge2IsBoundary) {
805 // P0 = B1 + (B1 - I1)
806 weights[4] += w0;
807 weights[4] += w0;
808 weights[8] -= w0;
809 } else {
810 // P0 = B1 + (B0 - I0)
811 weights[4] += w0;
812 weights[3] += w0;
813 weights[7] -= w0;
814 }
815
816 // P1 = B1 + (B2 - I1)
817 REAL w1 = weights[1];
818 weights[4] += w1;
819 weights[5] += w1;
820 weights[8] -= w1;
821
822 REAL w2 = weights[2];
823 if (edge1IsBoundary) {
824 // P2 = B2 + (B2 - I1)
825 weights[5] += w2;
826 weights[5] += w2;
827 weights[8] -= w2;
828 } else {
829 // P2 = B2 + (B3 - I2)
830 weights[5] += w2;
831 weights[6] += w2;
832 weights[9] -= w2;
833 }
834 // Clear weights for the phantom points:
835 weights[0] = weights[1] = weights[2] = 0.0f;
836 }
837 if (edge1IsBoundary) {
838 REAL w0 = weights[6];
839 if (edge0IsBoundary) {
840 // P0 = B1 + (B1 - I1)
841 weights[5] += w0;
842 weights[5] += w0;
843 weights[4] -= w0;
844 } else {
845 // P0 = B1 + (B0 - I0)
846 weights[5] += w0;
847 weights[2] += w0;
848 weights[1] -= w0;
849 }
850
851 // P1 = B1 + (B2 - I1)
852 REAL w1 = weights[9];
853 weights[5] += w1;
854 weights[8] += w1;
855 weights[4] -= w1;
856
857 REAL w2 = weights[11];
858 if (edge2IsBoundary) {
859 // P2 = B2 + (B2 - I1)
860 weights[8] += w2;
861 weights[8] += w2;
862 weights[4] -= w2;
863 } else {
864 // P2 = B2 + (B3 - I2)
865 weights[8] += w2;
866 weights[10] += w2;
867 weights[7] -= w2;
868 }
869 // Clear weights for the phantom points:
870 weights[6] = weights[9] = weights[11] = 0.0f;
871 }
872 if (edge2IsBoundary) {
873 REAL w0 = weights[10];
874 if (edge1IsBoundary) {
875 // P0 = B1 + (B1 - I1)
876 weights[8] += w0;
877 weights[8] += w0;
878 weights[5] -= w0;
879 } else {
880 // P0 = B1 + (B0 - I0)
881 weights[8] += w0;
882 weights[11] += w0;
883 weights[9] -= w0;
884 }
885
886 // P1 = B1 + (B2 - I1)
887 REAL w1 = weights[7];
888 weights[8] += w1;
889 weights[4] += w1;
890 weights[5] -= w1;
891
892 REAL w2 = weights[3];
893 if (edge0IsBoundary) {
894 // P2 = B2 + (B2 - I1)
895 weights[4] += w2;
896 weights[4] += w2;
897 weights[5] -= w2;
898 } else {
899 // P2 = B2 + (B3 - I2)
900 weights[4] += w2;
901 weights[0] += w2;
902 weights[1] -= w2;
903 }
904 // Clear weights for the phantom points:
905 weights[10] = weights[7] = weights[3] = 0.0f;
906 }
907
908 //
909 // Adjust weights for the 3 boundary points and the 2 interior
910 // points to account for the 2 phantom points adjacent to
911 // each boundary vertex:
912 //
913 if ((vBits & 1) != 0) {
914 // P0 = B1 + (B0 - I0)
915 REAL w0 = weights[3];
916 weights[4] += w0;
917 weights[7] += w0;
918 weights[8] -= w0;
919
920 // P1 = B1 + (B2 - I1)
921 REAL w1 = weights[0];
922 weights[4] += w1;
923 weights[1] += w1;
924 weights[5] -= w1;
925
926 // Clear weights for the phantom points:
927 weights[3] = weights[0] = 0.0f;
928 }
929 if ((vBits & 2) != 0) {
930 // P0 = B1 + (B0 - I0)
931 REAL w0 = weights[2];
932 weights[5] += w0;
933 weights[1] += w0;
934 weights[4] -= w0;
935
936 // P1 = B1 + (B2 - I1)
937 REAL w1 = weights[6];
938 weights[5] += w1;
939 weights[9] += w1;
940 weights[8] -= w1;
941
942 // Clear weights for the phantom points:
943 weights[2] = weights[6] = 0.0f;
944 }
945 if ((vBits & 4) != 0) {
946 // P0 = B1 + (B0 - I0)
947 REAL w0 = weights[11];
948 weights[8] += w0;
949 weights[9] += w0;
950 weights[5] -= w0;
951
952 // P1 = B1 + (B2 - I1)
953 REAL w1 = weights[10];
954 weights[8] += w1;
955 weights[7] += w1;
956 weights[4] -= w1;
957
958 // Clear weights for the phantom points:
959 weights[11] = weights[10] = 0.0f;
960 }
961 }
962
963 template <typename REAL>
964 void
boundBasisBoxSplineTri(int boundary,REAL wP[12],REAL wDs[12],REAL wDt[12],REAL wDss[12],REAL wDst[12],REAL wDtt[12])965 boundBasisBoxSplineTri(int boundary,
966 REAL wP[12], REAL wDs[12], REAL wDt[12],
967 REAL wDss[12], REAL wDst[12], REAL wDtt[12]) {
968
969 if (wP) {
970 adjustBoxSplineTriBoundaryWeights(boundary, wP);
971 }
972 if (wDs && wDt) {
973 adjustBoxSplineTriBoundaryWeights(boundary, wDs);
974 adjustBoxSplineTriBoundaryWeights(boundary, wDt);
975
976 if (wDss && wDst && wDtt) {
977 adjustBoxSplineTriBoundaryWeights(boundary, wDss);
978 adjustBoxSplineTriBoundaryWeights(boundary, wDst);
979 adjustBoxSplineTriBoundaryWeights(boundary, wDtt);
980 }
981 }
982 }
983 } // namespace
984
985 template <typename REAL>
EvalBasisBoxSplineTri(REAL s,REAL t,REAL wP[12],REAL wDs[12],REAL wDt[12],REAL wDss[12],REAL wDst[12],REAL wDtt[12])986 int EvalBasisBoxSplineTri(REAL s, REAL t,
987 REAL wP[12], REAL wDs[12], REAL wDt[12],
988 REAL wDss[12], REAL wDst[12], REAL wDtt[12]) {
989
990 REAL stMonomials[15];
991 evalBivariateMonomialsQuartic(s, t, stMonomials);
992
993 if (wP) {
994 evalBoxSplineTriDerivWeights<REAL>(stMonomials, 0, 0, wP);
995 }
996 if (wDs && wDt) {
997 evalBoxSplineTriDerivWeights(stMonomials, 1, 0, wDs);
998 evalBoxSplineTriDerivWeights(stMonomials, 0, 1, wDt);
999
1000 if (wDss && wDst && wDtt) {
1001 evalBoxSplineTriDerivWeights(stMonomials, 2, 0, wDss);
1002 evalBoxSplineTriDerivWeights(stMonomials, 1, 1, wDst);
1003 evalBoxSplineTriDerivWeights(stMonomials, 0, 2, wDtt);
1004 }
1005 }
1006 return 12;
1007 }
1008
1009
1010 //
1011 // Quartic Bezier triangle:
1012 //
1013 // The regular patch for Loop subdivision is a quartic triangular Box spline
1014 // with cubic boundaries. So we need a quartic Bezier patch to represent it
1015 // faithfully.
1016 //
1017 // The formulae for the 15 quartic basis functions are:
1018 //
1019 // 4! i j k
1020 // B (u,v,w) = ------- * (u * v * w )
1021 // ijk i!j!k!
1022 //
1023 // for each i + j + k = 4. The control points (P) and correspondingly
1024 // labeled p<i,j,k> are oriented as follows:
1025 //
1026 // P14 p040
1027 // P12 P13 p031 p130
1028 // P9 P10 P11 p022 p121 p220
1029 // P5 P6 P7 P8 p013 p112 p211 p310
1030 // P0 P1 P2 P3 P4 p004 p103 p202 p301 p400
1031 //
1032 namespace {
1033 template <typename REAL>
1034 void
evalBezierTriDerivWeights(REAL s,REAL t,int ds,int dt,REAL wB[])1035 evalBezierTriDerivWeights(REAL s, REAL t, int ds, int dt, REAL wB[]) {
1036
1037 REAL u = s;
1038 REAL v = t;
1039 REAL w = 1 - u - v;
1040
1041 REAL uu = u * u;
1042 REAL vv = v * v;
1043 REAL ww = w * w;
1044
1045 REAL uv = u * v;
1046 REAL vw = v * w;
1047 REAL uw = u * w;
1048
1049 int totalOrder = ds + dt;
1050 if (totalOrder == 0) {
1051 wB[0] = ww * ww;
1052 wB[1] = 4 * uw * ww;
1053 wB[2] = 6 * uw * uw;
1054 wB[3] = 4 * uw * uu;
1055 wB[4] = uu * uu;
1056 wB[5] = 4 * vw * ww;
1057 wB[6] = 12 * ww * uv;
1058 wB[7] = 12 * uu * vw;
1059 wB[8] = 4 * uv * uu;
1060 wB[9] = 6 * vw * vw;
1061 wB[10] = 12 * vv * uw;
1062 wB[11] = 6 * uv * uv;
1063 wB[12] = 4 * vw * vv;
1064 wB[13] = 4 * uv * vv;
1065 wB[14] = vv * vv;
1066 } else if (totalOrder == 1) {
1067 if (ds == 1) {
1068 wB[0] = -4 * ww * w;
1069 wB[1] = 4 * ww * (w - 3 * u);
1070 wB[2] = 12 * uw * (w - u);
1071 wB[3] = 4 * uu * (3 * w - u);
1072 wB[4] = 4 * uu * u;
1073 wB[5] = -12 * vw * w;
1074 wB[6] = 12 * vw * (w - 2 * u);
1075 wB[7] = 12 * uv * (2 * w - u);
1076 wB[8] = 12 * uv * u;
1077 wB[9] = -12 * vv * w;
1078 wB[10] = 12 * vv * (w - u);
1079 wB[11] = 12 * vv * u;
1080 wB[12] = -4 * vv * v;
1081 wB[13] = 4 * vv * v;
1082 wB[14] = 0;
1083 } else {
1084 wB[0] = -4 * ww * w;
1085 wB[1] = -12 * ww * u;
1086 wB[2] = -12 * uu * w;
1087 wB[3] = -4 * uu * u;
1088 wB[4] = 0;
1089 wB[5] = 4 * ww * (w - 3 * v);
1090 wB[6] = 12 * uw * (w - 2 * v);
1091 wB[7] = 12 * uu * (w - v);
1092 wB[8] = 4 * uu * u;
1093 wB[9] = 12 * vw * (w - v);
1094 wB[10] = 12 * uv * (2 * w - v);
1095 wB[11] = 12 * uv * u;;
1096 wB[12] = 4 * vv * (3 * w - v);
1097 wB[13] = 12 * vv * u;
1098 wB[14] = 4 * vv * v;
1099 }
1100 } else if (totalOrder == 2) {
1101 if (ds == 2) {
1102 wB[0] = 12 * ww;
1103 wB[1] = 24 * (uw - ww);
1104 wB[2] = 12 * (uu - 4 * uw + ww);
1105 wB[3] = 24 * (uw - uu);
1106 wB[4] = 12 * uu;
1107 wB[5] = 24 * vw;
1108 wB[6] = 24 * (uv - 2 * vw);
1109 wB[7] = 24 * (vw - 2 * uv);
1110 wB[8] = 24 * uv;
1111 wB[9] = 12 * vv;
1112 wB[10] = -24 * vv;
1113 wB[11] = 12 * vv;
1114 wB[12] = 0;
1115 wB[13] = 0;
1116 wB[14] = 0;
1117 } else if (dt == 2) {
1118 wB[0] = 12 * ww;
1119 wB[1] = 24 * uw;
1120 wB[2] = 12 * uu;
1121 wB[3] = 0;
1122 wB[4] = 0;
1123 wB[5] = 24 * (vw - ww);
1124 wB[6] = 24 * (uv - 2 * uw);
1125 wB[7] = -24 * uu;
1126 wB[8] = 0;
1127 wB[9] = 12 * (vv - 4 * vw + ww);
1128 wB[10] = 24 * (uw - 2 * uv);
1129 wB[11] = 12 * uu;
1130 wB[12] = 24 * (vw - vv);
1131 wB[13] = 24 * uv;
1132 wB[14] = 12 * vv;
1133 } else {
1134 wB[0] = 12 * ww;
1135 wB[3] = -12 * uu;
1136 wB[13] = 12 * vv;
1137 wB[11] = 24 * uv;
1138 wB[1] = 24 * uw - wB[0];
1139 wB[2] = -24 * uw - wB[3];
1140 wB[5] = 24 * vw - wB[0];
1141 wB[6] = -24 * vw + wB[11] - wB[1];
1142 wB[8] = - wB[3];
1143 wB[7] = -(wB[11] + wB[2]);
1144 wB[9] = wB[13] - wB[5] - wB[0];
1145 wB[10] = -(wB[9] + wB[11]);
1146 wB[12] = - wB[13];
1147 wB[4] = 0;
1148 wB[14] = 0;
1149 }
1150 } else {
1151 assert(totalOrder <= 2);
1152 }
1153 }
1154 } // end namespace
1155
1156 template <typename REAL>
1157 int
EvalBasisBezierTri(REAL s,REAL t,REAL wP[15],REAL wDs[15],REAL wDt[15],REAL wDss[15],REAL wDst[15],REAL wDtt[15])1158 EvalBasisBezierTri(REAL s, REAL t,
1159 REAL wP[15], REAL wDs[15], REAL wDt[15],
1160 REAL wDss[15], REAL wDst[15], REAL wDtt[15]) {
1161
1162 if (wP) {
1163 evalBezierTriDerivWeights<REAL>(s, t, 0, 0, wP);
1164 }
1165 if (wDs && wDt) {
1166 evalBezierTriDerivWeights(s, t, 1, 0, wDs);
1167 evalBezierTriDerivWeights(s, t, 0, 1, wDt);
1168
1169 if (wDss && wDst && wDtt) {
1170 evalBezierTriDerivWeights(s, t, 2, 0, wDss);
1171 evalBezierTriDerivWeights(s, t, 1, 1, wDst);
1172 evalBezierTriDerivWeights(s, t, 0, 2, wDtt);
1173 }
1174 }
1175 return 15;
1176 }
1177
1178
1179 //
1180 // Quartic Gregory triangle:
1181 //
1182 // The 18-point quaritic Gregory patch is an extension of the 15-point
1183 // quartic Bezier triangle with the 3 interior points of the Bezier patch
1184 // replaced with pairs of points (face points -- fi+ and fi-) that are
1185 // rationally combined.
1186 //
1187 // The point ordering of Gregory patches deviates considerably from the
1188 // BSpline and Bezier patches by grouping the 5 points at each corner and
1189 // ordering the groups by corner index. This is consistent with the cubic
1190 // Gregory quad patch.
1191 //
1192 // The 3 additional quartic boundary points are currently appended to these
1193 // 3 groups of 5 control points. In contrast to the 5 points associated
1194 // with each corner, these 3 points are more associated with the edge
1195 // between the corner vertices and are equally weighted between the two.
1196 //
1197 namespace {
1198 //
1199 // Expanding a set of 15 Bezier basis functions for the 6 (3 pairs) of
1200 // rational weights for the 18 Gregory basis functions:
1201 //
1202 template <typename REAL>
1203 void
convertBezierWeightsToGregory(REAL const wB[15],REAL const rG[6],REAL wG[18])1204 convertBezierWeightsToGregory(REAL const wB[15], REAL const rG[6], REAL wG[18]) {
1205
1206 wG[0] = wB[0];
1207 wG[1] = wB[1];
1208 wG[2] = wB[5];
1209 wG[3] = wB[6] * rG[0];
1210 wG[4] = wB[6] * rG[1];
1211
1212 wG[5] = wB[4];
1213 wG[6] = wB[8];
1214 wG[7] = wB[3];
1215 wG[8] = wB[7] * rG[2];
1216 wG[9] = wB[7] * rG[3];
1217
1218 wG[10] = wB[14];
1219 wG[11] = wB[12];
1220 wG[12] = wB[13];
1221 wG[13] = wB[10] * rG[4];
1222 wG[14] = wB[10] * rG[5];
1223
1224 wG[15] = wB[2];
1225 wG[16] = wB[11];
1226 wG[17] = wB[9];
1227 }
1228 } // end namespace
1229
1230 template <typename REAL>
1231 int
EvalBasisGregoryTri(REAL s,REAL t,REAL wP[18],REAL wDs[18],REAL wDt[18],REAL wDss[18],REAL wDst[18],REAL wDtt[18])1232 EvalBasisGregoryTri(REAL s, REAL t,
1233 REAL wP[18], REAL wDs[18], REAL wDt[18],
1234 REAL wDss[18], REAL wDst[18], REAL wDtt[18]) {
1235
1236 //
1237 // Bezier basis functions are denoted with B while the rational multipliers for the
1238 // interior points will be denoted G -- so we have B(s,t) and G(s,t) (though we
1239 // switch to barycentric (u,v,w) briefly to compute G)
1240 //
1241 REAL BP[15], BDs[15], BDt[15], BDss[15], BDst[15], BDtt[15];
1242
1243 REAL G[6] = { 1.0f, 0.0f, 1.0f, 0.0f, 1.0f, 0.0f };
1244 REAL u = s;
1245 REAL v = t;
1246 REAL w = 1 - u - v;
1247
1248 if ((u + v) > 0) {
1249 G[0] = u / (u + v);
1250 G[1] = v / (u + v);
1251 }
1252 if ((v + w) > 0) {
1253 G[2] = v / (v + w);
1254 G[3] = w / (v + w);
1255 }
1256 if ((w + u) > 0) {
1257 G[4] = w / (w + u);
1258 G[5] = u / (w + u);
1259 }
1260
1261 //
1262 // Compute Bezier basis functions and convert, adjusting interior points:
1263 //
1264 if (wP) {
1265 evalBezierTriDerivWeights<REAL>(s, t, 0, 0, BP);
1266 convertBezierWeightsToGregory(BP, G, wP);
1267 }
1268 if (wDs && wDt) {
1269 // TBD -- ifdef OPENSUBDIV_GREGORY_EVAL_TRUE_DERIVATIVES
1270
1271 evalBezierTriDerivWeights(s, t, 1, 0, BDs);
1272 evalBezierTriDerivWeights(s, t, 0, 1, BDt);
1273
1274 convertBezierWeightsToGregory(BDs, G, wDs);
1275 convertBezierWeightsToGregory(BDt, G, wDt);
1276
1277 if (wDss && wDst && wDtt) {
1278 evalBezierTriDerivWeights(s, t, 2, 0, BDss);
1279 evalBezierTriDerivWeights(s, t, 1, 1, BDst);
1280 evalBezierTriDerivWeights(s, t, 0, 2, BDtt);
1281
1282 convertBezierWeightsToGregory(BDss, G, wDss);
1283 convertBezierWeightsToGregory(BDst, G, wDst);
1284 convertBezierWeightsToGregory(BDtt, G, wDtt);
1285 }
1286 }
1287 return 18;
1288 }
1289
1290 //
1291 // Higher level basis evaluation functions that deal with parameterization and
1292 // boundary issues (reflected in PatchParam) for all patch types:
1293 //
1294 template <typename REAL>
1295 int
EvaluatePatchBasisNormalized(int patchType,PatchParam const & param,REAL s,REAL t,REAL wP[],REAL wDs[],REAL wDt[],REAL wDss[],REAL wDst[],REAL wDtt[])1296 EvaluatePatchBasisNormalized(int patchType, PatchParam const & param, REAL s, REAL t,
1297 REAL wP[], REAL wDs[], REAL wDt[],
1298 REAL wDss[], REAL wDst[], REAL wDtt[]) {
1299
1300 int boundaryMask = param.GetBoundary();
1301
1302 int nPoints = 0;
1303 if (patchType == PatchDescriptor::REGULAR) {
1304 nPoints = EvalBasisBSpline(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
1305 if (boundaryMask) {
1306 boundBasisBSpline(boundaryMask, wP, wDs, wDt, wDss, wDst, wDtt);
1307 }
1308 } else if (patchType == PatchDescriptor::LOOP) {
1309 nPoints = EvalBasisBoxSplineTri(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
1310 if (boundaryMask) {
1311 boundBasisBoxSplineTri(boundaryMask, wP, wDs, wDt, wDss, wDst, wDtt);
1312 }
1313 } else if (patchType == PatchDescriptor::GREGORY_BASIS) {
1314 nPoints = EvalBasisGregory(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
1315 } else if (patchType == PatchDescriptor::GREGORY_TRIANGLE) {
1316 nPoints = EvalBasisGregoryTri(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
1317 } else if (patchType == PatchDescriptor::QUADS) {
1318 nPoints = EvalBasisLinear(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
1319 } else if (patchType == PatchDescriptor::TRIANGLES) {
1320 nPoints = EvalBasisLinearTri(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
1321 } else {
1322 assert(0);
1323 }
1324 return nPoints;
1325 }
1326
1327 template <typename REAL>
1328 int
EvaluatePatchBasis(int patchType,PatchParam const & param,REAL s,REAL t,REAL wP[],REAL wDs[],REAL wDt[],REAL wDss[],REAL wDst[],REAL wDtt[])1329 EvaluatePatchBasis(int patchType, PatchParam const & param, REAL s, REAL t,
1330 REAL wP[], REAL wDs[], REAL wDt[],
1331 REAL wDss[], REAL wDst[], REAL wDtt[]) {
1332
1333 REAL derivSign = 1.0f;
1334
1335 if ((patchType == PatchDescriptor::LOOP) ||
1336 (patchType == PatchDescriptor::GREGORY_TRIANGLE) ||
1337 (patchType == PatchDescriptor::TRIANGLES)) {
1338 param.NormalizeTriangle(s, t);
1339 if (param.IsTriangleRotated()) {
1340 derivSign = -1.0f;
1341 }
1342 } else {
1343 param.Normalize(s, t);
1344 }
1345
1346 int nPoints = EvaluatePatchBasisNormalized(
1347 patchType, param, s, t, wP, wDs, wDt, wDss, wDst, wDtt);
1348
1349 if (wDs && wDt) {
1350 REAL d1Scale = derivSign * (REAL)(1 << param.GetDepth());
1351
1352 for (int i = 0; i < nPoints; ++i) {
1353 wDs[i] *= d1Scale;
1354 wDt[i] *= d1Scale;
1355 }
1356
1357 if (wDss && wDst && wDtt) {
1358 REAL d2Scale = derivSign * d1Scale * d1Scale;
1359
1360 for (int i = 0; i < nPoints; ++i) {
1361 wDss[i] *= d2Scale;
1362 wDst[i] *= d2Scale;
1363 wDtt[i] *= d2Scale;
1364 }
1365 }
1366 }
1367 return nPoints;
1368 }
1369
1370 //
1371 // Explicit float and double instantiations:
1372 //
1373 template int EvaluatePatchBasisNormalized<float>(int patchType, PatchParam const & param,
1374 float s, float t, float wP[], float wDs[], float wDt[], float wDss[], float wDst[], float wDtt[]);
1375 template int EvaluatePatchBasis<float>(int patchType, PatchParam const & param,
1376 float s, float t, float wP[], float wDs[], float wDt[], float wDss[], float wDst[], float wDtt[]);
1377
1378 template int EvaluatePatchBasisNormalized<double>(int patchType, PatchParam const & param,
1379 double s, double t, double wP[], double wDs[], double wDt[], double wDss[], double wDst[], double wDtt[]);
1380 template int EvaluatePatchBasis<double>(int patchType, PatchParam const & param,
1381 double s, double t, double wP[], double wDs[], double wDt[], double wDss[], double wDst[], double wDtt[]);
1382
1383 //
1384 // Most basis evaluation functions are implicitly instantiated above -- Bezier
1385 // require explicit instantiation as they are not invoked via a patch type:
1386 //
1387 template int EvalBasisBezier<float>(float s, float t,
1388 float wP[16], float wDs[16], float wDt[16], float wDss[16], float wDst[16], float wDtt[16]);
1389 template int EvalBasisBezierTri<float>(float s, float t,
1390 float wP[12], float wDs[12], float wDt[12], float wDss[12], float wDst[12], float wDtt[12]);
1391
1392 template int EvalBasisBezier<double>(double s, double t,
1393 double wP[16], double wDs[16], double wDt[16], double wDss[16], double wDst[16], double wDtt[16]);
1394 template int EvalBasisBezierTri<double>(double s, double t,
1395 double wP[12], double wDs[12], double wDt[12], double wDss[12], double wDst[12], double wDtt[12]);
1396
1397 } // end namespace internal
1398 } // end namespace Far
1399
1400 } // end namespace OPENSUBDIV_VERSION
1401 } // end namespace OpenSubdiv
1402