OpenSubdiv/opensubdiv/far/patchBasis.cpp
barry 484c2746d1 Minor refactoring of internal functions in far/patchBasis.*:
- individual basis functions now purely normalized with no PatchParam
    - two new higher level functions deal with patch type and PatchParam
    - updated Far::PatchTable and Osd evaluators that used old methods
2018-11-06 13:29:16 -08:00

1372 lines
50 KiB
C++

//
// Copyright 2013 Pixar
//
// Licensed under the Apache License, Version 2.0 (the "Apache License")
// with the following modification; you may not use this file except in
// compliance with the Apache License and the following modification to it:
// Section 6. Trademarks. is deleted and replaced with:
//
// 6. Trademarks. This License does not grant permission to use the trade
// names, trademarks, service marks, or product names of the Licensor
// and its affiliates, except as required to comply with Section 4(c) of
// the License and to reproduce the content of the NOTICE file.
//
// You may obtain a copy of the Apache License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the Apache License with the above modification is
// distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, either express or implied. See the Apache License for the specific
// language governing permissions and limitations under the Apache License.
//
#include "../far/patchBasis.h"
#include "../far/patchDescriptor.h"
#include <cassert>
#include <cstring>
#include <cmath>
#include <cstdio>
namespace OpenSubdiv {
namespace OPENSUBDIV_VERSION {
namespace Far {
namespace internal {
//
// Basis support for quadrilateral patches:
//
// Quadrilateral patches are parameterized in terms of (s,t) as follows:
//
// (1,0) *---------* (1,1)
// | 3 2 |
// t | |
// | |
// | 0 1 |
// (0,0) *---------* (1,0)
// s
//
//
// Simple bilinear quad:
//
template <typename REAL>
int
EvalBasisLinear(REAL s, REAL t,
REAL wP[4], REAL wDs[4], REAL wDt[4],
REAL wDss[4], REAL wDst[4], REAL wDtt[4]) {
REAL sC = 1.0f - s;
REAL tC = 1.0f - t;
if (wP) {
wP[0] = sC * tC;
wP[1] = s * tC;
wP[2] = s * t;
wP[3] = sC * t;
}
if (wDs && wDt) {
wDs[0] = -tC;
wDs[1] = tC;
wDs[2] = t;
wDs[3] = -t;
wDt[0] = -sC;
wDt[1] = -s;
wDt[2] = s;
wDt[3] = sC;
if (wDss && wDst && wDtt) {
for(int i = 0; i < 4; ++i) {
wDss[i] = 0.0f;
wDtt[i] = 0.0f;
}
wDst[0] = 1.0f;
wDst[1] = -1.0f;
wDst[2] = -1.0f;
wDst[3] = 1.0f;
}
}
return 4;
}
//
// Bicubic BSpline patch:
//
// 12-----13------14-----15
// | | | |
// | | | |
// 8------9------10-----11
// | | t | |
// | | s | |
// 4------5------6------7
// | | | |
// | | | |
// O------1------2------3
//
// The basis if a bicubic BSpline patch is a tensor product, which we make
// use of here by evaluating, differentiating and combining basis functions
// in each of the two parametric directions.
//
// Not all 16 points will be present. The boundary mask indicates boundary
// edges beyond which phantom points are implicitly extrapolated. Weights
// for missing points are set to zero while those contributing to their
// implicit extrapolation will be adjusted.
//
namespace {
//
// Cubic BSpline curve basis evaluation:
//
template <typename REAL>
void
evalBSplineCurve(REAL t, REAL wP[4], REAL wDP[4], REAL wDP2[4]) {
REAL const one6th = (REAL)(1.0 / 6.0);
REAL t2 = t * t;
REAL t3 = t * t2;
wP[0] = one6th * (1.0f - 3.0f*(t - t2) - t3);
wP[1] = one6th * (4.0f - 6.0f*t2 + 3.0f*t3);
wP[2] = one6th * (1.0f + 3.0f*(t + t2 - t3));
wP[3] = one6th * ( t3);
if (wDP) {
wDP[0] = -0.5f*t2 + t - 0.5f;
wDP[1] = 1.5f*t2 - 2.0f*t;
wDP[2] = -1.5f*t2 + t + 0.5f;
wDP[3] = 0.5f*t2;
}
if (wDP2) {
wDP2[0] = - t + 1.0f;
wDP2[1] = 3.0f * t - 2.0f;
wDP2[2] = -3.0f * t + 1.0f;
wDP2[3] = t;
}
}
//
// Weight adjustments to account for phantom end points:
//
template <typename REAL>
void
adjustBSplineBoundaryWeights(int boundary, REAL w[16]) {
if ((boundary & 1) != 0) {
for (int i = 0; i < 4; ++i) {
w[i + 8] -= w[i + 0];
w[i + 4] += w[i + 0] * 2.0f;
w[i + 0] = 0.0f;
}
}
if ((boundary & 2) != 0) {
for (int i = 0; i < 16; i += 4) {
w[i + 1] -= w[i + 3];
w[i + 2] += w[i + 3] * 2.0f;
w[i + 3] = 0.0f;
}
}
if ((boundary & 4) != 0) {
for (int i = 0; i < 4; ++i) {
w[i + 4] -= w[i + 12];
w[i + 8] += w[i + 12] * 2.0f;
w[i + 12] = 0.0f;
}
}
if ((boundary & 8) != 0) {
for (int i = 0; i < 16; i += 4) {
w[i + 2] -= w[i + 0];
w[i + 1] += w[i + 0] * 2.0f;
w[i + 0] = 0.0f;
}
}
}
template <typename REAL>
void
boundBasisBSpline(int boundary,
REAL wP[16], REAL wDs[16], REAL wDt[16],
REAL wDss[16], REAL wDst[16], REAL wDtt[16]) {
if (wP) {
adjustBSplineBoundaryWeights(boundary, wP);
}
if (wDs && wDt) {
adjustBSplineBoundaryWeights(boundary, wDs);
adjustBSplineBoundaryWeights(boundary, wDt);
if (wDss && wDst && wDtt) {
adjustBSplineBoundaryWeights(boundary, wDss);
adjustBSplineBoundaryWeights(boundary, wDst);
adjustBSplineBoundaryWeights(boundary, wDtt);
}
}
}
} // end namespace
template <typename REAL>
int
EvalBasisBSpline(REAL s, REAL t,
REAL wP[16], REAL wDs[16], REAL wDt[16],
REAL wDss[16], REAL wDst[16], REAL wDtt[16]) {
REAL sWeights[4], tWeights[4], dsWeights[4], dtWeights[4], dssWeights[4], dttWeights[4];
evalBSplineCurve(s, wP ? sWeights : 0, wDs ? dsWeights : 0, wDss ? dssWeights : 0);
evalBSplineCurve(t, wP ? tWeights : 0, wDt ? dtWeights : 0, wDtt ? dttWeights : 0);
if (wP) {
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
wP[4*i+j] = sWeights[j] * tWeights[i];
}
}
}
if (wDs && wDt) {
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
wDs[4*i+j] = dsWeights[j] * tWeights[i];
wDt[4*i+j] = sWeights[j] * dtWeights[i];
}
}
if (wDss && wDst && wDtt) {
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
wDss[4*i+j] = dssWeights[j] * tWeights[i];
wDst[4*i+j] = dsWeights[j] * dtWeights[i];
wDtt[4*i+j] = sWeights[j] * dttWeights[i];
}
}
}
}
return 16;
}
//
// Bicubic Bezier patch:
//
// 12-----13------14-----15
// | | | |
// | | | |
// 8------9------10-----11
// | | | |
// | | | |
// 4------5------6------7
// | t | | |
// | s | | |
// O------1------2------3
//
// As was the case with the BSpline patch, a bicubic Bezier patch can also
// make use of its tensor product property by evaluating, differentiating
// and combining basis functions in each of the two parametric directions.
//
namespace {
//
// Cubic Bezier curve basis evaluation:
//
template <typename REAL>
void
evalBezierCurve(REAL t, REAL wP[4], REAL wDP[4], REAL wDP2[4]) {
// The four uniform cubic Bezier basis functions (in terms of t and its
// complement tC) evaluated at t:
REAL t2 = t*t;
REAL tC = 1.0f - t;
REAL tC2 = tC * tC;
wP[0] = tC2 * tC;
wP[1] = tC2 * t * 3.0f;
wP[2] = t2 * tC * 3.0f;
wP[3] = t2 * t;
// Derivatives of the above four basis functions at t:
if (wDP) {
wDP[0] = -3.0f * tC2;
wDP[1] = 9.0f * t2 - 12.0f * t + 3.0f;
wDP[2] = -9.0f * t2 + 6.0f * t;
wDP[3] = 3.0f * t2;
}
// Second derivatives of the basis functions at t:
if (wDP2) {
wDP2[0] = 6.0f * tC;
wDP2[1] = 18.0f * t - 12.0f;
wDP2[2] = -18.0f * t + 6.0f;
wDP2[3] = 6.0f * t;
}
}
} // end namespace
template <typename REAL>
int
EvalBasisBezier(REAL s, REAL t,
REAL wP[16], REAL wDs[16], REAL wDt[16],
REAL wDss[16], REAL wDst[16], REAL wDtt[16]) {
REAL sWeights[4], tWeights[4], dsWeights[4], dtWeights[4], dssWeights[4], dttWeights[4];
evalBezierCurve(s, wP ? sWeights : 0, wDs ? dsWeights : 0, wDss ? dssWeights : 0);
evalBezierCurve(t, wP ? tWeights : 0, wDt ? dtWeights : 0, wDtt ? dttWeights : 0);
if (wP) {
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
wP[4*i+j] = sWeights[j] * tWeights[i];
}
}
}
if (wDs && wDt) {
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
wDs[4*i+j] = dsWeights[j] * tWeights[i];
wDt[4*i+j] = sWeights[j] * dtWeights[i];
}
}
if (wDss && wDst && wDtt) {
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
wDss[4*i+j] = dssWeights[j] * tWeights[i];
wDst[4*i+j] = dsWeights[j] * dtWeights[i];
wDtt[4*i+j] = sWeights[j] * dttWeights[i];
}
}
}
}
return 16;
}
//
// Cubic Gregory patch:
//
// P3 e3- e2+ P2
// 15------17-------11--------10
// | | | |
// | | | |
// | | f3- | f2+ |
// | 19 13 |
// e3+ 16-----18 14-----12 e2-
// | f3+ f2- |
// | |
// | |
// | f0- f1+ |
// e0- 2------4 8------6 e1+
// | 3 9 |
// | | f0+ | f1- |
// | t | | |
// | s | | |
// O--------1--------7--------5
// P0 e0+ e1- P1
//
// The 20-point cubic Gregory patch is an extension of the 16-point bicubic
// Bezier patch with the 4 interior points of the Bezier patch replaced with
// pairs of points (face points -- fi+ and fi-) that are rationally combined.
//
// The point ordering of the Gregory patch deviates considerably from the
// BSpline and Bezier patches by grouping the 5 points at each corner and
// ordering the groups by corner index.
//
template <typename REAL>
int
EvalBasisGregory(REAL s, REAL t,
REAL point[20], REAL wDs[20], REAL wDt[20],
REAL wDss[20], REAL wDst[20], REAL wDtt[20]) {
// Indices of boundary and interior points and their corresponding Bezier points
// (this can be reduced with more direct indexing and unrolling of loops):
//
static int const boundaryGregory[12] = { 0, 1, 7, 5, 2, 6, 16, 12, 15, 17, 11, 10 };
static int const boundaryBezSCol[12] = { 0, 1, 2, 3, 0, 3, 0, 3, 0, 1, 2, 3 };
static int const boundaryBezTRow[12] = { 0, 0, 0, 0, 1, 1, 2, 2, 3, 3, 3, 3 };
static int const interiorGregory[8] = { 3, 4, 8, 9, 13, 14, 18, 19 };
static int const interiorBezSCol[8] = { 1, 1, 2, 2, 2, 2, 1, 1 };
static int const interiorBezTRow[8] = { 1, 1, 1, 1, 2, 2, 2, 2 };
//
// Bezier basis functions are denoted with B while the rational multipliers for the
// interior points will be denoted G -- so we have B(s), B(t) and G(s,t):
//
// Directional Bezier basis functions B at s and t:
REAL Bs[4], Bds[4], Bdss[4];
REAL Bt[4], Bdt[4], Bdtt[4];
evalBezierCurve(s, Bs, wDs ? Bds : 0, wDss ? Bdss : 0);
evalBezierCurve(t, Bt, wDt ? Bdt : 0, wDtt ? Bdtt : 0);
// Rational multipliers G at s and t:
REAL sC = 1.0f - s;
REAL tC = 1.0f - t;
// Use <= here to avoid compiler warnings -- the sums should always be non-negative:
REAL df0 = s + t; df0 = (df0 <= 0.0f) ? (REAL)1.0f : (1.0f / df0);
REAL df1 = sC + t; df1 = (df1 <= 0.0f) ? (REAL)1.0f : (1.0f / df1);
REAL df2 = sC + tC; df2 = (df2 <= 0.0f) ? (REAL)1.0f : (1.0f / df2);
REAL df3 = s + tC; df3 = (df3 <= 0.0f) ? (REAL)1.0f : (1.0f / df3);
REAL G[8] = { s*df0, t*df0, t*df1, sC*df1, sC*df2, tC*df2, tC*df3, s*df3 };
// Combined weights for boundary and interior points:
for (int i = 0; i < 12; ++i) {
point[boundaryGregory[i]] = Bs[boundaryBezSCol[i]] * Bt[boundaryBezTRow[i]];
}
for (int i = 0; i < 8; ++i) {
point[interiorGregory[i]] = Bs[interiorBezSCol[i]] * Bt[interiorBezTRow[i]] * G[i];
}
//
// For derivatives, the basis functions for the interior points are rational and ideally
// require appropriate differentiation, i.e. product rule for the combination of B and G
// and the quotient rule for the rational G itself. As initially proposed by Loop et al
// though, the approximation using the 16 Bezier points arising from the G(s,t) has
// proved adequate (and is what the GPU shaders use) so we continue to use that here.
//
// An implementation of the true derivatives is provided and conditionally compiled for
// those that require it, e.g.:
//
// dclyde's note: skipping half of the product rule like this does seem to change the
// result a lot in my tests. This is not a runtime bottleneck for cloth sims anyway
// so I'm just using the accurate version.
//
if (wDs && wDt) {
bool find_second_partials = wDs && wDst && wDtt;
// Combined weights for boundary points -- simple tensor products:
for (int i = 0; i < 12; ++i) {
int iDst = boundaryGregory[i];
int tRow = boundaryBezTRow[i];
int sCol = boundaryBezSCol[i];
wDs[iDst] = Bds[sCol] * Bt[tRow];
wDt[iDst] = Bdt[tRow] * Bs[sCol];
if (find_second_partials) {
wDss[iDst] = Bdss[sCol] * Bt[tRow];
wDst[iDst] = Bds[sCol] * Bdt[tRow];
wDtt[iDst] = Bs[sCol] * Bdtt[tRow];
}
}
#ifndef OPENSUBDIV_GREGORY_EVAL_TRUE_DERIVATIVES
// Approximation to the true Gregory derivatives by differentiating the Bezier patch
// unique to the given (s,t), i.e. having F = (g^+ * f^+) + (g^- * f^-) as its four
// interior points:
//
// Combined weights for interior points -- tensor products with G+ or G-:
for (int i = 0; i < 8; ++i) {
int iDst = interiorGregory[i];
int tRow = interiorBezTRow[i];
int sCol = interiorBezSCol[i];
wDs[iDst] = Bds[sCol] * Bt[tRow] * G[i];
wDt[iDst] = Bdt[tRow] * Bs[sCol] * G[i];
if (find_second_partials) {
wDss[iDst] = Bdss[sCol] * Bt[tRow] * G[i];
wDst[iDst] = Bds[sCol] * Bdt[tRow] * G[i];
wDtt[iDst] = Bs[sCol] * Bdtt[tRow] * G[i];
}
}
#else
// True Gregory derivatives using appropriate differentiation of composite functions:
//
// Note that for G(s,t) = N(s,t) / D(s,t), all N' and D' are trivial constants (which
// simplifies things for higher order derivatives). And while each pair of functions
// G (i.e. the G+ and G- corresponding to points f+ and f-) must sum to 1 to ensure
// Bezier equivalence (when f+ = f-), the pairs of G' must similarly sum to 0. So we
// can potentially compute only one of the pair and negate the result for the other
// (and with 4 or 8 computations involving these constants, this is all very SIMD
// friendly...) but for now we treat all 8 independently for simplicity.
//
//REAL N[8] = { s, t, t, sC, sC, tC, tC, s };
REAL D[8] = { df0, df0, df1, df1, df2, df2, df3, df3 };
static REAL const Nds[8] = { 1.0f, 0.0f, 0.0f, -1.0f, -1.0f, 0.0f, 0.0f, 1.0f };
static REAL const Ndt[8] = { 0.0f, 1.0f, 1.0f, 0.0f, 0.0f, -1.0f, -1.0f, 0.0f };
static REAL const Dds[8] = { 1.0f, 1.0f, -1.0f, -1.0f, -1.0f, -1.0f, 1.0f, 1.0f };
static REAL const Ddt[8] = { 1.0f, 1.0f, 1.0f, 1.0f, -1.0f, -1.0f, -1.0f, -1.0f };
// Combined weights for interior points -- combinations of B, B', G and G':
for (int i = 0; i < 8; ++i) {
int iDst = interiorGregory[i];
int tRow = interiorBezTRow[i];
int sCol = interiorBezSCol[i];
// Quotient rule for G' (re-expressed in terms of G to simplify (and D = 1/D)):
REAL Gds = (Nds[i] - Dds[i] * G[i]) * D[i];
REAL Gdt = (Ndt[i] - Ddt[i] * G[i]) * D[i];
// Product rule combining B and B' with G and G':
wDs[iDst] = (Bds[sCol] * G[i] + Bs[sCol] * Gds) * Bt[tRow];
wDt[iDst] = (Bdt[tRow] * G[i] + Bt[tRow] * Gdt) * Bs[sCol];
if (find_second_partials) {
REAL Dsqr_inv = D[i]*D[i];
REAL Gdss = 2.0f * Dds[i] * Dsqr_inv * (G[i] * Dds[i] - Nds[i]);
REAL Gdst = Dsqr_inv * (2.0f * G[i] * Dds[i] * Ddt[i] - Nds[i] * Ddt[i] - Ndt[i] * Dds[i]);
REAL Gdtt = 2.0f * Ddt[i] * Dsqr_inv * (G[i] * Ddt[i] - Ndt[i]);
wDss[iDst] = (Bdss[sCol] * G[i] + 2.0f * Bds[sCol] * Gds + Bs[sCol] * Gdss) * Bt[tRow];
wDst[iDst] = Bt[tRow] * (Bs[sCol] * Gdst + Bds[sCol] * Gdt) +
Bdt[tRow] * (Bds[sCol] * G[i] + Bs[sCol] * Gds);
wDtt[iDst] = (Bdtt[tRow] * G[i] + 2.0f * Bdt[tRow] * Gdt + Bt[tRow] * Gdtt) * Bs[sCol];
}
}
#endif
}
return 20;
}
//
// Basis support for triangular patches:
//
// Triangular patches may be evaluated in barycentric (trivariate) or
// bivariate form, depending on the complexity of their basis functions.
// The parametric orientation for a triangle is as follows:
//
// (1,0)
// *
// . .
// t . 2 .
// . .
// . 0 1 .
// (0,0) *---------* (1,0)
// s
//
// With the origin (0,0) -- barycentric (0,0,w = 1) -- oriented at the
// corner V0, the corners V0, V1, and V2 correspond to barycentric
// coordinates W, U and V. This is consistent with GPU tessellation
// shaders, but not with many publications where the corners correspond
// more intuitively to U, V and W.
//
//
// Simple linear triangle:
//
template <typename REAL>
int
EvalBasisLinearTri(REAL s, REAL t,
REAL wP[3], REAL wDs[3], REAL wDt[3],
REAL wDss[3], REAL wDst[3], REAL wDtt[3]) {
if (wP) {
wP[0] = 1.0f - s - t;
wP[1] = s;
wP[2] = t;
}
if (wDs && wDt) {
wDs[0] = -1.0f;
wDs[1] = 1.0f;
wDs[2] = 0.0f;
wDt[0] = -1.0f;
wDt[1] = 0.0f;
wDt[2] = 1.0f;
if (wDss && wDst && wDtt) {
wDss[0] = wDss[1] = wDss[2] = 0.0f;
wDst[0] = wDst[1] = wDst[2] = 0.0f;
wDtt[0] = wDtt[1] = wDtt[2] = 0.0f;
}
}
return 3;
}
//
// Quartic Box spline triangle:
//
// Points for the quartic triangular Box spline (representing regular
// patches for Loop subdivision) are as follows:
//
// 10-----11
// . . . .
// . . . .
// 7-----8-----9
// . . . . . .
// . . . . . .
// 3-----4-----5-----6
// . . . . . .
// . . . . . /
// 0-----1-----2
//
// Stam provided the basis functions for these patches in terms of barycentric
// coordinates (u,v,w) (see Stam's "Evaluation of Loop Subdivision Surfaces").
// Unfortunately, unlike the basis functions for a quartic Bezier triangle,
// they are not very compact -- 3 functions involving 9 quartic terms and 3
// others involving 15 quartic terms. (In contrast, the maximum number of
// terms in bivariate form is 15.)
//
// Since we also need to differentiate with respect to u and v, we eliminate w
// and use the coefficient matrix C multiplied by the set of monomials M
// evaluated at (u,v), i.e. the full set of basis functions is:
//
// B(u,v) = C * M(u,v)
//
// where
//
// M(u,v) = { 1, u,v, uu,uv,vv, uuu,uuv,uvv,vvv, uuuu,uuuv,uuvv,uvvv,vvvv }
//
// and the 12 x 15 matrix C is as follows, scaled by a common factor of 1/12:
//
// { 1, -2,-4, 0, 6, 6, 2, 0, -6, -4, -1, -2, 0, 2, 1 },
// { 1, 2,-2, 0, -6, 0, -4, 0, 6, 2, 2, 4, 0, -2, -1 },
// { 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, -1, -2, 0, 0, 0 },
// { 1, -4,-2, 6, 6, 0, -4, -6, 0, 2, 1, 2, 0, -2, -1 },
// { 6, 0, 0, -12,-12,-12, 8, 12, 12, 8, -1, -2, 0, -2, -1 },
// { 1, 4, 2, 6, 6, 0, -4, -6,-12, -4, -1, -2, 0, 4, 2 },
// { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0 },
// { 1, -2, 2, 0, -6, 0, 2, 6, 0, -4, -1, -2, 0, 4, 2 },
// { 1, 2, 4, 0, 6, 6, -4,-12, -6, -4, 2, 4, 0, -2, -1 },
// { 0, 0, 0, 0, 0, 0, 2, 6, 6, 2, -1, -2, 0, -2, -1 },
// { 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, -2, -1 },
// { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1 }
//
// Differentiating the monomials and refactoring yields a unique set of
// coefficients for each of the derivatives, which we multiply by M(u,v).
//
namespace {
template <typename REAL>
inline void
evalBivariateMonomialsQuartic(REAL s, REAL t, REAL M[]) {
M[0] = 1.0;
M[1] = s;
M[2] = t;
M[3] = s * s;
M[4] = s * t;
M[5] = t * t;
M[6] = M[3] * s;
M[7] = M[4] * s;
M[8] = M[4] * t;
M[9] = M[5] * t;
M[10] = M[6] * s;
M[11] = M[7] * s;
M[12] = M[3] * M[5];
M[13] = M[8] * t;
M[14] = M[9] * t;
}
template <typename REAL>
void
evalBoxSplineTriDerivWeights(REAL const stMonomials[], int ds, int dt, REAL w[]) {
REAL const * M = stMonomials;
REAL S = 1.0;
int totalOrder = ds + dt;
if (totalOrder == 0) {
S *= (REAL) (1.0 / 12.0);
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]);
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]);
w[2] = S * ( 2*M[6] - M[10] - 2*M[11] );
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]);
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]);
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]);
w[6] = S * ( M[10] + 2*M[11] );
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]);
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]);
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]);
w[10] = S * ( 2*M[9] - 2*M[13] - M[14]);
w[11] = S * ( 2*M[13] + M[14]);
} else if (totalOrder == 1) {
S *= (REAL) (1.0 / 6.0);
if (ds) {
w[0] = S * (-1 + 3*M[2] + 3*M[3] - 3*M[5] - 2*M[6] - 3*M[7] + M[9]);
w[1] = S * ( 1 - 3*M[2] - 6*M[3] + 3*M[5] + 4*M[6] + 6*M[7] - M[9]);
w[2] = S * ( 3*M[3] - 2*M[6] - 3*M[7] );
w[3] = S * (-2 + 6*M[1] + 3*M[2] - 6*M[3] - 6*M[4] + 2*M[6] + 3*M[7] - M[9]);
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]);
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]);
w[6] = S * ( 2*M[6] + 3*M[7] );
w[7] = S * (-1 - 3*M[2] + 3*M[3] + 6*M[4] - 2*M[6] - 3*M[7] + 2*M[9]);
w[8] = S * ( 1 + 3*M[2] - 6*M[3] -12*M[4] - 3*M[5] + 4*M[6] + 6*M[7] - M[9]);
w[9] = S * ( 3*M[3] + 6*M[4] + 3*M[5] - 2*M[6] - 3*M[7] - M[9]);
w[10] = S * ( - M[9]);
w[11] = S * ( M[9]);
} else {
w[0] = S * (-2 + 3*M[1] + 6*M[2] - 6*M[4] - 6*M[5] - M[6] + 3*M[8] + 2*M[9]);
w[1] = S * (-1 - 3*M[1] + 6*M[4] + 3*M[5] + 2*M[6] - 3*M[8] - 2*M[9]);
w[2] = S * ( - M[6] );
w[3] = S * (-1 + 3*M[1] - 3*M[3] + 3*M[5] + M[6] - 3*M[8] - 2*M[9]);
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]);
w[5] = S * ( 1 + 3*M[1] - 3*M[3] -12*M[4] - 6*M[5] - M[6] + 6*M[8] + 4*M[9]);
w[6] = S * ( + M[6] );
w[7] = S * ( 1 - 3*M[1] + 3*M[3] - 6*M[5] - M[6] + 6*M[8] + 4*M[9]);
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]);
w[9] = S * ( + 3*M[3] + 6*M[4] + 3*M[5] - M[6] - 3*M[8] - 2*M[9]);
w[10] = S * ( 3*M[5] - 3*M[8] - 2*M[9]);
w[11] = S * ( 3*M[8] + 2*M[9]);
}
} else if (totalOrder == 2) {
if (ds == 2) {
w[0] = S * ( + M[1] - M[3] - M[4]);
w[1] = S * ( - 2*M[1] + 2*M[3] + 2*M[4]);
w[2] = S * ( M[1] - M[3] - M[4]);
w[3] = S * ( 1 - 2*M[1] - M[2] + M[3] + M[4]);
w[4] = S * (-2 + 4*M[1] + 2*M[2] - M[3] - M[4]);
w[5] = S * ( 1 - 2*M[1] - M[2] - M[3] - M[4]);
w[6] = S * ( M[3] + M[4]);
w[7] = S * ( + M[1] + M[2] - M[3] - M[4]);
w[8] = S * ( - 2*M[1] - 2*M[2] + 2*M[3] + 2*M[4]);
w[9] = S * ( M[1] + M[2] - M[3] - M[4]);
w[10] = 0;
w[11] = 0;
} else if (dt == 2) {
w[0] = S * ( 1 - M[1] - 2*M[2] + M[4] + M[5]);
w[1] = S * ( + M[1] + M[2] - M[4] - M[5]);
w[2] = 0;
w[3] = S * ( + M[2] - M[4] - M[5]);
w[4] = S * (-2 + 2*M[1] + 4*M[2] - M[4] - M[5]);
w[5] = S * ( - 2*M[1] - 2*M[2] + 2*M[4] + 2*M[5]);
w[6] = 0;
w[7] = S * ( - 2*M[2] + 2*M[4] + 2*M[5]);
w[8] = S * ( 1 - M[1] - 2*M[2] - M[4] - M[5]);
w[9] = S * ( + M[1] + M[2] - M[4] - M[5]);
w[10] = S * ( M[2] - M[4] - M[5]);
w[11] = S * ( M[4] + M[5]);
} else {
S *= (REAL) (1.0 / 2.0);
w[0] = S * ( 1 - 2*M[2] - M[3] + M[5]);
w[1] = S * (-1 + 2*M[2] + 2*M[3] - M[5]);
w[2] = S * ( - M[3] );
w[3] = S * ( 1 - 2*M[1] + M[3] - M[5]);
w[4] = S * (-2 + 4*M[1] + 4*M[2] - M[3] - M[5]);
w[5] = S * ( 1 - 2*M[1] - 4*M[2] - M[3] + 2*M[5]);
w[6] = S * ( + M[3] );
w[7] = S * (-1 + 2*M[1] - M[3] + 2*M[5]);
w[8] = S * ( 1 - 4*M[1] - 2*M[2] + 2*M[3] - M[5]);
w[9] = S * ( + 2*M[1] + 2*M[2] - M[3] - M[5]);
w[10] = S * ( - M[5]);
w[11] = S * ( M[5]);
}
} else {
assert(totalOrder <= 2);
}
}
template <typename REAL>
void
adjustBoxSplineTriBoundaryWeights(int boundaryMask, REAL weights[]) {
if (boundaryMask == 0) return;
//
// Determine boundary edges and vertices from the lower 3 and upper
// 2 bits of the 5-bit mask:
//
bool edgeIsBoundary[3 + 2]; // +2 filled in to avoid +1 and +2 mod 3
bool vertexIsBoundary[3];
bool lowerBits[3];
lowerBits[0] = (boundaryMask & 0x1) != 0;
lowerBits[1] = (boundaryMask & 0x2) != 0;
lowerBits[2] = (boundaryMask & 0x4) != 0;
int upperBits = (boundaryMask >> 3) & 0x3;
if (upperBits == 0) {
// Boundary edges only:
for (int i = 0; i < 3; ++i) {
edgeIsBoundary[i] = lowerBits[i];
vertexIsBoundary[i] = false;
}
} else if (upperBits == 1) {
// Boundary vertices only:
for (int i = 0; i < 3; ++i) {
vertexIsBoundary[i] = lowerBits[i];
edgeIsBoundary[i] = false;
}
} else if (upperBits == 2) {
// Boundary edge and opposite boundary vertex:
edgeIsBoundary[0] = vertexIsBoundary[2] = lowerBits[0];
edgeIsBoundary[1] = vertexIsBoundary[0] = lowerBits[1];
edgeIsBoundary[2] = vertexIsBoundary[1] = lowerBits[2];
}
// Wrap the 2 additional values to avoid modulo 3 in the edge tests:
edgeIsBoundary[3] = edgeIsBoundary[0];
edgeIsBoundary[4] = edgeIsBoundary[1];
//
// Adjust weights for the 4 boundary points (eB) and 3 interior points
// (eI) to account for the 3 phantom points (eP) adjacent to each
// boundary edge:
//
int const eP[3][3] = { { 0, 1, 2 }, { 6, 9, 11 }, { 10, 7, 3 } };
int const eB[3][4] = { {3, 4, 5, 6}, {2, 5, 8, 10 }, {11, 8, 4, 0} };
int const eI[3][3] = { { 7, 8, 9 }, { 1, 4, 7 }, { 9, 5, 1 } };
for (int i = 0; i < 3; ++i) {
if (edgeIsBoundary[i]) {
int const * iPhantom = eP[i];
int const * iBoundary = eB[i];
int const * iInterior = eI[i];
// Adjust weights for points contributing to phantom point
// P0 -- extrapolated according to the presence of adj edge:
REAL w0 = weights[iPhantom[0]];
if (edgeIsBoundary[i + 2]) {
// P0 = B1 + (B1 - I1)
weights[iBoundary[1]] += w0;
weights[iBoundary[1]] += w0;
weights[iInterior[1]] -= w0;
} else {
// P0 = B1 + (B0 - I0)
weights[iBoundary[1]] += w0;
weights[iBoundary[0]] += w0;
weights[iInterior[0]] -= w0;
}
// Adjust weights for points contributing to phantom point
// P1 = B1 + (B2 - I1)
REAL w1 = weights[iPhantom[1]];
weights[iBoundary[1]] += w1;
weights[iBoundary[2]] += w1;
weights[iInterior[1]] -= w1;
// Adjust weights for points contributing to phantom point
// P2 -- extrapolated according to the presence of adj edge:
REAL w2 = weights[iPhantom[2]];
if (edgeIsBoundary[i + 1]) {
// P2 = B2 + (B2 - I1)
weights[iBoundary[2]] += w2;
weights[iBoundary[2]] += w2;
weights[iInterior[1]] -= w2;
} else {
// P2 = B2 + (B3 - I2)
weights[iBoundary[2]] += w2;
weights[iBoundary[3]] += w2;
weights[iInterior[2]] -= w2;
}
// Clear weights for the phantom points:
weights[iPhantom[0]] = 0.0f;
weights[iPhantom[1]] = 0.0f;
weights[iPhantom[2]] = 0.0f;
}
}
//
// Adjust weights for the 3 boundary points (vB) and the 2 interior
// points (vI) to account for the 2 phantom points (vP) adjacent to
// each boundary vertex:
//
int const vP[3][2] = { { 3, 0 }, { 2, 6 }, { 11, 10 } };
int const vB[3][3] = { { 7, 4, 1 }, { 1, 5, 9 }, { 9, 8, 7 } };
int const vI[3][2] = { { 8, 5 }, { 4, 8 }, { 5, 4 } };
for (int i = 0; i < 3; ++i) {
if (vertexIsBoundary[i]) {
int const * iPhantom = vP[i];
int const * iBoundary = vB[i];
int const * iInterior = vI[i];
// Adjust weights for points contributing to phantom point
// P0 = B1 + (B0 - I0)
REAL w0 = weights[iPhantom[0]];
weights[iBoundary[1]] += w0;
weights[iBoundary[0]] += w0;
weights[iInterior[0]] -= w0;
// Adjust weights for points contributing to phantom point
// P1 = B1 + (B2 - I1)
REAL w1 = weights[iPhantom[1]];
weights[iBoundary[1]] += w1;
weights[iBoundary[2]] += w1;
weights[iInterior[1]] -= w1;
// Clear weights for the phantom points:
weights[iPhantom[0]] = 0.0f;
weights[iPhantom[1]] = 0.0f;
}
}
}
template <typename REAL>
void
boundBasisBoxSplineTri(int boundary,
REAL wP[12], REAL wDs[12], REAL wDt[12],
REAL wDss[12], REAL wDst[12], REAL wDtt[12]) {
if (wP) {
adjustBoxSplineTriBoundaryWeights(boundary, wP);
}
if (wDs && wDt) {
adjustBoxSplineTriBoundaryWeights(boundary, wDs);
adjustBoxSplineTriBoundaryWeights(boundary, wDt);
if (wDss && wDst && wDtt) {
adjustBoxSplineTriBoundaryWeights(boundary, wDss);
adjustBoxSplineTriBoundaryWeights(boundary, wDst);
adjustBoxSplineTriBoundaryWeights(boundary, wDtt);
}
}
}
} // namespace
template <typename REAL>
int EvalBasisBoxSplineTri(REAL s, REAL t,
REAL wP[12], REAL wDs[12], REAL wDt[12],
REAL wDss[12], REAL wDst[12], REAL wDtt[12]) {
REAL stMonomials[15];
evalBivariateMonomialsQuartic(s, t, stMonomials);
if (wP) {
evalBoxSplineTriDerivWeights<REAL>(stMonomials, 0, 0, wP);
}
if (wDs && wDt) {
evalBoxSplineTriDerivWeights(stMonomials, 1, 0, wDs);
evalBoxSplineTriDerivWeights(stMonomials, 0, 1, wDt);
if (wDss && wDst && wDtt) {
evalBoxSplineTriDerivWeights(stMonomials, 2, 0, wDss);
evalBoxSplineTriDerivWeights(stMonomials, 1, 1, wDst);
evalBoxSplineTriDerivWeights(stMonomials, 0, 2, wDtt);
}
}
return 12;
}
//
// Hybrid (cubic-quartic) Bezier triangle:
//
// The regular patch for Loop subdivision is a quartic triangular Box spline
// with cubic boundaries. So we need a quartic Bezier patch to represent it
// faithfully, but we use a cubic-quartic hybrid to keep the representation
// of boundaries as cubic -- useful for a number of purposes, in addition to
// reducing the number points required from 15 to 12.
//
// Ultimately this patch is quartic and its basis functions are of maximum
// quartic degree. The formulae for the 15 true quartic basis functions are:
//
// 4! i j k
// B (u,v,w) = ------- * (u * v * w )
// ijk i!j!k!
//
// for each i + j + k = 4, and the quartic points and corresponding p<i,j,k>
// are oriented as follows:
//
// Q14 p040
// Q12 Q13 p031 p130
// Q9 Q10 Q11 p022 p121 p220
// Q5 Q6 Q7 Q8 p013 p112 p211 p310
// Q0 Q1 Q2 Q3 Q4 p004 p103 p202 p301 p400
//
// The points for the corresponding hybrid patch are oriented and numbered:
//
// H11
// H8 H10
// H9
// H4 H5 H6 H7
// H0 H1 H2 H3
//
// Their corresponding basis functions h(u,v,w) are derived by combining the
// quartic basis functions according to degree elevation of their boundary
// curves. This leads to the 12 basis functions:
//
// h[0] = w^3
// h[3] = u^3
// h[11] = v^3
//
// h[1] = 3 * u * w^2 * (u + w)
// h[2] = 3 * u^2 * w * (u + w)
//
// h[7] = 3 * u^2 * v * (u + v)
// h[10] = 3 * u * v^2 * (u + v)
//
// h[8] = 3 * v^2 * w * (w + v)
// h[10] = 3 * v * w^2 * (v + w)
//
// h[5] = 12 * u * v * w^2;
// h[6] = 12 * u^2 * v * w;
// h[9] = 12 * u * v^2 * w;
//
// These remain compact with at most two trivariate terms, and so relatively
// easy to differentiate in this form while keeping the number of terms low.
//
namespace {
template <typename REAL>
void
evalBezierTriDerivWeights(REAL s, REAL t, int ds, int dt, REAL wB[]) {
REAL u = s;
REAL v = t;
REAL w = 1 - u - v;
REAL u2 = u * u;
REAL v2 = v * v;
REAL w2 = w * w;
REAL uv = u * v;
REAL vw = v * w;
REAL uw = u * w;
int totalOrder = ds + dt;
if (totalOrder == 0) {
wB[0] = w*w2;
wB[3] = u*u2;
wB[11] = v*v2;
wB[1] = 3 * uw * (uw + w2);
wB[2] = 3 * uw * (uw + u2);
wB[7] = 3 * uv * (uv + u2);
wB[10] = 3 * uv * (uv + v2);
wB[8] = 3 * vw * (vw + v2);
wB[4] = 3 * vw * (vw + w2);
wB[5] = 12 * w2 * uv;
wB[6] = 12 * u2 * vw;
wB[9] = 12 * v2 * uw;
} else if (totalOrder == 1) {
if (ds) {
wB[0] = -3 * w2;
wB[3] = 3 * u2;
wB[11] = 0;
wB[1] = 3 * w * (w2 - uw - 2*u2);
wB[2] = -3 * u * (u2 - uw - 2*w2);
wB[7] = 9 * u2*v + 6 * u*v2;
wB[10] = 3 * v*v2 + 6 * u*v2;
wB[8] = -3 * v*v2 - 6 * v2*w;
wB[4] = -9 * v*w2 - 6 * v2*w;
wB[5] = 12 * vw * (w - 2*u);
wB[6] = 12 * uv * (2*w - u);
wB[9] = 12 * v2 * (w - u);
} else {
wB[0] = -3 * w2;
wB[3] = 0;
wB[11] = 3 * v2;
wB[1] = -9 * u*w2 - 6 * u2*w;
wB[2] = -3 * u*u2 - 6 * u2*w;
wB[7] = 3 * u*u2 + 6 * u2*v;
wB[10] = 9 * u*v2 + 6 * u2*v;
wB[8] = -3 * v * (v2 - vw - 2*w2);
wB[4] = 3 * w * (w2 - vw - 2*v2);
wB[5] = 12 * uw * (w - 2*v);
wB[6] = 12 * u2 * (w - v);
wB[9] = 12 * uv * (2*w - v);
}
} else if (totalOrder == 2) {
if (ds == 2) {
wB[0] = 6 * w;
wB[3] = 6 * u;
wB[11] = 0;
wB[1] = 6 * (u2 - uw - 2*w2);
wB[2] = 6 * (w2 - uw - 2*u2);
wB[7] = 6 * v2 + 18 * uv;
wB[10] = 6 * v2;
wB[8] = 6 * v2;
wB[4] = 6 * v2 + 18 * vw;
wB[5] = 24 * (uv - 2*vw);
wB[6] = 24 * (vw - 2*uv);
wB[9] = -24 * v2;
} else if (dt == 2) {
wB[0] = 6 * w;
wB[3] = 0;
wB[11] = 6 * v;
wB[1] = 6 * u2 + 18 * uw;
wB[2] = 6 * u2;
wB[7] = 6 * u2;
wB[10] = 6 * u2 + 18 * uv;
wB[8] = 6 * (w2 - vw - 2*v2);
wB[4] = 6 * (v2 - vw - 2*w2);
wB[5] = 24 * (uv - 2*uw);
wB[6] = -24 * u2;
wB[9] = 24 * (uw - 2*uv);
} else {
wB[0] = 6 * w;
wB[3] = 0;
wB[11] = 0;
wB[1] = 6 * (u2 + uw - 1.5f*w2);
wB[2] = -3 * (u2 + 4*uw);
wB[7] = 9 * u2 + 12 * uv;
wB[10] = 9 * v2 + 12 * uv;
wB[8] = -3 * (v2 + 4*vw);
wB[4] = 6 * (v2 + vw - 1.5f*w2);
wB[5] = 24 * (uv - vw - uw + 0.5f*w2);
wB[6] = -24 * (uv - uw + 0.5f*u2);
wB[9] = -24 * (uv - vw + 0.5f*v2);
}
} else {
assert(totalOrder <= 2);
}
}
} // end namespace
template <typename REAL>
int
EvalBasisBezierTri(REAL s, REAL t,
REAL wP[12], REAL wDs[12], REAL wDt[12],
REAL wDss[12], REAL wDst[12], REAL wDtt[12]) {
if (wP) {
evalBezierTriDerivWeights<REAL>(s, t, 0, 0, wP);
}
if (wDs && wDt) {
evalBezierTriDerivWeights(s, t, 1, 0, wDs);
evalBezierTriDerivWeights(s, t, 0, 1, wDt);
if (wDss && wDst && wDtt) {
evalBezierTriDerivWeights(s, t, 2, 0, wDss);
evalBezierTriDerivWeights(s, t, 1, 1, wDst);
evalBezierTriDerivWeights(s, t, 0, 2, wDtt);
}
}
return 12;
}
//
// Hybrid (cubic-quartic) Gregory triangle:
//
// As with the Bezier triangle, and consistent with Loop, Schaefer at al (in
// ("Approximating Subdivision Surfaces with Gregory Patches for Hardware
// Tessellation") we use a cubic-quartic hybrid Gregory patch. Like the
// quad Gregory patch, this patch uses Bezier basis functions (from the
// cubic-quartic hybrid above) and rational multipliers to blend pairs of
// interior points (face points).
//
namespace {
//
// Expanding a set of 12 Bezier basis functions for the 6 (3 pairs) of
// rational weights for the 15 Gregory basis functions:
//
template <typename REAL>
void
convertBezierWeightsToGregory(REAL const wB[12], REAL const rG[6], REAL wG[15]) {
wG[0] = wB[0];
wG[1] = wB[1];
wG[2] = wB[4];
wG[3] = wB[5] * rG[0];
wG[4] = wB[5] * rG[1];
wG[5] = wB[3];
wG[6] = wB[7];
wG[7] = wB[2];
wG[8] = wB[6] * rG[2];
wG[9] = wB[6] * rG[3];
wG[10] = wB[11];
wG[11] = wB[8];
wG[12] = wB[10];
wG[13] = wB[9] * rG[4];
wG[14] = wB[9] * rG[5];
}
} // end namespace
template <typename REAL>
int
EvalBasisGregoryTri(REAL s, REAL t,
REAL wP[15], REAL wDs[15], REAL wDt[15],
REAL wDss[15], REAL wDst[15], REAL wDtt[15]) {
//
// Bezier basis functions are denoted with B while the rational multipliers for the
// interior points will be denoted G -- so we have B(s,t) and G(s,t) (though we
// switch to barycentric (u,v,w) briefly to compute G)
//
REAL BP[12], BDs[12], BDt[12], BDss[12], BDst[12], BDtt[12];
REAL G[6] = { 1.0f, 0.0f, 1.0f, 0.0f, 1.0f, 0.0f };
REAL u = s;
REAL v = t;
REAL w = 1 - u - v;
if ((u + v) > 0) {
G[0] = u / (u + v);
G[1] = v / (u + v);
}
if ((v + w) > 0) {
G[2] = v / (v + w);
G[3] = w / (v + w);
}
if ((w + u) > 0) {
G[4] = w / (w + u);
G[5] = u / (w + u);
}
//
// Compute Bezier basis functions and convert, adjusting interior points:
//
if (wP) {
evalBezierTriDerivWeights<REAL>(s, t, 0, 0, BP);
convertBezierWeightsToGregory(BP, G, wP);
}
if (wDs && wDt) {
// TBD -- ifdef OPENSUBDIV_GREGORY_EVAL_TRUE_DERIVATIVES
evalBezierTriDerivWeights(s, t, 1, 0, BDs);
evalBezierTriDerivWeights(s, t, 0, 1, BDt);
convertBezierWeightsToGregory(BDs, G, wDs);
convertBezierWeightsToGregory(BDt, G, wDt);
if (wDss && wDst && wDtt) {
evalBezierTriDerivWeights(s, t, 2, 0, BDss);
evalBezierTriDerivWeights(s, t, 1, 1, BDst);
evalBezierTriDerivWeights(s, t, 0, 2, BDtt);
convertBezierWeightsToGregory(BDss, G, wDss);
convertBezierWeightsToGregory(BDst, G, wDst);
convertBezierWeightsToGregory(BDtt, G, wDtt);
}
}
return 15;
}
//
// Higher level basis evaluation functions that deal with parameterization and
// boundary issues (reflected in PatchParam) for all patch types:
//
template <typename REAL>
int
EvaluatePatchBasisNormalized(int patchType, PatchParam const & param, REAL s, REAL t,
REAL wP[], REAL wDs[], REAL wDt[],
REAL wDss[], REAL wDst[], REAL wDtt[]) {
int boundaryMask = param.GetBoundary();
int nPoints = 0;
if (patchType == PatchDescriptor::REGULAR) {
nPoints = EvalBasisBSpline(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
if (boundaryMask) {
boundBasisBSpline(boundaryMask, wP, wDs, wDt, wDss, wDst, wDtt);
}
} else if (patchType == PatchDescriptor::LOOP) {
nPoints = EvalBasisBoxSplineTri(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
if (boundaryMask) {
boundBasisBoxSplineTri(boundaryMask, wP, wDs, wDt, wDss, wDst, wDtt);
}
} else if (patchType == PatchDescriptor::GREGORY_BASIS) {
nPoints = EvalBasisGregory(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
} else if (patchType == PatchDescriptor::GREGORY_TRIANGLE) {
nPoints = EvalBasisGregoryTri(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
} else if (patchType == PatchDescriptor::QUADS) {
nPoints = EvalBasisLinear(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
} else if (patchType == PatchDescriptor::TRIANGLES) {
nPoints = EvalBasisLinearTri(s, t, wP, wDs, wDt, wDss, wDst, wDtt);
} else {
assert(0);
}
return nPoints;
}
template <typename REAL>
int
EvaluatePatchBasis(int patchType, PatchParam const & param, REAL s, REAL t,
REAL wP[], REAL wDs[], REAL wDt[],
REAL wDss[], REAL wDst[], REAL wDtt[]) {
REAL derivSign = 1.0f;
if ((patchType == PatchDescriptor::LOOP) ||
(patchType == PatchDescriptor::GREGORY_TRIANGLE) ||
(patchType == PatchDescriptor::TRIANGLES)) {
param.NormalizeTriangle(s, t);
if (param.IsTriangleRotated()) {
derivSign = -1.0f;
}
} else {
param.Normalize(s, t);
}
int nPoints = EvaluatePatchBasisNormalized(
patchType, param, s, t, wP, wDs, wDt, wDss, wDst, wDtt);
if (wDs && wDt) {
REAL d1Scale = derivSign * (REAL)(1 << param.GetDepth());
for (int i = 0; i < nPoints; ++i) {
wDs[i] *= d1Scale;
wDt[i] *= d1Scale;
}
if (wDss && wDst && wDtt) {
REAL d2Scale = derivSign * d1Scale * d1Scale;
for (int i = 0; i < nPoints; ++i) {
wDss[i] *= d2Scale;
wDst[i] *= d2Scale;
wDtt[i] *= d2Scale;
}
}
}
return nPoints;
}
//
// Explicit float and double instantiations:
//
template int EvaluatePatchBasisNormalized<float>(int patchType, PatchParam const & param,
float s, float t, float wP[], float wDs[], float wDt[], float wDss[], float wDst[], float wDtt[]);
template int EvaluatePatchBasis<float>(int patchType, PatchParam const & param,
float s, float t, float wP[], float wDs[], float wDt[], float wDss[], float wDst[], float wDtt[]);
template int EvaluatePatchBasisNormalized<double>(int patchType, PatchParam const & param,
double s, double t, double wP[], double wDs[], double wDt[], double wDss[], double wDst[], double wDtt[]);
template int EvaluatePatchBasis<double>(int patchType, PatchParam const & param,
double s, double t, double wP[], double wDs[], double wDt[], double wDss[], double wDst[], double wDtt[]);
//
// Most basis evaluation functions are implicitly instantiated above -- Bezier
// require explicit instantiation as they are not invoked via a patch type:
//
template int EvalBasisBezier<float>(float s, float t,
float wP[16], float wDs[16], float wDt[16], float wDss[16], float wDst[16], float wDtt[16]);
template int EvalBasisBezierTri<float>(float s, float t,
float wP[12], float wDs[12], float wDt[12], float wDss[12], float wDst[12], float wDtt[12]);
template int EvalBasisBezier<double>(double s, double t,
double wP[16], double wDs[16], double wDt[16], double wDss[16], double wDst[16], double wDtt[16]);
template int EvalBasisBezierTri<double>(double s, double t,
double wP[12], double wDs[12], double wDt[12], double wDss[12], double wDst[12], double wDtt[12]);
} // end namespace internal
} // end namespace Far
} // end namespace OPENSUBDIV_VERSION
} // end namespace OpenSubdiv