AuroraMiddleware/OpenSubdiv
1
0
Fork 0
mirror of https://github.com/PixarAnimationStudios/OpenSubdiv synced 2024-11-27 14:00:10 +00:00
Code Issues Releases Wiki Activity

Merge pull request #977 from barfowl/real_templates_internal

Browse Source 
Internal improvements to support double precision
This commit is contained in:
David G Yu 2018-07-24 13:42:05 -07:00 committed by GitHub
commit b3ab407553
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
13 changed files with 596 additions and 464 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
opensubdiv/far/bilinearPatchBuilder.cpp
View File

@ -73,10 +73,11 @@ BilinearPatchBuilder::patchTypeFromBasis(BasisType basis) const {
return patchTypeFromBasisArray[(int)basis];
}
template <typename REAL>
int
BilinearPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
BilinearPatchBuilder::convertSourcePatch(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const {
SparseMatrix<REAL> & matrix) const {
assert("Conversion from Bilinear patches to other bases not yet supported" == 0);
@ -87,6 +88,19 @@ BilinearPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
return -1;
}
int
BilinearPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const {
return convertSourcePatch(sourcePatch, patchType, matrix);
}
int
BilinearPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<double> & matrix) const {
return convertSourcePatch(sourcePatch, patchType, matrix);
}
} // end namespace Far
} // end namespace OPENSUBDIV_VERSION

8
opensubdiv/far/bilinearPatchBuilder.h
View File

@ -53,6 +53,14 @@ protected:
virtual int convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const;
virtual int convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<double> & matrix) const;
private:
template <typename REAL>
int convertSourcePatch(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<REAL> & matrix) const;
};
} // end namespace Far

283
opensubdiv/far/catmarkPatchBuilder.cpp
View File

@ -22,12 +22,12 @@
// language governing permissions and limitations under the Apache License.
//
#include "../far/catmarkPatchBuilder.h"
#include "../vtr/stackBuffer.h"
#include <cassert>
#include <cmath>
#include <cstdio>
#include <cassert>
#include "../far/catmarkPatchBuilder.h"
#include "../vtr/stackBuffer.h"
namespace OpenSubdiv {
namespace OPENSUBDIV_VERSION {
@ -53,51 +53,63 @@ namespace Far {
// retrieved from Sdc::Scheme to ensure they conform, so future factoring
// of the formulae is still necessary.
//
// XXXX (barfowl) - /internal computations are not implemented in terms
// of <REAL>, only the interface supports <REAL>. Need to eventually
// remove remaining computations using "float" and define a mechanism
// to invoke the appropriate precision math functions, e.g. calling
// something like Math<float>::Cos(), etc.
// Regarding support for multiple precision, like Sdc, some intermediate
// calculations are performed in double and cast to float. Historically
// this conversion code has been purely float and later extended to
// support template <typename REAL>. For math functions such as cos(),
// sin(), etc., we rely on overloading via <cmath> through the use of
// std::cos(), std::sin(), etc.
//
template <typename REAL>
struct CatmarkLimits {
public:
typedef REAL Weight;
static void ComputeInteriorPointWeights(int valence, int faceInRing,
REAL* pWeights, REAL* epWeights, REAL* emWeights);
Weight* pWeights, Weight* epWeights, Weight* emWeights);
static void ComputeBoundaryPointWeights(int valence, int faceInRing,
REAL* pWeights, REAL* epWeights, REAL* emWeights);
Weight* pWeights, Weight* epWeights, Weight* emWeights);
private:
//
// Lookup table and formula for the scale factor applied to limit
// tangents that arises from eigen values of the subdivision matrix
//
static inline float computeCoefficient(int valence) {
// precomputed coefficient table up to valence 29
static float efTable[] = {
0, 0, 0,
0.812816f, 0.500000f, 0.363644f, 0.287514f,
0.238688f, 0.204544f, 0.179229f, 0.159657f,
0.144042f, 0.131276f, 0.120632f, 0.111614f,
0.103872f, 0.09715f, 0.0912559f, 0.0860444f,
0.0814022f, 0.0772401f, 0.0734867f, 0.0700842f,
0.0669851f, 0.0641504f, 0.0615475f, 0.0591488f,
0.0569311f, 0.0548745f, 0.0529621f
};
assert(valence > 0);
if (valence < 30) return efTable[valence];
float t = 2.0f * float(M_PI) / float(valence);
return 1.0f / (valence * (cosf(t) + 5.0f +
sqrtf((cosf(t) + 9) * (cosf(t) + 1)))/16.0f);
}
static double computeCoefficient(int valence);
};
//
// Lookup table and formula for the scale factor applied to limit
// tangents that arises from eigen values of the subdivision matrix.
// Historically 30 values have been stored -- up to valence 29.
//
template <typename REAL>
inline double
CatmarkLimits<REAL>::computeCoefficient(int valence) {
static double const efTable[] = {
0.0, 0.0, 0.0,
8.1281572906372312e-01, 0.5, 3.6364406329142801e-01,
2.8751379706077085e-01, 2.3868786685851678e-01, 2.0454364190756097e-01,
1.7922903958061159e-01, 1.5965737079986253e-01, 1.4404233443011302e-01,
1.3127568415883017e-01, 1.2063172212675841e-01, 1.1161437506676930e-01,
1.0387245516114274e-01, 9.7150019090724835e-02, 9.1255917505950648e-02,
8.6044378511602668e-02, 8.1402211336798411e-02, 7.7240129516184072e-02,
7.3486719751997026e-02, 7.0084157479797987e-02, 6.6985104030725440e-02,
6.4150420569810074e-02, 6.1547457638637268e-02, 5.9148757447233989e-02,
5.6931056818776957e-02, 5.4874512279256417e-02, 5.2962091433796134e-02
};
assert(valence > 0);
if (valence < 30) return efTable[valence];
double invValence = 1.0 / valence;
double cosT = std::cos(2.0 * M_PI * invValence);
double divisor = (cosT + 5.0) + std::sqrt((cosT + 9.0) * (cosT + 1.0));
return (16.0 * invValence / divisor);
}
template <typename REAL>
void
CatmarkLimits<REAL>::ComputeInteriorPointWeights(int valence, int faceInRing,
REAL* pWeights, REAL* epWeights, REAL* emWeights) {
Weight* pWeights, Weight* epWeights, Weight* emWeights) {
//
// For the limit tangents of an interior vertex, the second tangent is a
@ -123,14 +135,14 @@ CatmarkLimits<REAL>::ComputeInteriorPointWeights(int valence, int faceInRing,
//
bool computeEdgePoints = epWeights && emWeights;
float fValence = float(valence);
float oneOverValence = 1.0f / fValence;
float oneOverValPlus5 = 1.0f / (fValence + 5.0f);
double fValence = (double) valence;
double oneOverValence = 1.0f / fValence;
double oneOverValPlus5 = 1.0f / (fValence + 5.0f);
float pCoeff = oneOverValence * oneOverValPlus5;
float tanCoeff = computeCoefficient(valence) * 0.5f * oneOverValPlus5;
double pCoeff = oneOverValence * oneOverValPlus5;
double tanCoeff = computeCoefficient(valence) * 0.5f * oneOverValPlus5;
float faceAngle = 2.0f * float(M_PI) * oneOverValence;
double faceAngle = 2.0 * M_PI * oneOverValence;
//
// Assign position weights directly while accumulating an intermediate set
@ -139,28 +151,28 @@ CatmarkLimits<REAL>::ComputeInteriorPointWeights(int valence, int faceInRing,
// have to deal with the off-by-one offset within the loop:
//
int weightWidth = 1 + 2 * valence;
StackBuffer<REAL, 64, true> tanWeights(weightWidth);
std::memset(&tanWeights[0], 0, weightWidth * sizeof(REAL));
StackBuffer<Weight, 64, true> tanWeights(weightWidth);
std::memset(&tanWeights[0], 0, weightWidth * sizeof(Weight));
pWeights[0] = fValence * oneOverValPlus5;
pWeights[0] = (Weight) (fValence * oneOverValPlus5);
REAL *pW = &pWeights[1];
REAL *tW = &tanWeights[1];
Weight *pW = &pWeights[1];
Weight *tW = &tanWeights[1];
for (int i = 0; i < valence; ++i) {
pW[2*i] = pCoeff * 4.0f;
pW[2*i + 1] = pCoeff;
pW[2*i] = (Weight) (pCoeff * 4.0f);
pW[2*i + 1] = (Weight) pCoeff;
if (computeEdgePoints) {
int iPrev = (i + valence - 1) % valence;
int iNext = (i + 1) % valence;
float cosICoeff = tanCoeff * cosf(float(i) * faceAngle);
double cosICoeff = tanCoeff * std::cos(faceAngle * i);
tW[2*iPrev] += cosICoeff * 2.0f;
tW[2*iPrev + 1] += cosICoeff;
tW[2*i] += cosICoeff * 4.0f;
tW[2*i + 1] += cosICoeff;
tW[2*iNext] += cosICoeff * 2.0f;
tW[2*iPrev] += (Weight) (cosICoeff * 2.0f);
tW[2*iPrev + 1] += (Weight) cosICoeff;
tW[2*i] += (Weight) (cosICoeff * 4.0f);
tW[2*i + 1] += (Weight) cosICoeff;
tW[2*iNext] += (Weight) (cosICoeff * 2.0f);
}
}
@ -189,10 +201,10 @@ CatmarkLimits<REAL>::ComputeInteriorPointWeights(int valence, int faceInRing,
template <typename REAL>
void
CatmarkLimits<REAL>::ComputeBoundaryPointWeights(int valence, int faceInRing,
REAL* pWeights, REAL* epWeights, REAL* emWeights) {
Weight* pWeights, Weight* epWeights, Weight* emWeights) {
int numFaces = valence - 1;
float faceAngle = float(M_PI) / float(numFaces);
int numFaces = valence - 1;
double faceAngle = M_PI / numFaces;
int weightWidth = 2 * valence;
@ -201,11 +213,11 @@ CatmarkLimits<REAL>::ComputeBoundaryPointWeights(int valence, int faceInRing,
//
// Position weights are trivial:
//
std::memset(&pWeights[0], 0, weightWidth * sizeof(float));
std::memset(&pWeights[0], 0, weightWidth * sizeof(Weight));
pWeights[0] = 4.0f / 6.0f;
pWeights[1] = 1.0f / 6.0f;
pWeights[N] = 1.0f / 6.0f;
pWeights[0] = (Weight) (4.0 / 6.0);
pWeights[1] = (Weight) (1.0 / 6.0);
pWeights[N] = (Weight) (1.0 / 6.0);
if ((epWeights == 0) && (emWeights == 0)) return;
@ -215,32 +227,36 @@ CatmarkLimits<REAL>::ComputeBoundaryPointWeights(int valence, int faceInRing,
// by two non-zero weights, so allocate and compute weights for the
// interior tangent:
//
float tBoundaryWeight_1 = 1.0f / 6.0f;
float tBoundaryWeight_N = -1.0f / 6.0f;
double tBoundaryCoeff_1 = ( 1.0 / 6.0);
double tBoundaryCoeff_N = (-1.0 / 6.0);
StackBuffer<REAL, 64, true> tanWeights(weightWidth);
StackBuffer<Weight, 64, true> tanWeights(weightWidth);
{
float k = float(numFaces);
float theta = faceAngle;
float c = cosf(theta);
float s = sinf(theta);
float div3 = 1.0f / 3.0f;
float div3kc = 1.0f / (3.0f*k+c);
float gamma = -4.0f * s * div3kc;
float alpha_0k = -((1.0f+2.0f*c) * sqrtf(1.0f+c)) * div3kc / sqrtf(1.0f-c);
float beta_0 = s * div3kc;
double k = (double) numFaces;
double theta = faceAngle;
double c = std::cos(theta);
double s = std::sin(theta);
double div3 = 1.0 / 3.0;
double div3kc = 1.0f / (3.0f * k + c);
double gamma = -4.0f * s * div3kc;
double alpha_0k = -((1.0f + 2.0f * c) * std::sqrt(1.0f + c)) * div3kc
/ std::sqrt(1.0f - c);
double beta_0 = s * div3kc;
tanWeights[0] = gamma * div3;
tanWeights[1] = alpha_0k * div3;
tanWeights[2] = beta_0 * div3;
tanWeights[N] = alpha_0k * div3;
tanWeights[0] = (Weight) (gamma * div3);
tanWeights[1] = (Weight) (alpha_0k * div3);
tanWeights[2] = (Weight) (beta_0 * div3);
tanWeights[N] = (Weight) (alpha_0k * div3);
for (int i = 1; i < valence - 1; ++i) {
float alpha = 4.0f * sinf(float(i)*theta) * div3kc;
float beta = (sinf(float(i)*theta) + sinf(float(i+1)*theta)) * div3kc;
double sinThetaI = std::sin(theta * i);
double sinThetaIplus1 = std::sin(theta * (i+1));
tanWeights[1 + 2*i] = alpha * div3;
tanWeights[1 + 2*i + 1] = beta * div3;
double alpha = 4.0f * sinThetaI * div3kc;
double beta = (sinThetaI + sinThetaIplus1) * div3kc;
tanWeights[1 + 2*i] = (Weight) (alpha * div3);
tanWeights[1 + 2*i + 1] = (Weight) (beta * div3);
}
}
@ -249,23 +265,23 @@ CatmarkLimits<REAL>::ComputeBoundaryPointWeights(int valence, int faceInRing,
//
if (faceInRing == 0) {
// Ep is on boundary edge and has only two weights: w[1] and w[N]
std::memset(&epWeights[0], 0, weightWidth * sizeof(float));
std::memset(&epWeights[0], 0, weightWidth * sizeof(Weight));
epWeights[0] = 2.0f / 3.0f;
epWeights[1] = 1.0f / 3.0f;
epWeights[0] = (Weight) (2.0 / 3.0);
epWeights[1] = (Weight) (1.0 / 3.0);
} else {
// Ep is on interior edge and has all weights
int iEdgeNext = faceInRing;
float faceAngleNext = faceAngle * float(iEdgeNext);
float cosAngleNext = cosf(faceAngleNext);
float sinAngleNext = sinf(faceAngleNext);
double faceAngleNext = faceAngle * iEdgeNext;
double cosAngleNext = std::cos(faceAngleNext);
double sinAngleNext = std::sin(faceAngleNext);
for (int i = 0; i < weightWidth; ++i) {
epWeights[i] = tanWeights[i] * sinAngleNext;
epWeights[i] = (Weight)(tanWeights[i] * sinAngleNext);
}
epWeights[0] += pWeights[0];
epWeights[1] += pWeights[1] + tBoundaryWeight_1 * cosAngleNext;
epWeights[N] += pWeights[N] + tBoundaryWeight_N * cosAngleNext;
epWeights[1] += pWeights[1] + (Weight)(tBoundaryCoeff_1 * cosAngleNext);
epWeights[N] += pWeights[N] + (Weight)(tBoundaryCoeff_N * cosAngleNext);
}
//
@ -273,23 +289,23 @@ CatmarkLimits<REAL>::ComputeBoundaryPointWeights(int valence, int faceInRing,
//
if (faceInRing == (numFaces - 1)) {
// Em is on boundary edge and has only two weights: w[1] and w[N]
std::memset(&emWeights[0], 0, weightWidth * sizeof(float));
std::memset(&emWeights[0], 0, weightWidth * sizeof(Weight));
emWeights[0] = 2.0f / 3.0f;
emWeights[N] = 1.0f / 3.0f;
emWeights[0] = (Weight) (2.0 / 3.0);
emWeights[N] = (Weight) (1.0 / 3.0);
} else {
// Em is on interior edge and has all weights
int iEdgePrev = (faceInRing + 1) % valence;
float faceAnglePrev = faceAngle * float(iEdgePrev);
float cosAnglePrev = cosf(faceAnglePrev);
float sinAnglePrev = sinf(faceAnglePrev);
double faceAnglePrev = faceAngle * iEdgePrev;
double cosAnglePrev = std::cos(faceAnglePrev);
double sinAnglePrev = std::sin(faceAnglePrev);
for (int i = 0; i < weightWidth; ++i) {
emWeights[i] = tanWeights[i] * sinAnglePrev;
emWeights[i] = (Weight)(tanWeights[i] * sinAnglePrev);
}
emWeights[0] += pWeights[0];
emWeights[1] += pWeights[1] + tBoundaryWeight_1 * cosAnglePrev;
emWeights[N] += pWeights[N] + tBoundaryWeight_N * cosAnglePrev;
emWeights[1] += pWeights[1] + (Weight)(tBoundaryCoeff_1 * cosAnglePrev);
emWeights[N] += pWeights[N] + (Weight)(tBoundaryCoeff_N * cosAnglePrev);
}
}
@ -711,12 +727,11 @@ GregoryConverter<REAL>::Initialize(SourcePatch const & sourcePatch) {
corner.cosFaceAngle = 0.0f;
corner.sinFaceAngle = 1.0f;
} else {
// XXXX (barfowl) - use of sine/cosine here needs to respect <REAL>
corner.faceAngle =
(corner.isBoundary ? REAL(M_PI) : (2.0f * REAL(M_PI)))
/ REAL(corner.numFaces);
corner.cosFaceAngle = cosf(corner.faceAngle);
corner.sinFaceAngle = sinf(corner.faceAngle);
corner.cosFaceAngle = std::cos(corner.faceAngle);
corner.sinFaceAngle = std::sin(corner.faceAngle);
}
corner.ringPoints.SetSize(sourcePatch.GetCornerRingSize(cIndex));
@ -1047,11 +1062,11 @@ GregoryConverter<REAL>::computeIrregularEdgePoints(int cIndex,
p.Append(cIndex, 1.0f);
// Approximating these for now, pending future investigation...
ep.Append(cIndex, 2.0f / 3.0f);
ep.Append((cIndex+1) & 0x3, 1.0f / 3.0f);
ep.Append(cIndex, (REAL)(2.0 / 3.0));
ep.Append((cIndex+1) & 0x3, (REAL)(1.0 / 3.0));
em.Append(cIndex, 2.0f / 3.0f);
em.Append((cIndex+3) & 0x3, 1.0f / 3.0f);
em.Append(cIndex, (REAL)(2.0 / 3.0));
em.Append((cIndex+3) & 0x3, (REAL)(1.0 / 3.0));
} else if (! corner.isBoundary) {
//
// The irregular interior case:
@ -1066,15 +1081,15 @@ GregoryConverter<REAL>::computeIrregularEdgePoints(int cIndex,
//
// The irregular/smooth corner case:
//
p.Append(cIndex, 4.0f / 6.0f);
p.Append((cIndex+1) & 0x3, 1.0f / 6.0f);
p.Append((cIndex+3) & 0x3, 1.0f / 6.0f);
p.Append(cIndex, (REAL)(4.0 / 6.0));
p.Append((cIndex+1) & 0x3, (REAL)(1.0 / 6.0));
p.Append((cIndex+3) & 0x3, (REAL)(1.0 / 6.0));
ep.Append(cIndex, 2.0f / 3.0f);
ep.Append((cIndex+1) & 0x3, 1.0f / 3.0f);
ep.Append(cIndex, (REAL)(2.0 / 3.0));
ep.Append((cIndex+1) & 0x3, (REAL)(1.0 / 3.0));
em.Append(cIndex, 2.0f / 3.0f);
em.Append((cIndex+3) & 0x3, 1.0f / 3.0f);
em.Append(cIndex, (REAL)(2.0 / 3.0));
em.Append((cIndex+3) & 0x3, (REAL)(1.0 / 3.0));
}
assert(matrix.GetRowSize(5*cIndex + 0) == p.GetSize());
@ -1288,17 +1303,17 @@ GregoryConverter<REAL>::assignRegularFacePoints(int cIndex, Matrix & matrix) con
// Assign regular Fp and/or Fm:
if (corner.fpIsRegular) {
fp.Append(cIndex, 4.0f / 9.0f);
fp.Append(cPrev, 2.0f / 9.0f);
fp.Append(cNext, 2.0f / 9.0f);
fp.Append(cOpp, 1.0f / 9.0f);
fp.Append(cIndex, (REAL)(4.0 / 9.0));
fp.Append(cPrev, (REAL)(2.0 / 9.0));
fp.Append(cNext, (REAL)(2.0 / 9.0));
fp.Append(cOpp, (REAL)(1.0 / 9.0));
assert(matrix.GetRowSize(5*cIndex + 3) == fp.GetSize());
}
if (corner.fmIsRegular) {
fm.Append(cIndex, 4.0f / 9.0f);
fm.Append(cPrev, 2.0f / 9.0f);
fm.Append(cNext, 2.0f / 9.0f);
fm.Append(cOpp, 1.0f / 9.0f);
fm.Append(cIndex, (REAL)(4.0 / 9.0));
fm.Append(cPrev, (REAL)(2.0 / 9.0));
fm.Append(cNext, (REAL)(2.0 / 9.0));
fm.Append(cOpp, (REAL)(1.0 / 9.0));
assert(matrix.GetRowSize(5*cIndex + 4) == fm.GetSize());
}
}
@ -1692,7 +1707,7 @@ BSplineConverter<REAL>::convertIrregularCorner(int irregularCorner,
wX0[p3inRing] = 4.0f;
// Combine weights for all X[] in one iteration through the ring:
const REAL oneThird = 1.0f / 3.0f;
const REAL oneThird = (REAL) (1.0 / 3.0);
for (int i = 0; i < ringSizePlusCorner; ++i) {
wX1[i] = (36.0f * wEp[i] - wX1[i]) * oneThird;
wX2[i] = (36.0f * wEm[i] - wX2[i]) * oneThird;
@ -1898,10 +1913,11 @@ namespace {
};
}
template <typename REAL>
int
CatmarkPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
CatmarkPatchBuilder::convertSourcePatch(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const {
SparseMatrix<REAL> & matrix) const {
assert(_schemeType == Sdc::SCHEME_CATMARK);
@ -1912,17 +1928,30 @@ CatmarkPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
//
if (patchType == PatchDescriptor::GREGORY_BASIS) {
GregoryConverter<float>(sourcePatch, matrix);
GregoryConverter<REAL>(sourcePatch, matrix);
} else if (patchType == PatchDescriptor::REGULAR) {
BSplineConverter<float>(sourcePatch, matrix);
BSplineConverter<REAL>(sourcePatch, matrix);
} else if (patchType == PatchDescriptor::QUADS) {
LinearConverter<float>(sourcePatch, matrix);
LinearConverter<REAL>(sourcePatch, matrix);
} else {
assert("Unknown or unsupported patch type" == 0);
}
return matrix.GetNumRows();
}
int
CatmarkPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const {
return convertSourcePatch(sourcePatch, patchType, matrix);
}
int
CatmarkPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<double> & matrix) const {
return convertSourcePatch(sourcePatch, patchType, matrix);
}
CatmarkPatchBuilder::CatmarkPatchBuilder(
TopologyRefiner const& refiner, Options const& options) :
PatchBuilder(refiner, options) {

8
opensubdiv/far/catmarkPatchBuilder.h
View File

@ -53,9 +53,17 @@ protected:
virtual int convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const;
virtual int convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<double> & matrix) const;
private:
typedef SparseMatrix<float> ConversionMatrix;
template <typename REAL>
int convertSourcePatch(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<REAL> & matrix) const;
};
} // end namespace Far

21
opensubdiv/far/loopPatchBuilder.cpp
View File

@ -72,10 +72,11 @@ LoopPatchBuilder::patchTypeFromBasis(BasisType basis) const {
return patchTypeFromBasisArray[(int)basis];
}
template <typename REAL>
int
LoopPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const {
LoopPatchBuilder::convertSourcePatch(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<REAL> & matrix) const {
assert("Non-linear patches for Loop not yet supported" == 0);
@ -86,6 +87,20 @@ LoopPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
return -1;
}
int
LoopPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const {
return convertSourcePatch(sourcePatch, patchType, matrix);
}
int
LoopPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<double> & matrix) const {
return convertSourcePatch(sourcePatch, patchType, matrix);
}
} // end namespace Far
} // end namespace OPENSUBDIV_VERSION

8
opensubdiv/far/loopPatchBuilder.h
View File

@ -53,6 +53,14 @@ protected:
virtual int convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const;
virtual int convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<double> & matrix) const;
private:
template <typename REAL>
int convertSourcePatch(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<REAL> & matrix) const;
};
} // end namespace Far

460
opensubdiv/far/patchBasis.cpp
View File

@ -33,104 +33,79 @@ namespace OPENSUBDIV_VERSION {
namespace Far {
namespace internal {
enum SplineBasis {
BASIS_BILINEAR,
BASIS_BEZIER,
BASIS_BSPLINE,
BASIS_BOX_SPLINE
};
template <SplineBasis BASIS>
class Spline {
public:
// curve weights
static void GetWeights(float t, float point[], float deriv[], float deriv2[]);
// box-spline weights
static void GetWeights(float v, float w, float point[]);
// patch weights
static void GetPatchWeights(PatchParam const & param,
float s, float t, float point[], float deriv1[], float deriv2[], float deriv11[], float deriv12[], float deriv22[]);
// adjust patch weights for boundary (and corner) edges
static void AdjustBoundaryWeights(PatchParam const & param,
float sWeights[4], float tWeights[4]);
};
template <>
inline void Spline<BASIS_BEZIER>::GetWeights(
float t, float point[4], float deriv[4], float deriv2[4]) {
namespace {
//
// Evaluation functions for curves used to assemble tensor product bases:
//
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:
float t2 = t*t;
float tC = 1.0f - t;
float tC2 = tC * tC;
REAL t2 = t*t;
REAL tC = 1.0f - t;
REAL tC2 = tC * tC;
assert(point);
point[0] = tC2 * tC;
point[1] = tC2 * t * 3.0f;
point[2] = t2 * tC * 3.0f;
point[3] = t2 * t;
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 (deriv) {
deriv[0] = -3.0f * tC2;
deriv[1] = 9.0f * t2 - 12.0f * t + 3.0f;
deriv[2] = -9.0f * t2 + 6.0f * t;
deriv[3] = 3.0f * t2;
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 (deriv2) {
deriv2[0] = 6.0f * tC;
deriv2[1] = 18.0f * t - 12.0f;
deriv2[2] = -18.0f * t + 6.0f;
deriv2[3] = 6.0f * 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;
}
}
template <>
inline void Spline<BASIS_BSPLINE>::GetWeights(
float t, float point[4], float deriv[4], float deriv2[4]) {
template <typename REAL>
void evalBSplineCurve(REAL t, REAL wP[4], REAL wDP[4], REAL wDP2[4]) {
// The four uniform cubic B-Spline basis functions evaluated at t:
float const one6th = 1.0f / 6.0f;
REAL const one6th = (REAL)(1.0 / 6.0);
float t2 = t * t;
float t3 = t * t2;
REAL t2 = t * t;
REAL t3 = t * t2;
assert(point);
point[0] = one6th * (1.0f - 3.0f*(t - t2) - t3);
point[1] = one6th * (4.0f - 6.0f*t2 + 3.0f*t3);
point[2] = one6th * (1.0f + 3.0f*(t + t2 - t3));
point[3] = one6th * ( t3);
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);
// Derivatives of the above four basis functions at t:
if (deriv) {
deriv[0] = -0.5f*t2 + t - 0.5f;
deriv[1] = 1.5f*t2 - 2.0f*t;
deriv[2] = -1.5f*t2 + t + 0.5f;
deriv[3] = 0.5f*t2;
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;
}
// Second derivatives of the basis functions at t:
if (deriv2) {
deriv2[0] = - t + 1.0f;
deriv2[1] = 3.0f * t - 2.0f;
deriv2[2] = -3.0f * t + 1.0f;
deriv2[3] = t;
if (wDP2) {
wDP2[0] = - t + 1.0f;
wDP2[1] = 3.0f * t - 2.0f;
wDP2[2] = -3.0f * t + 1.0f;
wDP2[3] = t;
}
}
}
template <>
inline void Spline<BASIS_BOX_SPLINE>::GetWeights(
float v, float w, float point[12]) {
template <typename REAL>
void GetBoxSplineWeights(PatchParam const & param, REAL s, REAL t, REAL wP[12]) {
float u = 1.0f - v - w;
float u = s;
float v = t;
float w = 1.0f - u - v;
//
// The 12 basis functions of the quartic box spline (unscaled by their common
@ -153,180 +128,208 @@ inline void Spline<BASIS_BOX_SPLINE>::GetWeights(
float w4 = w*w3;
// And now the basis functions:
point[ 0] = u4 + 2.0f*u3*v;
point[ 1] = u4 + 2.0f*u3*w;
point[ 8] = w4 + 2.0f*w3*u;
point[11] = w4 + 2.0f*w3*v;
point[ 9] = v4 + 2.0f*v3*w;
point[ 5] = v4 + 2.0f*v3*u;
wP[ 0] = u4 + 2.0f*u3*v;
wP[ 1] = u4 + 2.0f*u3*w;
wP[ 8] = w4 + 2.0f*w3*u;
wP[11] = w4 + 2.0f*w3*v;
wP[ 9] = v4 + 2.0f*v3*w;
wP[ 5] = v4 + 2.0f*v3*u;
point[ 2] = u4 + 2.0f*u3*w + 6.0f*u3*v + 6.0f*u2*v*w + 12.0f*u2*v2 +
wP[ 2] = u4 + 2.0f*u3*w + 6.0f*u3*v + 6.0f*u2*v*w + 12.0f*u2*v2 +
v4 + 2.0f*v3*w + 6.0f*v3*u + 6.0f*v2*u*w;
point[ 4] = w4 + 2.0f*w3*v + 6.0f*w3*u + 6.0f*w2*u*v + 12.0f*w2*u2 +
wP[ 4] = w4 + 2.0f*w3*v + 6.0f*w3*u + 6.0f*w2*u*v + 12.0f*w2*u2 +
u4 + 2.0f*u3*v + 6.0f*u3*w + 6.0f*u2*v*w;
point[10] = v4 + 2.0f*v3*u + 6.0f*v3*w + 6.0f*v2*w*u + 12.0f*v2*w2 +
wP[10] = v4 + 2.0f*v3*u + 6.0f*v3*w + 6.0f*v2*w*u + 12.0f*v2*w2 +
w4 + 2.0f*w3*u + 6.0f*w3*v + 6.0f*w3*u*v;
point[ 3] = v4 + 6*v3*w + 8*v3*u + 36*v2*w*u + 24*v2*u2 + 24*v*u3 +
wP[ 3] = v4 + 6*v3*w + 8*v3*u + 36*v2*w*u + 24*v2*u2 + 24*v*u3 +
w4 + 6*w3*v + 8*w3*u + 36*w2*v*u + 24*w2*u2 + 24*w*u3 + 6*u4 + 60*u2*v*w + 12*v2*w2;
point[ 6] = w4 + 6*w3*u + 8*w3*v + 36*w2*u*v + 24*w2*v2 + 24*w*v3 +
wP[ 6] = w4 + 6*w3*u + 8*w3*v + 36*w2*u*v + 24*w2*v2 + 24*w*v3 +
u4 + 6*u3*w + 8*u3*v + 36*u2*v*w + 24*u2*v2 + 24*u*v3 + 6*v4 + 60*v2*w*u + 12*w2*u2;
point[ 7] = u4 + 6*u3*v + 8*u3*w + 36*u2*v*w + 24*u2*w2 + 24*u*w3 +
wP[ 7] = u4 + 6*u3*v + 8*u3*w + 36*u2*v*w + 24*u2*w2 + 24*u*w3 +
v4 + 6*v3*u + 8*v3*w + 36*v2*u*w + 24*v2*w2 + 24*v*w3 + 6*w4 + 60*w2*u*v + 12*u2*v2;
for (int i = 0; i < 12; ++i) {
point[i] *= 1.0f / 12.0f;
wP[i] *= 1.0f / 12.0f;
}
}
template <>
inline void Spline<BASIS_BILINEAR>::GetPatchWeights(PatchParam const & param,
float s, float t, float point[4], float derivS[4], float derivT[4], float derivSS[4], float derivST[4], float derivTT[4]) {
template <typename REAL>
void GetBilinearWeights(PatchParam const & param, REAL s, REAL t,
REAL wP[4], REAL wDs[4], REAL wDt[4],
REAL wDss[4], REAL wDst[4], REAL wDtt[4]) {
param.Normalize(s,t);
float sC = 1.0f - s,
tC = 1.0f - t;
REAL sC = 1.0f - s;
REAL tC = 1.0f - t;
if (point) {
point[0] = sC * tC;
point[1] = s * tC;
point[2] = s * t;
point[3] = sC * t;
if (wP) {
wP[0] = sC * tC;
wP[1] = s * tC;
wP[2] = s * t;
wP[3] = sC * t;
}
if (derivS && derivT) {
float dScale = (float)(1 << param.GetDepth());
derivS[0] = -tC * dScale;
derivS[1] = tC * dScale;
derivS[2] = t * dScale;
derivS[3] = -t * dScale;
if (wDs && wDt) {
REAL dScale = (REAL)(1 << param.GetDepth());
derivT[0] = -sC * dScale;
derivT[1] = -s * dScale;
derivT[2] = s * dScale;
derivT[3] = sC * dScale;
wDs[0] = -tC * dScale;
wDs[1] = tC * dScale;
wDs[2] = t * dScale;
wDs[3] = -t * dScale;
if (derivSS && derivST && derivTT) {
float d2Scale = dScale * dScale;
wDt[0] = -sC * dScale;
wDt[1] = -s * dScale;
wDt[2] = s * dScale;
wDt[3] = sC * dScale;
if (wDss && wDst && wDtt) {
REAL d2Scale = dScale * dScale;
for(int i=0;i<4;i++) {
derivSS[i] = 0;
derivTT[i] = 0;
wDss[i] = 0.0f;
wDtt[i] = 0.0f;
}
derivST[0] = d2Scale;
derivST[1] = -d2Scale;
derivST[2] = -d2Scale;
derivST[3] = d2Scale;
wDst[0] = d2Scale;
wDst[1] = -d2Scale;
wDst[2] = -d2Scale;
wDst[3] = d2Scale;
}
}
}
template <SplineBasis BASIS>
void Spline<BASIS>::AdjustBoundaryWeights(PatchParam const & param,
float sWeights[4], float tWeights[4]) {
//
// BSpline patch evaluation -- involves adjustments to weights when boundary
// points are missing and implicitly extrapolated.
//
namespace {
template <typename REAL>
void adjustBSplineBoundaryWeights(PatchParam const & param, REAL sWeights[4], REAL tWeights[4]) {
int boundary = param.GetBoundary();
if (boundary & 1) {
if ((boundary & 1) != 0) {
tWeights[2] -= tWeights[0];
tWeights[1] += 2*tWeights[0];
tWeights[0] = 0;
tWeights[1] += tWeights[0] * 2.0f;
tWeights[0] = 0.0f;
}
if (boundary & 2) {
if ((boundary & 2) != 0) {
sWeights[1] -= sWeights[3];
sWeights[2] += 2*sWeights[3];
sWeights[3] = 0;
sWeights[2] += sWeights[3] * 2.0f;
sWeights[3] = 0.0f;
}
if (boundary & 4) {
if ((boundary & 4) != 0) {
tWeights[1] -= tWeights[3];
tWeights[2] += 2*tWeights[3];
tWeights[3] = 0;
tWeights[2] += tWeights[3] * 2.0f;
tWeights[3] = 0.0f;
}
if (boundary & 8) {
if ((boundary & 8) != 0) {
sWeights[2] -= sWeights[0];
sWeights[1] += 2*sWeights[0];
sWeights[0] = 0;
sWeights[1] += sWeights[0] * 2.0f;
sWeights[0] = 0.0f;
}
}
}
template <SplineBasis BASIS>
void Spline<BASIS>::GetPatchWeights(PatchParam const & param,
float s, float t, float point[16], float derivS[16], float derivT[16], float derivSS[16], float derivST[16], float derivTT[16]) {
template <typename REAL>
void GetBSplineWeights(PatchParam const & param, REAL s, REAL t,
REAL wP[16], REAL wDs[16], REAL wDt[16],
REAL wDss[16], REAL wDst[16], REAL wDtt[16]) {
float sWeights[4], tWeights[4], dsWeights[4], dtWeights[4], dssWeights[4], dttWeights[4];
REAL sWeights[4], tWeights[4], dsWeights[4], dtWeights[4], dssWeights[4], dttWeights[4];
param.Normalize(s,t);
Spline<BASIS>::GetWeights(s, point ? sWeights : 0, derivS ? dsWeights : 0, derivSS ? dssWeights : 0);
Spline<BASIS>::GetWeights(t, point ? tWeights : 0, derivT ? dtWeights : 0, derivTT ? dttWeights : 0);
evalBSplineCurve(s, wP ? sWeights : 0, wDs ? dsWeights : 0, wDss ? dssWeights : 0);
evalBSplineCurve(t, wP ? tWeights : 0, wDt ? dtWeights : 0, wDtt ? dttWeights : 0);
if (point) {
// Compute the tensor product weight of the (s,t) basis function
// corresponding to each control vertex:
AdjustBoundaryWeights(param, sWeights, tWeights);
if (wP) {
adjustBSplineBoundaryWeights(param, sWeights, tWeights);
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
point[4*i+j] = sWeights[j] * tWeights[i];
wP[4*i+j] = sWeights[j] * tWeights[i];
}
}
}
if (derivS && derivT) {
// Compute the tensor product weight of the differentiated (s,t) basis
// function corresponding to each control vertex (scaled accordingly):
if (wDs && wDt) {
REAL dScale = (REAL)(1 << param.GetDepth());
float dScale = (float)(1 << param.GetDepth());
AdjustBoundaryWeights(param, dsWeights, dtWeights);
adjustBSplineBoundaryWeights(param, dsWeights, dtWeights);
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
derivS[4*i+j] = dsWeights[j] * tWeights[i] * dScale;
derivT[4*i+j] = sWeights[j] * dtWeights[i] * dScale;
wDs[4*i+j] = dsWeights[j] * tWeights[i] * dScale;
wDt[4*i+j] = sWeights[j] * dtWeights[i] * dScale;
}
}
if (derivSS && derivST && derivTT) {
// Compute the tensor product weight of appropriate differentiated
// (s,t) basis functions for each control vertex (scaled accordingly):
float d2Scale = dScale * dScale;
if (wDss && wDst && wDtt) {
REAL d2Scale = dScale * dScale;
AdjustBoundaryWeights(param, dssWeights, dttWeights);
adjustBSplineBoundaryWeights(param, dssWeights, dttWeights);
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
derivSS[4*i+j] = dssWeights[j] * tWeights[i] * d2Scale;
derivST[4*i+j] = dsWeights[j] * dtWeights[i] * d2Scale;
derivTT[4*i+j] = sWeights[j] * dttWeights[i] * d2Scale;
wDss[4*i+j] = dssWeights[j] * tWeights[i] * d2Scale;
wDst[4*i+j] = dsWeights[j] * dtWeights[i] * d2Scale;
wDtt[4*i+j] = sWeights[j] * dttWeights[i] * d2Scale;
}
}
}
}
}
void GetBilinearWeights(PatchParam const & param,
float s, float t, float point[4], float deriv1[4], float deriv2[4], float deriv11[4], float deriv12[4], float deriv22[4]) {
Spline<BASIS_BILINEAR>::GetPatchWeights(param, s, t, point, deriv1, deriv2, deriv11, deriv12, deriv22);
template <typename REAL>
void GetBezierWeights(PatchParam const & param, 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];
param.Normalize(s,t);
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) {
REAL dScale = (REAL)(1 << param.GetDepth());
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
wDs[4*i+j] = dsWeights[j] * tWeights[i] * dScale;
wDt[4*i+j] = sWeights[j] * dtWeights[i] * dScale;
}
}
if (wDss && wDst && wDtt) {
REAL d2Scale = dScale * dScale;
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
wDss[4*i+j] = dssWeights[j] * tWeights[i] * d2Scale;
wDst[4*i+j] = dsWeights[j] * dtWeights[i] * d2Scale;
wDtt[4*i+j] = sWeights[j] * dttWeights[i] * d2Scale;
}
}
}
}
}
void GetBezierWeights(PatchParam const param,
float s, float t, float point[16], float deriv1[16], float deriv2[16], float deriv11[16], float deriv12[16], float deriv22[16]) {
Spline<BASIS_BEZIER>::GetPatchWeights(param, s, t, point, deriv1, deriv2, deriv11, deriv12, deriv22);
}
template <typename REAL>
void GetGregoryWeights(PatchParam const & param, REAL s, REAL t,
REAL point[20], REAL wDs[20], REAL wDt[20],
REAL wDss[20], REAL wDst[20], REAL wDtt[20]) {
void GetBSplineWeights(PatchParam const & param,
float s, float t, float point[16], float deriv1[16], float deriv2[16], float deriv11[16], float deriv12[16], float deriv22[16]) {
Spline<BASIS_BSPLINE>::GetPatchWeights(param, s, t, point, deriv1, deriv2, deriv11, deriv12, deriv22);
}
void GetGregoryWeights(PatchParam const & param,
float s, float t, float point[20], float deriv1[20], float deriv2[20], float deriv11[20], float deriv12[20], float deriv22[20]) {
//
// P3 e3- e2+ P2
// 15------17-------11--------10
@ -364,25 +367,25 @@ void GetGregoryWeights(PatchParam const & param,
// 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:
float Bs[4], Bds[4], Bdss[4];
float Bt[4], Bdt[4], Bdtt[4];
REAL Bs[4], Bds[4], Bdss[4];
REAL Bt[4], Bdt[4], Bdtt[4];
param.Normalize(s,t);
Spline<BASIS_BEZIER>::GetWeights(s, Bs, deriv1 ? Bds : 0, deriv11 ? Bdss : 0);
Spline<BASIS_BEZIER>::GetWeights(t, Bt, deriv2 ? Bdt : 0, deriv22 ? Bdtt : 0);
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:
float sC = 1.0f - s;
float tC = 1.0f - t;
REAL sC = 1.0f - s;
REAL tC = 1.0f - t;
// Use <= here to avoid compiler warnings -- the sums should always be non-negative:
float df0 = s + t; df0 = (df0 <= 0.0f) ? 1.0f : (1.0f / df0);
float df1 = sC + t; df1 = (df1 <= 0.0f) ? 1.0f : (1.0f / df1);
float df2 = sC + tC; df2 = (df2 <= 0.0f) ? 1.0f : (1.0f / df2);
float df3 = s + tC; df3 = (df3 <= 0.0f) ? 1.0f : (1.0f / df3);
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);
float G[8] = { s*df0, t*df0, t*df1, sC*df1, sC*df2, tC*df2, tC*df3, s*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) {
@ -403,11 +406,11 @@ void GetGregoryWeights(PatchParam const & param,
// unclear if the approximations will hold up under surface analysis involving higher
// order differentiation.
//
if (deriv1 && deriv2) {
bool find_second_partials = deriv1 && deriv12 && deriv22;
if (wDs && wDt) {
bool find_second_partials = wDs && wDst && wDtt;
// Remember to include derivative scaling in all assignments below:
float dScale = (float)(1 << param.GetDepth());
float d2Scale = dScale * dScale;
REAL dScale = (REAL)(1 << param.GetDepth());
REAL d2Scale = dScale * dScale;
// Combined weights for boundary points -- simple (scaled) tensor products:
for (int i = 0; i < 12; ++i) {
@ -415,13 +418,13 @@ void GetGregoryWeights(PatchParam const & param,
int tRow = boundaryBezTRow[i];
int sCol = boundaryBezSCol[i];
deriv1[iDst] = Bds[sCol] * Bt[tRow] * dScale;
deriv2[iDst] = Bdt[tRow] * Bs[sCol] * dScale;
wDs[iDst] = Bds[sCol] * Bt[tRow] * dScale;
wDt[iDst] = Bdt[tRow] * Bs[sCol] * dScale;
if (find_second_partials) {
deriv11[iDst] = Bdss[sCol] * Bt[tRow] * d2Scale;
deriv12[iDst] = Bds[sCol] * Bdt[tRow] * d2Scale;
deriv22[iDst] = Bs[sCol] * Bdtt[tRow] * d2Scale;
wDss[iDst] = Bdss[sCol] * Bt[tRow] * d2Scale;
wDst[iDst] = Bds[sCol] * Bdt[tRow] * d2Scale;
wDtt[iDst] = Bs[sCol] * Bdtt[tRow] * d2Scale;
}
}
@ -438,13 +441,13 @@ void GetGregoryWeights(PatchParam const & param,
int tRow = interiorBezTRow[i];
int sCol = interiorBezSCol[i];
deriv1[iDst] = Bds[sCol] * Bt[tRow] * G[i] * dScale;
deriv2[iDst] = Bdt[tRow] * Bs[sCol] * G[i] * dScale;
wDs[iDst] = Bds[sCol] * Bt[tRow] * G[i] * dScale;
wDt[iDst] = Bdt[tRow] * Bs[sCol] * G[i] * dScale;
if (find_second_partials) {
deriv11[iDst] = Bdss[sCol] * Bt[tRow] * G[i] * d2Scale;
deriv12[iDst] = Bds[sCol] * Bdt[tRow] * G[i] * d2Scale;
deriv22[iDst] = Bs[sCol] * Bdtt[tRow] * G[i] * d2Scale;
wDss[iDst] = Bdss[sCol] * Bt[tRow] * G[i] * d2Scale;
wDst[iDst] = Bds[sCol] * Bdt[tRow] * G[i] * d2Scale;
wDtt[iDst] = Bs[sCol] * Bdtt[tRow] * G[i] * d2Scale;
}
}
#else
@ -458,14 +461,14 @@ void GetGregoryWeights(PatchParam const & param,
// (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.
//
//float N[8] = { s, t, t, sC, sC, tC, tC, s };
float D[8] = { df0, df0, df1, df1, df2, df2, df3, df3 };
//REAL N[8] = { s, t, t, sC, sC, tC, tC, s };
REAL D[8] = { df0, df0, df1, df1, df2, df2, df3, df3 };
static float const Nds[8] = { 1.0f, 0.0f, 0.0f, -1.0f, -1.0f, 0.0f, 0.0f, 1.0f };
static float const Ndt[8] = { 0.0f, 1.0f, 1.0f, 0.0f, 0.0f, -1.0f, -1.0f, 0.0f };
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 float const Dds[8] = { 1.0f, 1.0f, -1.0f, -1.0f, -1.0f, -1.0f, 1.0f, 1.0f };
static float const Ddt[8] = { 1.0f, 1.0f, 1.0f, 1.0f, -1.0f, -1.0f, -1.0f, -1.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 -- (scaled) combinations of B, B', G and G':
for (int i = 0; i < 8; ++i) {
@ -474,29 +477,54 @@ void GetGregoryWeights(PatchParam const & param,
int sCol = interiorBezSCol[i];
// Quotient rule for G' (re-expressed in terms of G to simplify (and D = 1/D)):
float Gds = (Nds[i] - Dds[i] * G[i]) * D[i];
float Gdt = (Ndt[i] - Ddt[i] * G[i]) * D[i];
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' (and scaled):
deriv1[iDst] = (Bds[sCol] * G[i] + Bs[sCol] * Gds) * Bt[tRow] * dScale;
deriv2[iDst] = (Bdt[tRow] * G[i] + Bt[tRow] * Gdt) * Bs[sCol] * dScale;
wDs[iDst] = (Bds[sCol] * G[i] + Bs[sCol] * Gds) * Bt[tRow] * dScale;
wDt[iDst] = (Bdt[tRow] * G[i] + Bt[tRow] * Gdt) * Bs[sCol] * dScale;
if (find_second_partials) {
float Dsqr_inv = D[i]*D[i];
REAL Dsqr_inv = D[i]*D[i];
float Gdss = 2.0f * Dds[i] * Dsqr_inv * (G[i] * Dds[i] - Nds[i]);
float Gdst = Dsqr_inv * (2.0f * G[i] * Dds[i] * Ddt[i] - Nds[i] * Ddt[i] - Ndt[i] * Dds[i]);
float Gdtt = 2.0f * Ddt[i] * Dsqr_inv * (G[i] * Ddt[i] - Ndt[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]);
deriv11[iDst] = (Bdss[sCol] * G[i] + 2.0f * Bds[sCol] * Gds + Bs[sCol] * Gdss) * Bt[tRow] * d2Scale;
deriv12[iDst] = (Bt[tRow] * (Bs[sCol] * Gdst + Bds[sCol] * Gdt) + Bdt[tRow] * (Bds[sCol] * G[i] + Bs[sCol] * Gds)) * d2Scale;
deriv22[iDst] = (Bdtt[tRow] * G[i] + 2.0f * Bdt[tRow] * Gdt + Bt[tRow] * Gdtt) * Bs[sCol] * d2Scale;
wDss[iDst] = (Bdss[sCol] * G[i] + 2.0f * Bds[sCol] * Gds + Bs[sCol] * Gdss) * Bt[tRow] * d2Scale;
wDst[iDst] = (Bt[tRow] * (Bs[sCol] * Gdst + Bds[sCol] * Gdt) + Bdt[tRow] * (Bds[sCol] * G[i] + Bs[sCol] * Gds)) * d2Scale;
wDtt[iDst] = (Bdtt[tRow] * G[i] + 2.0f * Bdt[tRow] * Gdt + Bt[tRow] * Gdtt) * Bs[sCol] * d2Scale;
}
}
#endif
}
}
//
// Explicit float and double instantiations:
//
template void GetBilinearWeights<float>(PatchParam const & patchParam, float s, float t,
float wP[4], float wDs[4], float wDt[4], float wDss[4], float wDst[4], float wDtt[4]);
template void GetBezierWeights<float>(PatchParam const & patchParam, float s, float t,
float wP[16], float wDs[16], float wDt[16], float wDss[16], float wDst[16], float wDtt[16]);
template void GetBSplineWeights<float>(PatchParam const & patchParam, float s, float t,
float wP[16], float wDs[16], float wDt[16], float wDss[16], float wDst[16], float wDtt[16]);
template void GetGregoryWeights<float>(PatchParam const & patchParam, float s, float t,
float wP[20], float wDs[20], float wDt[20], float wDss[20], float wDst[20], float wDtt[20]);
// Cannot enable these until PatchParam::Normalize() et al are extended...
/*
template void GetBilinearWeights<double>(PatchParam const & patchParam, double s, double t,
double wP[4], double wDs[4], double wDt[4], double wDss[4], double wDst[4], double wDtt[4]);
template void GetBezierWeights<double>(PatchParam const & patchParam, double s, double t,
double wP[16], double wDs[16], double wDt[16], double wDss[16], double wDst[16], double wDtt[16]);
template void GetBSplineWeights<double>(PatchParam const & patchParam, double s, double t,
double wP[16], double wDs[16], double wDt[16], double wDss[16], double wDst[16], double wDtt[16]);
template void GetGregoryWeights<double>(PatchParam const & patchParam, double s, double t,
double wP[20], double wDs[20], double wDt[20], double wDss[20], double wDst[20], double wDtt[20]);
*/
} // end namespace internal
} // end namespace Far

20
opensubdiv/far/patchBasis.h
View File

@ -45,17 +45,21 @@ namespace internal {
// So this interface will be changing in future.
//
void GetBilinearWeights(PatchParam const & patchParam,
float s, float t, float wP[4], float wDs[4], float wDt[4], float wDss[4] = 0, float wDst[4] = 0, float wDtt[4] = 0);
template <typename REAL>
void GetBilinearWeights(PatchParam const & patchParam, REAL s, REAL t,
REAL wP[4], REAL wDs[4], REAL wDt[4], REAL wDss[4] = 0, REAL wDst[4] = 0, REAL wDtt[4] = 0);
void GetBezierWeights(PatchParam const & patchParam,
float s, float t, float wP[16], float wDs[16], float wDt[16], float wDss[16] = 0, float wDst[16] = 0, float wDtt[16] = 0);
template <typename REAL>
void GetBezierWeights(PatchParam const & patchParam, REAL s, REAL t,
REAL wP[16], REAL wDs[16], REAL wDt[16], REAL wDss[16] = 0, REAL wDst[16] = 0, REAL wDtt[16] = 0);
void GetBSplineWeights(PatchParam const & patchParam,
float s, float t, float wP[16], float wDs[16], float wDt[16], float wDss[16] = 0, float wDst[16] = 0, float wDtt[16] = 0);
template <typename REAL>
void GetBSplineWeights(PatchParam const & patchParam, REAL s, REAL t,
REAL wP[16], REAL wDs[16], REAL wDt[16], REAL wDss[16] = 0, REAL wDst[16] = 0, REAL wDtt[16] = 0);
void GetGregoryWeights(PatchParam const & patchParam,
float s, float t, float wP[20], float wDs[20], float wDt[20], float wDss[20] = 0, float wDst[20] = 0, float wDtt[20] = 0);
template <typename REAL>
void GetGregoryWeights(PatchParam const & patchParam, REAL s, REAL t,
REAL wP[20], REAL wDs[20], REAL wDt[20], REAL wDss[20] = 0, REAL wDst[20] = 0, REAL wDtt[20] = 0);
} // end namespace internal

13
opensubdiv/far/patchBuilder.cpp
View File

@ -798,11 +798,16 @@ PatchBuilder::GetIrregularPatchSourcePoints(
cornerSpans, sourcePatch, sourcePoints, fvarChannel);
}
//
// Template conversion methods for the matrix type -- explicit instantiation
// for float and double is required and follows the definition:
//
template <typename REAL>
int
PatchBuilder::GetIrregularPatchConversionMatrix(
int levelIndex, Index faceIndex,
Level::VSpan const cornerSpans[],
SparseMatrix<float> & conversionMatrix) const {
SparseMatrix<REAL> & conversionMatrix) const {
assertTriangularPatchesNotYetSupportedHere();
@ -813,6 +818,12 @@ PatchBuilder::GetIrregularPatchConversionMatrix(
return convertToPatchType(
sourcePatch, GetIrregularPatchType(), conversionMatrix);
}
template int PatchBuilder::GetIrregularPatchConversionMatrix<float>(
int levelIndex, Index faceIndex, Level::VSpan const cornerSpans[],
SparseMatrix<float> & conversionMatrix) const;
template int PatchBuilder::GetIrregularPatchConversionMatrix<double>(
int levelIndex, Index faceIndex, Level::VSpan const cornerSpans[],
SparseMatrix<double> & conversionMatrix) const;
bool

7
opensubdiv/far/patchBuilder.h
View File

@ -226,9 +226,10 @@ public:
Index patchPoints[],
int fvc = -1) const;
template <typename REAL>
int GetIrregularPatchConversionMatrix(int level, Index face,
Vtr::internal::Level::VSpan const cornerSpans[],
SparseMatrix<float> & matrix) const;
SparseMatrix<REAL> & matrix) const;
int GetIrregularPatchSourcePoints(int level, Index face,
Vtr::internal::Level::VSpan const cornerSpans[],
@ -293,9 +294,13 @@ protected:
//
virtual PatchDescriptor::Type patchTypeFromBasis(BasisType basis) const = 0;
// Note overloading of the conversion for SparseMatrix<REAL>:
virtual int convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<float> & matrix) const = 0;
virtual int convertToPatchType(SourcePatch const & sourcePatch,
PatchDescriptor::Type patchType,
SparseMatrix<double> & matrix) const = 0;
protected:
TopologyRefiner const& _refiner;

141
opensubdiv/far/patchTableFactory.cpp
View File

@ -123,7 +123,6 @@ private:
// Struct comprising a collection of topological properties for a patch that
// may be shared (between vertex and face-varying patches):
typedef SparseMatrix<float> Matrix;
struct PatchInfo {
PatchInfo() : isRegular(false), isRegSingleCrease(false),
regBoundaryMask(0), regSharpness(0.0f),
@ -135,7 +134,8 @@ private:
float regSharpness;
Level::VSpan irregCornerSpans[4];
int paramBoundaryMask;
Matrix matrix;
SparseMatrix<float> matrix;
};
private:
@ -160,11 +160,13 @@ private:
struct Options {
Options() : shareLocalPoints(false),
reuseSourcePoints(false),
createStencilTable(true) { }
createStencilTable(true),
createVaryingTable(false) { }
unsigned int shareLocalPoints : 1;
unsigned int reuseSourcePoints : 1;
unsigned int createStencilTable : 1;
unsigned int createVaryingTable : 1;
};
public:
@ -178,11 +180,11 @@ private:
int GetNumLocalPoints() const { return _numLocalPoints; }
int AppendLocalPatchPoints(int levelIndex, Index faceIndex,
Matrix const & conversionMatrix,
PatchDescriptor::Type patchType,
Index const sourcePoints[],
int sourcePointOffset,
Index patchPoints[]);
SparseMatrix<float> const & conversionMatrix,
PatchDescriptor::Type patchType,
Index const sourcePoints[],
int sourcePointOffset,
Index patchPoints[]);
StencilTable * AcquireStencilTable() {
return acquireStencilTable(&_stencilTable);
@ -190,14 +192,14 @@ private:
private:
// Internal methods:
void appendLocalPointStencil(Matrix const & conversionMatrix,
int stencilRow,
Index const sourcePoints[],
int sourcePointOffset);
void appendLocalPointStencil(SparseMatrix<float> const & conversionMatrix,
int stencilRow,
Index const sourcePoints[],
int sourcePointOffset);
void appendLocalPointStencils(Matrix const & conversionMatrix,
Index const sourcePoints[],
int sourcePointOffset);
void appendLocalPointStencils(SparseMatrix<float> const & conversionMatrix,
Index const sourcePoints[],
int sourcePointOffset);
// Methods for local point Varying stencils
// XXXX -- hope to get rid of these...
@ -206,11 +208,6 @@ private:
Index const sourcePoints[],
int sourcePointOffset);
void appendLocalPointVaryingStencils(int const * varyingIndices,
Matrix const & conversionMatrix,
Index const sourcePoints[],
int sourcePointOffset);
StencilTable * acquireStencilTable(StencilTable **stencilTableMember);
Index findSharedCornerPoint(int levelIndex, Index valueIndex,
@ -1043,16 +1040,18 @@ PatchTableBuilder::populateAdaptivePatches() {
}
//
// Initializing StencilTable and other helpers:
// Initializing StencilTable and other helpers. Note Varying local points
// are managed by the LocalPointHelper for Vertex patches (as Varying local
// points are tightly coupled to their local points)
//
LocalPointHelper * vertexLocalPointHelper = 0;
LocalPointHelper * varyingLocalPointHelper = 0;
StackBuffer<LocalPointHelper*,4> fvarLocalPointHelpers;
if (_requiresLocalPoints) {
LocalPointHelper::Options opts;
opts.createStencilTable = true;
opts.createVaryingTable = _requiresVaryingLocalPoints;
opts.shareLocalPoints = _options.shareEndCapPatchPoints;
opts.reuseSourcePoints = (_patchBuilder->GetIrregularPatchType() ==
_patchBuilder->GetNativePatchType() );
@ -1061,6 +1060,9 @@ PatchTableBuilder::populateAdaptivePatches() {
_refiner, opts, -1, estimateLocalPointCount(opts, -1));
if (_requiresFVarPatches) {
opts.createStencilTable = true;
opts.createVaryingTable = false;
fvarLocalPointHelpers.SetSize((int)_fvarChannelIndices.size());
for (int fvc = 0; fvc < (int)_fvarChannelIndices.size(); ++fvc) {
@ -1172,15 +1174,12 @@ PatchTableBuilder::populateAdaptivePatches() {
// Finalizing and destroying StencilTable and other helpers:
//
if (_requiresLocalPoints) {
if (_requiresVaryingPatches) {
// XXXX (barfowl) - note this is still managed by the vertex helper
_table->_localPointVaryingStencils =
vertexLocalPointHelper->AcquireStencilTableVarying();
delete varyingLocalPointHelper;
}
_table->_localPointStencils =
vertexLocalPointHelper->AcquireStencilTable();
if (_requiresVaryingLocalPoints) {
_table->_localPointVaryingStencils =
vertexLocalPointHelper->AcquireStencilTableVarying();
}
delete vertexLocalPointHelper;
if (_requiresFVarPatches) {
@ -1307,7 +1306,7 @@ PatchTableBuilder::LocalPointHelper::LocalPointHelper(
if (_options.createStencilTable) {
_stencilTable = new StencilTable(0);
if (_stencilTable && (_fvarChannel < 0)) {
if (_options.createVaryingTable) {
_stencilTableVarying = new StencilTable(0);
}
}
@ -1415,10 +1414,10 @@ PatchTableBuilder::LocalPointHelper::findSharedEdgePoint(int levelIndex,
void
PatchTableBuilder::LocalPointHelper::appendLocalPointStencil(
Matrix const & conversionMatrix,
int stencilRow,
Index const sourcePoints[],
int sourcePointOffset) {
SparseMatrix<float> const & conversionMatrix,
int stencilRow,
Index const sourcePoints[],
int sourcePointOffset) {
int stencilSize = conversionMatrix.GetRowSize(stencilRow);
ConstArray<int> matrixColumns = conversionMatrix.GetRowColumns(stencilRow);
@ -1434,9 +1433,9 @@ PatchTableBuilder::LocalPointHelper::appendLocalPointStencil(
void
PatchTableBuilder::LocalPointHelper::appendLocalPointStencils(
Matrix const & conversionMatrix,
Index const sourcePoints[],
int sourcePointOffset) {
SparseMatrix<float> const & conversionMatrix,
Index const sourcePoints[],
int sourcePointOffset) {
//
// Resize the StencilTable members to accomodate all rows and elements from
@ -1519,20 +1518,6 @@ PatchTableBuilder::LocalPointHelper::appendLocalPointVaryingStencil(
_stencilTableVarying->_weights.push_back(1.0f);
}
void
PatchTableBuilder::LocalPointHelper::appendLocalPointVaryingStencils(
int const * varyingIndices, Matrix const & conversionMatrix,
Index const sourcePoints[], int sourcePointOffset) {
for (int i = 0; i < conversionMatrix.GetNumRows(); ++i) {
Index varyingPoint = sourcePoints[varyingIndices[i]] + sourcePointOffset;
_stencilTableVarying->_sizes.push_back(1);
_stencilTableVarying->_indices.push_back(varyingPoint);
_stencilTableVarying->_weights.push_back(1.0f);
}
}
namespace {
typedef short ShareBits;
@ -1565,11 +1550,11 @@ namespace {
int
PatchTableBuilder::LocalPointHelper::AppendLocalPatchPoints(
int levelIndex, Index faceIndex,
Matrix const & matrix,
PatchDescriptor::Type patchType,
Index const sourcePoints[],
int sourcePointOffset,
Index patchPoints[]) {
SparseMatrix<float> const & matrix,
PatchDescriptor::Type patchType,
Index const sourcePoints[],
int sourcePointOffset,
Index patchPoints[]) {
//
// If sharing local points, verify the type of this patch supports
@ -1585,8 +1570,13 @@ PatchTableBuilder::LocalPointHelper::AppendLocalPatchPoints(
bool shareLocalPointsForThisPatch = (shareBitsPerPoint != 0);
int const * varyingIndices =
_stencilTableVarying ? GetVaryingIndicesPerType(patchType) : 0;
int const * varyingIndices = 0;
if (_stencilTableVarying) {
varyingIndices = GetVaryingIndicesPerType(patchType);
}
bool applyVertexStencils = (_stencilTable != 0);
bool applyVaryingStencils = (varyingIndices != 0);
//
// When point-sharing is not enabled, all patch points are generally
@ -1597,10 +1587,15 @@ PatchTableBuilder::LocalPointHelper::AppendLocalPatchPoints(
//
if (!shareLocalPointsForThisPatch) {
if (!_options.reuseSourcePoints) {
appendLocalPointStencils(matrix, sourcePoints, sourcePointOffset);
if (varyingIndices) {
appendLocalPointVaryingStencils(
varyingIndices, matrix, sourcePoints, sourcePointOffset);
if (applyVertexStencils) {
appendLocalPointStencils(
matrix, sourcePoints, sourcePointOffset);
if (applyVaryingStencils) {
for (int i = 0; i < numPatchPoints; ++i) {
appendLocalPointVaryingStencil(
varyingIndices, i, sourcePoints, sourcePointOffset);
}
}
}
for (int i = 0; i < numPatchPoints; ++i) {
patchPoints[i] = nextNewLocalPoint ++;
@ -1612,12 +1607,13 @@ PatchTableBuilder::LocalPointHelper::AppendLocalPatchPoints(
+ sourcePointOffset;
continue;
}
appendLocalPointStencil(
if (applyVertexStencils) {
appendLocalPointStencil(
matrix, i, sourcePoints, sourcePointOffset);
if (varyingIndices) {
appendLocalPointVaryingStencil(
varyingIndices, i, sourcePoints, sourcePointOffset);
if (applyVaryingStencils) {
appendLocalPointVaryingStencil(
varyingIndices, i, sourcePoints, sourcePointOffset);
}
}
patchPoints[i] = nextNewLocalPoint ++;
}
@ -1682,10 +1678,13 @@ PatchTableBuilder::LocalPointHelper::AppendLocalPatchPoints(
}
}
if (patchPoint == nextNewLocalPoint) {
appendLocalPointStencil(matrix, i, sourcePoints, sourcePointOffset);
if (varyingIndices) {
appendLocalPointVaryingStencil(
varyingIndices, i, sourcePoints, sourcePointOffset);
if (applyVertexStencils) {
appendLocalPointStencil(
matrix, i, sourcePoints, sourcePointOffset);
if (applyVaryingStencils) {
appendLocalPointVaryingStencil(
varyingIndices, i, sourcePoints, sourcePointOffset);
}
}
nextNewLocalPoint ++;
}

26
opensubdiv/sdc/catmarkScheme.h
View File

@ -121,7 +121,7 @@ Scheme<SCHEME_CATMARK>::assignSmoothMaskForEdge(EDGE const& edge, MASK& mask) co
//
// This mimics the implementation in Hbr in terms of order of operations.
//
const Weight CATMARK_SMOOTH_TRI_EDGE_WEIGHT = 0.470f;
const Weight CATMARK_SMOOTH_TRI_EDGE_WEIGHT = (Weight) 0.470;
Weight f0Weight = face0IsTri ? CATMARK_SMOOTH_TRI_EDGE_WEIGHT : 0.25f;
Weight f1Weight = face1IsTri ? CATMARK_SMOOTH_TRI_EDGE_WEIGHT : 0.25f;
@ -230,8 +230,8 @@ Scheme<SCHEME_CATMARK>::assignCreaseLimitMask(VERTEX const& vertex, MASK& posMas
posMask.SetNumFaceWeights(0);
posMask.SetFaceWeightsForFaceCenters(false);
Weight vWeight = 2.0f / 3.0f;
Weight eWeight = 1.0f / 6.0f;
Weight vWeight = (Weight)(2.0 / 3.0);
Weight eWeight = (Weight)(1.0 / 6.0);
posMask.VertexWeight(0) = vWeight;
for (int i = 0; i < valence; ++i) {
@ -261,9 +261,9 @@ Scheme<SCHEME_CATMARK>::assignSmoothLimitMask(VERTEX const& vertex, MASK& posMas
// Specialize for the regular case:
if (valence == 4) {
Weight fWeight = 1.0f / 36.0f;
Weight eWeight = 1.0f / 9.0f;
Weight vWeight = 4.0f / 9.0f;
Weight fWeight = (Weight)(1.0 / 36.0);
Weight eWeight = (Weight)(1.0 / 9.0);
Weight vWeight = (Weight)(4.0 / 9.0);
posMask.VertexWeight(0) = vWeight;
@ -385,14 +385,14 @@ Scheme<SCHEME_CATMARK>::assignCreaseLimitTangentMasks(VERTEX const& vertex,
if (interiorEdgeCount == 1) {
// The regular case -- uniform B-spline cross-tangent:
tan2Mask.VertexWeight(0) = -4.0f / 6.0f;
tan2Mask.VertexWeight(0) = (Weight)(-4.0 / 6.0);
tan2Mask.EdgeWeight(creaseEnds[0]) = -1.0f / 6.0f;
tan2Mask.EdgeWeight(creaseEnds[0] + 1) = 4.0f / 6.0f;
tan2Mask.EdgeWeight(creaseEnds[1]) = -1.0f / 6.0f;
tan2Mask.EdgeWeight(creaseEnds[0]) = (Weight)(-1.0 / 6.0);
tan2Mask.EdgeWeight(creaseEnds[0] + 1) = (Weight)( 4.0 / 6.0);
tan2Mask.EdgeWeight(creaseEnds[1]) = (Weight)(-1.0 / 6.0);
tan2Mask.FaceWeight(creaseEnds[0]) = 1.0f / 6.0f;
tan2Mask.FaceWeight(creaseEnds[0] + 1) = 1.0f / 6.0f;
tan2Mask.FaceWeight(creaseEnds[0]) = (Weight)(1.0 / 6.0);
tan2Mask.FaceWeight(creaseEnds[0] + 1) = (Weight)(1.0 / 6.0);
} else if (interiorEdgeCount > 1) {
// The irregular case -- formulae from Biermann et al:
@ -483,7 +483,7 @@ Scheme<SCHEME_CATMARK>::assignSmoothLimitTangentMasks(VERTEX const& vertex,
double cosTheta = std::cos(theta);
double cosHalfTheta = std::cos(theta * 0.5f);
double lambda = (5.0f / 16.0f) + (1.0f / 16.0f) *
double lambda = (5.0 / 16.0) + (1.0 / 16.0) *
(cosTheta + cosHalfTheta * std::sqrt(2.0f * (9.0f + cosTheta)));
double edgeWeightScale = 4.0f;

47
opensubdiv/sdc/loopScheme.h
View File

@ -29,6 +29,7 @@
#include "../sdc/scheme.h"
#include <cassert>
#include <cmath>
namespace OpenSubdiv {
namespace OPENSUBDIV_VERSION {
@ -194,13 +195,16 @@ Scheme<SCHEME_LOOP>::assignSmoothMaskForVertex(VERTEX const& vertex, MASK& mask)
if (valence != 6) {
// From HbrLoopSubdivision<T>::Subdivide(mesh, vertex):
// - could use some lookup tables here for common irregular valence (5, 7, 8)
// or all of these cosf() calls will be adding up...
// or all of these cosine calls will be adding up...
Weight invValence = 1.0f / (Weight) valence;
Weight beta = 0.25f * cosf((Weight)M_PI * 2.0f * invValence) + 0.375f;
double dValence = (double) valence;
double invValence = 1.0f / dValence;
double cosTheta = std::cos(M_PI * 2.0f * invValence);
eWeight = (0.625f - (beta * beta)) * invValence;;
vWeight = 1.0f - (eWeight * (Weight)valence);
double beta = 0.25f * cosTheta + 0.375f;
eWeight = (Weight) ((0.625f - (beta * beta)) * invValence);
vWeight = (Weight) (1.0f - (eWeight * dValence));
}
mask.VertexWeight(0) = vWeight;
@ -252,8 +256,8 @@ Scheme<SCHEME_LOOP>::assignCreaseLimitMask(VERTEX const& vertex, MASK& posMask,
// is based on an alternate refinement mask for the edge -- (3/8, 5/8) versus
// the usual (1/2, 1/2) -- and will not produce the B-spline curve desired.
//
Weight vWeight = 4.0f / 6.0f;
Weight eWeight = 1.0f / 6.0f;
Weight vWeight = (Weight) (4.0 / 6.0);
Weight eWeight = (Weight) (1.0 / 6.0);
posMask.VertexWeight(0) = vWeight;
for (int i = 0; i < valence; ++i) {
@ -271,7 +275,6 @@ Scheme<SCHEME_LOOP>::assignSmoothLimitMask(VERTEX const& vertex, MASK& posMask)
typedef typename MASK::Weight Weight;
int valence = vertex.GetNumFaces();
assert(valence != 2);
posMask.SetNumVertexWeights(1);
posMask.SetNumEdgeWeights(valence);
@ -280,7 +283,7 @@ Scheme<SCHEME_LOOP>::assignSmoothLimitMask(VERTEX const& vertex, MASK& posMask)
// Specialize for the regular case: 1/12 per edge-vert, 1/2 for the vert itself:
if (valence == 6) {
Weight eWeight = 1.0f / 12.0f;
Weight eWeight = (Weight) (1.0 / 12.0);
Weight vWeight = 0.5f;
posMask.VertexWeight(0) = vWeight;
@ -293,13 +296,15 @@ Scheme<SCHEME_LOOP>::assignSmoothLimitMask(VERTEX const& vertex, MASK& posMask)
posMask.EdgeWeight(5) = eWeight;
} else {
Weight invValence = 1.0f / valence;
double dValence = (double) valence;
double invValence = 1.0f / dValence;
double cosTheta = std::cos(M_PI * 2.0f * invValence);
Weight beta = 0.25f * cosf((Weight)M_PI * 2.0f * invValence) + 0.375f;
beta = (0.625f - (beta * beta)) * invValence;;
double beta = 0.25f * cosTheta + 0.375f;
double gamma = (0.625f - (beta * beta)) * invValence;
Weight eWeight = 1.0f / (valence + 3.0f / (8.0f * beta));
Weight vWeight = (Weight)(1.0f - (eWeight * valence));
Weight eWeight = (Weight) (1.0f / (dValence + 3.0f / (8.0f * gamma)));
Weight vWeight = (Weight) (1.0f - (eWeight * dValence));
posMask.VertexWeight(0) = vWeight;
for (int i = 0; i < valence; ++i) {
@ -472,7 +477,7 @@ Scheme<SCHEME_LOOP>::assignCreaseLimitTangentMasks(VERTEX const& vertex,
if (interiorEdgeCount == 2) {
// See note above regarding scale factor of (sin(60 degs) == sqrt(3)/2:
static Weight const Root3 = (Weight) 1.73205080756887729352f;
static Weight const Root3 = (Weight) 1.73205080756887729352;
static Weight const Root3by2 = (Weight) (Root3 * 0.5);
tan2Mask.VertexWeight(0) = -Root3;
@ -489,16 +494,15 @@ Scheme<SCHEME_LOOP>::assignCreaseLimitTangentMasks(VERTEX const& vertex,
double theta = M_PI / (interiorEdgeCount + 1);
Weight cWeight = -3.0f * (Weight) std::sin(theta);
Weight eWeightCoeff = -3.0f * (2.0f * (Weight) std::cos(theta) - 2.0f);
tan2Mask.VertexWeight(0) = 0.0f;
Weight cWeight = (Weight) (-3.0f * std::sin(theta));
tan2Mask.EdgeWeight(creaseEnds[0]) = cWeight;
tan2Mask.EdgeWeight(creaseEnds[1]) = cWeight;
double eCoeff = -3.0f * 2.0f * (std::cos(theta) - 1.0f);
for (int i = 1; i <= interiorEdgeCount; ++i) {
tan2Mask.EdgeWeight(creaseEnds[0] + i) = eWeightCoeff * (Weight) std::sin(i * theta);
tan2Mask.EdgeWeight(creaseEnds[0] + i) = (Weight) (eCoeff * std::sin(i * theta));
}
} else if (interiorEdgeCount == 1) {
// See notes above regarding scale factor of 3.0:
@ -531,7 +535,6 @@ Scheme<SCHEME_LOOP>::assignSmoothLimitTangentMasks(VERTEX const& vertex,
typedef typename MASK::Weight Weight;
int valence = vertex.GetNumFaces();
assert(valence != 2);
tan1Mask.SetNumVertexWeights(1);
tan1Mask.SetNumEdgeWeights(valence);
@ -547,7 +550,7 @@ Scheme<SCHEME_LOOP>::assignSmoothLimitTangentMasks(VERTEX const& vertex,
tan2Mask.VertexWeight(0) = 0.0f;
if (valence == 6) {
static Weight const Root3by2 = (Weight)(0.5f * 1.73205080756887729352f);
static Weight const Root3by2 = (Weight)(0.5 * 1.73205080756887729352);
tan1Mask.EdgeWeight(0) = 1.0f;
tan1Mask.EdgeWeight(1) = 0.5f;
@ -563,7 +566,7 @@ Scheme<SCHEME_LOOP>::assignSmoothLimitTangentMasks(VERTEX const& vertex,
tan2Mask.EdgeWeight(4) = -Root3by2;
tan2Mask.EdgeWeight(5) = -Root3by2;
} else {
Weight alpha = (Weight) (2.0f * M_PI / valence);
double alpha = 2.0f * M_PI / valence;
for (int i = 0; i < valence; ++i) {
double alphaI = alpha * i;
tan1Mask.EdgeWeight(i) = (Weight) std::cos(alphaI);