AuroraMiddleware/OpenSubdiv
1
0
Fork 0
mirror of https://github.com/PixarAnimationStudios/OpenSubdiv synced 2025-01-08 07:40:17 +00:00
Code Issues Releases Wiki Activity

Merge pull request #991 from barfowl/val2_gregory_simplify

Browse Source 
Simplify construction of Gregory patches with val-2 interior corners
This commit is contained in:
David G Yu 2018-09-13 12:29:06 -07:00 committed by GitHub
commit 79476c87cb
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 254 additions and 239 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

493
opensubdiv/far/catmarkPatchBuilder.cpp
View File

@ -310,88 +310,44 @@ CatmarkLimits<REAL>::ComputeBoundaryPointWeights(int valence, int faceInRing,
}
//
// SparseMatrixPoint
// SparseMatrixRow
//
// This is a utility class representing a row of a SparseMatrix -- which
// in turn corresponds to a point of a resulting patch.
// in turn corresponds to a point of a resulting patch. Instances of this
// class are intended to encapsulate the contributions of a point and be
// passed to functions as such.
//
// This interface was originally transitional (supporting a migration away
// from the former GregoryBasis::Point class) and its unclear if it will
// persist. Its needs have been simplified and given the usual pre-sizing
// of a sparse row, a simple Assign() method may be all that is necessary.
// (Consider moving this to PatchBuilder as a protected class or maybe a
// public class within SparseMatrix itself, e.g. SparseMatrix<REAL>::Row.)
//
// Use of the AddOrAppend() method is highly discouraged as it requires
// iteration and testing of potentially all entries, before appending if
// the specified entry does not exist. It is currently only used to
// simplify the awkward case of valence-2 interior vertices -- where the
// 1-ring of neighboring vertices overlaps and has duplicate entries.
// Alternative to simplify this without resorting to AddOrAppend() are
// under consideration.
//
template <typename REAL>
class SparseMatrixPoint {
public:
typedef Index index_type;
typedef REAL weight_type;
namespace {
template <typename REAL>
class SparseMatrixRow {
public:
SparseMatrixRow(SparseMatrix<REAL> & matrix, int row) :
_size(matrix.GetRowSize(row)),
_indices(matrix.SetRowColumns(row).begin()),
_weights(matrix.SetRowElements(row).begin()) { }
typedef SparseMatrix<weight_type> matrix_type;
public:
SparseMatrixPoint(matrix_type & matrix, int row, int size = -1);
int GetSize() const { return _size; }
int GetSize() const { return _size; }
int GetCapacity() const { return _indices.size(); }
void Append( index_type index, weight_type weight);
void AddOrAppend( index_type index, weight_type weight);
void Assign(int rowEntry, index_type index, weight_type weight);
void Copy(SparseMatrixPoint<weight_type> const & other);
public:
int _size;
Array<index_type> _indices;
Array<weight_type> _weights;
};
template <typename REAL>
inline
SparseMatrixPoint<REAL>::SparseMatrixPoint(matrix_type & matrix, int row, int size) {
_indices = matrix.SetRowColumns(row);
_weights = matrix.SetRowElements(row);
_size = (size < 0) ? _weights.size() : size;
}
template <typename REAL>
inline void
SparseMatrixPoint<REAL>::Assign(int rowEntry, index_type index, weight_type weight) {
_indices[rowEntry] = index;
_weights[rowEntry] = weight;
}
template <typename REAL>
inline void
SparseMatrixPoint<REAL>::Append(index_type index, weight_type weight) {
assert(GetSize() < GetCapacity());
_indices[_size] = index;
_weights[_size] = weight;
_size ++;
}
template <typename REAL>
inline void
SparseMatrixPoint<REAL>::AddOrAppend(index_type index, weight_type weight) {
for (int i = 0; i < GetSize(); ++i) {
if (_indices[i] == index) {
_weights[i] += weight;
return;
void Assign(int rowEntry, Index index, REAL weight) {
_indices[rowEntry] = index;
_weights[rowEntry] = weight;
}
}
Append(index, weight);
}
template <typename REAL>
inline void
SparseMatrixPoint<REAL>::Copy(SparseMatrixPoint const & other) {
assert(GetCapacity() == other.GetCapacity());
_size = other._size;
std::memcpy(&_indices[0], &other._indices[0], _size * sizeof(index_type));
std::memcpy(&_weights[0], &other._weights[0], _size * sizeof(weight_type));
}
void Copy(SparseMatrixRow<REAL> const & other) {
assert(GetSize() == other.GetSize());
std::memcpy(_indices, other._indices, _size * sizeof(Index));
std::memcpy(_weights, other._weights, _size * sizeof(REAL));
}
public:
int _size;
Index * _indices;
REAL * _weights;
};
} // end namespace
//
@ -466,7 +422,7 @@ namespace {
template <typename REAL>
void
_addSparsePointToFullRow(REAL * fullRow,
SparseMatrixPoint<REAL> const & p,
SparseMatrixRow<REAL> const & p,
REAL s, int * indexMask) {
for (int i = 0; i < p.GetSize(); ++i) {
@ -524,6 +480,85 @@ namespace {
prefix, M.GetNumRows(), M.GetNumColumns(), fullSize,
sparseSize, nonZeroSize, (REAL)nonZeroSize * 100.0f / (REAL)fullSize);
}
//
// The valence-2 interior case poses problems for the way patch points
// are computed as combinations of source points and stored as a row in
// a SparseMatrix. An interior vertex of valence-2 causes duplicate
// vertices to appear in the 1-rings of its neighboring vertices and we
// want the entries of a SparseMatrix row to be unique.
//
// For the most part, this does not pose a problem while the matrix (set
// of patch points) is being constructed, so we leave those duplicate
// entries in place and deal with them as a post-process here.
//
// The SourcePatch is also sensitive to the presence of such valence-2
// vertices for its own reasons (it needs to identifiy a unique set of
// source points from a set of corner rings), so a simple query of its
// corners indicates when this post-process is necessary. (And since
// this case is a rare occurrence, efficiency is not a major concern.)
//
template <typename REAL>
void _removeValence2Duplicates(SparseMatrix<REAL> & M) {
// This will later be determined by the PatchBuilder member:
int regFaceSize = 4;
SparseMatrix<REAL> T;
T.Resize(M.GetNumRows(), M.GetNumColumns(), M.GetNumElements());
int nRows = M.GetNumRows();
for (int row = 0; row < nRows; ++row) {
int srcRowSize = M.GetRowSize(row);
int const * srcIndices = M.GetRowColumns(row).begin();
REAL const * srcWeights = M.GetRowElements(row).begin();
// Scan the entries to see if there are duplicates -- copy
// the row if not, otherwise, need to compress it:
bool cornerUsed[4] = { false, false, false, false };
int srcDupCount = 0;
for (int i = 0; i < srcRowSize; ++i) {
int srcIndex = srcIndices[i];
if (srcIndex < regFaceSize) {
srcDupCount += (int) cornerUsed[srcIndex];
cornerUsed[srcIndex] = true;
}
}
// Size this row for the destination and copy or compress:
T.SetRowSize(row, srcRowSize - srcDupCount);
int* dstIndices = T.SetRowColumns(row).begin();
REAL* dstWeights = T.SetRowElements(row).begin();
if (srcDupCount) {
REAL * cornerDstPtr[4] = { 0, 0, 0, 0 };
for (int i = 0; i < srcRowSize; ++i) {
int srcIndex = *srcIndices++;
REAL srcWeight = *srcWeights++;
if (srcIndex < regFaceSize) {
if (cornerDstPtr[srcIndex]) {
*cornerDstPtr[srcIndex] += srcWeight;
continue;
} else {
cornerDstPtr[srcIndex] = dstWeights;
}
}
*dstIndices++ = srcIndex;
*dstWeights++ = srcWeight;
}
} else {
std::memcpy(&dstIndices[0], &srcIndices[0], srcRowSize * sizeof(Index));
std::memcpy(&dstWeights[0], &srcWeights[0], srcRowSize * sizeof(REAL));
}
}
M.Swap(T);
}
} // end namespace for SparseMatrix utilities
@ -568,7 +603,7 @@ class GregoryConverter {
public:
typedef REAL Weight;
typedef SparseMatrix<Weight> Matrix;
typedef SparseMatrixPoint<Weight> Point;
typedef SparseMatrixRow<Weight> Point;
public:
GregoryConverter() : _numSourcePoints(0) { }
GregoryConverter(SourcePatch const & sourcePatch);
@ -577,6 +612,7 @@ public:
void Initialize(SourcePatch const & sourcePatch);
bool IsIsolatedInteriorPatch() const { return _isIsolatedInteriorPatch; }
bool HasVal2InteriorCorner() const { return _hasVal2InteriorCorner; }
int GetIsolatedInteriorCorner() const { return _isolatedCorner; }
int GetIsolatedInteriorValence() const { return _isolatedValence; }
@ -605,10 +641,6 @@ private:
unsigned int fpIsCopied : 1;
unsigned int fmIsCopied : 1;
unsigned int val2InRing : 1;
unsigned int epAdjToVal2 : 1;
unsigned int emAdjToVal2 : 1;
// Other values stored for repeated use:
int valence;
int numFaces;
@ -660,6 +692,7 @@ private:
int _maxValence;
bool _isIsolatedInteriorPatch;
bool _hasVal2InteriorCorner;
int _isolatedCorner;
int _isolatedValence;
@ -705,6 +738,7 @@ GregoryConverter<REAL>::Initialize(SourcePatch const & sourcePatch) {
int irregularCorner = -1;
int irregularValence = -1;
int sharpCount = 0;
int val2IntCount = 0;
for (int cIndex = 0; cIndex < 4; ++cIndex) {
SourcePatch::Corner srcCorner = sourcePatch._corners[cIndex];
@ -717,7 +751,6 @@ GregoryConverter<REAL>::Initialize(SourcePatch const & sourcePatch) {
corner.numFaces = srcCorner._numFaces;
corner.faceInRing = srcCorner._patchFace;
corner.isVal2Int = srcCorner._val2Interior;
corner.val2InRing = srcCorner._val2Adjacent;
corner.valence = corner.numFaces + corner.isBoundary;
corner.isRegular = ((corner.numFaces << corner.isBoundary) == 4)
@ -745,6 +778,7 @@ GregoryConverter<REAL>::Initialize(SourcePatch const & sourcePatch) {
irregularValence = corner.valence;
}
sharpCount += corner.isSharp;
val2IntCount += corner.isVal2Int;
}
// Make a second pass to assign tags dependent on adjacent corners
@ -754,9 +788,6 @@ GregoryConverter<REAL>::Initialize(SourcePatch const & sourcePatch) {
int cNext = (cIndex + 1) & 0x3;
int cPrev = (cIndex + 3) & 0x3;
corner.epAdjToVal2 = _corners[cNext].isVal2Int;
corner.emAdjToVal2 = _corners[cPrev].isVal2Int;
//
// Identify if the face points are regular or shared/copied from
// one of the pair:
@ -793,6 +824,7 @@ GregoryConverter<REAL>::Initialize(SourcePatch const & sourcePatch) {
_isolatedCorner = irregularCorner;
_isolatedValence = irregularValence;
}
_hasVal2InteriorCorner = (val2IntCount > 0);
}
template <typename REAL>
@ -838,6 +870,9 @@ GregoryConverter<REAL>::Convert(Matrix & matrix) const {
computeIrregularFacePoints(cIndex, matrix, weightBuffer, indexBuffer);
}
}
if (_hasVal2InteriorCorner) {
_removeValence2Duplicates(matrix);
}
}
template <typename REAL>
@ -904,9 +939,9 @@ GregoryConverter<REAL>::resizeMatrixUnisolated(Matrix & matrix) const {
// First, the corner and pair of edge points:
if (corner.isRegular) {
if (! corner.isBoundary) {
rowSize[0] = 9 - corner.val2InRing;
rowSize[1] = 6 - corner.epAdjToVal2;
rowSize[2] = 6 - corner.emAdjToVal2;
rowSize[0] = 9;
rowSize[1] = 6;
rowSize[2] = 6;
} else {
rowSize[0] = 3;
rowSize[1] = corner.faceInRing ? 6 : 2;
@ -918,12 +953,12 @@ GregoryConverter<REAL>::resizeMatrixUnisolated(Matrix & matrix) const {
rowSize[1] = 2;
rowSize[2] = 2;
} else if (! corner.isBoundary) {
int ringSize = 1 + 2 * corner.valence - corner.val2InRing;
int ringSize = 1 + 2 * corner.valence;
rowSize[0] = ringSize;
rowSize[1] = ringSize;
rowSize[2] = ringSize;
} else if (corner.numFaces > 1) {
int ringSize = 1 + corner.valence + corner.numFaces - corner.val2InRing;
int ringSize = 1 + corner.valence + corner.numFaces;
rowSize[0] = 3;
rowSize[1] = (corner.faceInRing > 0) ? ringSize : 2;
rowSize[2] = (corner.faceInRing < (corner.numFaces - 1)) ? ringSize : 2;
@ -960,31 +995,25 @@ template <typename REAL>
void
GregoryConverter<REAL>::assignRegularEdgePoints(int cIndex, Matrix & matrix) const {
// Declare with 0 size for use with Append()
Point p (matrix, 5*cIndex + 0, 0);
Point ep(matrix, 5*cIndex + 1, 0);
Point em(matrix, 5*cIndex + 2, 0);
Point p (matrix, 5*cIndex + 0);
Point ep(matrix, 5*cIndex + 1);
Point em(matrix, 5*cIndex + 2);
CornerTopology const & corner = _corners[cIndex];
int const * cRing = corner.ringPoints;
if (! corner.isBoundary) {
p.Append(cIndex, (REAL) (4.0 / 9.0));
p.Append(cRing[0], (REAL) (1.0 / 9.0));
p.Append(cRing[2], (REAL) (1.0 / 9.0));
p.Append(cRing[4], (REAL) (1.0 / 9.0));
p.Append(cRing[6], (REAL) (1.0 / 9.0));
p.Append(cRing[1], (REAL) (1.0 / 36.0));
if (!corner.val2InRing) {
p.Append(cRing[3], (REAL) (1.0 / 36.0));
p.Append(cRing[5], (REAL) (1.0 / 36.0));
p.Append(cRing[7], (REAL) (1.0 / 36.0));
} else {
p.AddOrAppend(cRing[3], (REAL) (1.0 / 36.0));
p.AddOrAppend(cRing[5], (REAL) (1.0 / 36.0));
p.AddOrAppend(cRing[7], (REAL) (1.0 / 36.0));
}
p.Assign(0, cIndex, (REAL) (4.0 / 9.0));
p.Assign(1, cRing[0], (REAL) (1.0 / 9.0));
p.Assign(2, cRing[2], (REAL) (1.0 / 9.0));
p.Assign(3, cRing[4], (REAL) (1.0 / 9.0));
p.Assign(4, cRing[6], (REAL) (1.0 / 9.0));
p.Assign(5, cRing[1], (REAL) (1.0 / 36.0));
p.Assign(6, cRing[3], (REAL) (1.0 / 36.0));
p.Assign(7, cRing[5], (REAL) (1.0 / 36.0));
p.Assign(8, cRing[7], (REAL) (1.0 / 36.0));
assert(p.GetSize() == 9);
// Identify the edges along Ep and Em and those opposite them:
int iEdgeEp = 2 * corner.faceInRing;
@ -992,50 +1021,44 @@ GregoryConverter<REAL>::assignRegularEdgePoints(int cIndex, Matrix & matrix) con
int iEdgeOp = 2 * ((corner.faceInRing + 2) & 0x3);
int iEdgeOm = 2 * ((corner.faceInRing + 3) & 0x3);
ep.Append(cIndex, (REAL) (4.0 / 9.0));
ep.Append(cRing[iEdgeEp], (REAL) (2.0 / 9.0));
ep.Append(cRing[iEdgeEm], (REAL) (1.0 / 9.0));
ep.Append(cRing[iEdgeOm], (REAL) (1.0 / 9.0));
ep.Append(cRing[iEdgeEp + 1], (REAL) (1.0 / 18.0));
if (corner.epAdjToVal2) {
ep.AddOrAppend(cRing[iEdgeOm + 1], (REAL) (1.0 / 18.0));
} else {
ep.Append(cRing[iEdgeOm + 1], (REAL) (1.0 / 18.0));
}
ep.Assign(0, cIndex, (REAL) (4.0 / 9.0));
ep.Assign(1, cRing[iEdgeEp], (REAL) (2.0 / 9.0));
ep.Assign(2, cRing[iEdgeEm], (REAL) (1.0 / 9.0));
ep.Assign(3, cRing[iEdgeOm], (REAL) (1.0 / 9.0));
ep.Assign(4, cRing[iEdgeEp + 1], (REAL) (1.0 / 18.0));
ep.Assign(5, cRing[iEdgeOm + 1], (REAL) (1.0 / 18.0));
assert(ep.GetSize() == 6);
em.Append(cIndex, (REAL) (4.0 / 9.0));
em.Append(cRing[iEdgeEm], (REAL) (2.0 / 9.0));
em.Append(cRing[iEdgeEp], (REAL) (1.0 / 9.0));
em.Append(cRing[iEdgeOp], (REAL) (1.0 / 9.0));
em.Append(cRing[iEdgeEp + 1], (REAL) (1.0 / 18.0));
if (corner.emAdjToVal2) {
em.AddOrAppend(cRing[iEdgeEm + 1], (REAL) (1.0 / 18.0));
} else {
em.Append(cRing[iEdgeEm + 1], (REAL) (1.0 / 18.0));
}
em.Assign(0, cIndex, (REAL) (4.0 / 9.0));
em.Assign(1, cRing[iEdgeEm], (REAL) (2.0 / 9.0));
em.Assign(2, cRing[iEdgeEp], (REAL) (1.0 / 9.0));
em.Assign(3, cRing[iEdgeOp], (REAL) (1.0 / 9.0));
em.Assign(4, cRing[iEdgeEp + 1], (REAL) (1.0 / 18.0));
em.Assign(5, cRing[iEdgeEm + 1], (REAL) (1.0 / 18.0));
assert(em.GetSize() == 6);
} else {
// Decide which point corresponds to interior vs exterior tangent:
Point & eInterior = corner.faceInRing ? ep : em;
Point & eBoundary = corner.faceInRing ? em : ep;
int iBoundary = corner.faceInRing ? 4 : 0;
p.Append(cIndex, (REAL) (2.0 / 3.0));
p.Append(cRing[0], (REAL) (1.0 / 6.0));
p.Append(cRing[4], (REAL) (1.0 / 6.0));
p.Assign(0, cIndex, (REAL) (2.0 / 3.0));
p.Assign(1, cRing[0], (REAL) (1.0 / 6.0));
p.Assign(2, cRing[4], (REAL) (1.0 / 6.0));
assert(p.GetSize() == 3);
eBoundary.Append(cIndex, (REAL) (2.0 / 3.0));
eBoundary.Append(cRing[iBoundary], (REAL) (1.0 / 3.0));
eBoundary.Assign(0, cIndex, (REAL) (2.0 / 3.0));
eBoundary.Assign(1, cRing[iBoundary], (REAL) (1.0 / 3.0));
assert(eBoundary.GetSize() == 2);
eInterior.Append(cIndex, (REAL) (4.0 / 9.0));
eInterior.Append(cRing[2], (REAL) (2.0 / 9.0));
eInterior.Append(cRing[0], (REAL) (1.0 / 9.0));
eInterior.Append(cRing[4], (REAL) (1.0 / 9.0));
eInterior.Append(cRing[1], (REAL) (1.0 / 18.0));
eInterior.Append(cRing[3], (REAL) (1.0 / 18.0));
eInterior.Assign(0, cIndex, (REAL) (4.0 / 9.0));
eInterior.Assign(1, cRing[2], (REAL) (2.0 / 9.0));
eInterior.Assign(2, cRing[0], (REAL) (1.0 / 9.0));
eInterior.Assign(3, cRing[4], (REAL) (1.0 / 9.0));
eInterior.Assign(4, cRing[1], (REAL) (1.0 / 18.0));
eInterior.Assign(5, cRing[3], (REAL) (1.0 / 18.0));
assert(eInterior.GetSize() == 6);
}
assert(matrix.GetRowSize(5*cIndex + 0) == p.GetSize());
assert(matrix.GetRowSize(5*cIndex + 1) == ep.GetSize());
assert(matrix.GetRowSize(5*cIndex + 2) == em.GetSize());
}
template <typename REAL>
@ -1043,10 +1066,9 @@ void
GregoryConverter<REAL>::computeIrregularEdgePoints(int cIndex,
Matrix & matrix, Weight *weightBuffer) const {
// Declare with 0 size for use with Append()
Point p (matrix, 5*cIndex + 0, 0);
Point ep(matrix, 5*cIndex + 1, 0);
Point em(matrix, 5*cIndex + 2, 0);
Point p (matrix, 5*cIndex + 0);
Point ep(matrix, 5*cIndex + 1);
Point em(matrix, 5*cIndex + 2);
//
// The corner and edge points P, Ep and Em are completely determined
@ -1059,14 +1081,17 @@ GregoryConverter<REAL>::computeIrregularEdgePoints(int cIndex,
//
// The sharp case -- both interior and boundary...
//
p.Append(cIndex, 1.0f);
p.Assign(0, cIndex, 1.0f);
assert(p.GetSize() == 1);
// Approximating these for now, pending future investigation...
ep.Append(cIndex, (REAL)(2.0 / 3.0));
ep.Append((cIndex+1) & 0x3, (REAL)(1.0 / 3.0));
ep.Assign(0, cIndex, (REAL)(2.0 / 3.0));
ep.Assign(1, (cIndex+1) & 0x3, (REAL)(1.0 / 3.0));
assert(ep.GetSize() == 2);
em.Append(cIndex, (REAL)(2.0 / 3.0));
em.Append((cIndex+3) & 0x3, (REAL)(1.0 / 3.0));
em.Assign(0, cIndex, (REAL)(2.0 / 3.0));
em.Assign(1, (cIndex+3) & 0x3, (REAL)(1.0 / 3.0));
assert(em.GetSize() == 2);
} else if (! corner.isBoundary) {
//
// The irregular interior case:
@ -1081,20 +1106,19 @@ GregoryConverter<REAL>::computeIrregularEdgePoints(int cIndex,
//
// The irregular/smooth corner case:
//
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));
p.Assign(0, cIndex, (REAL)(4.0 / 6.0));
p.Assign(1, (cIndex+1) & 0x3, (REAL)(1.0 / 6.0));
p.Assign(2, (cIndex+3) & 0x3, (REAL)(1.0 / 6.0));
assert(p.GetSize() == 3);
ep.Append(cIndex, (REAL)(2.0 / 3.0));
ep.Append((cIndex+1) & 0x3, (REAL)(1.0 / 3.0));
ep.Assign(0, cIndex, (REAL)(2.0 / 3.0));
ep.Assign(1, (cIndex+1) & 0x3, (REAL)(1.0 / 3.0));
assert(ep.GetSize() == 2);
em.Append(cIndex, (REAL)(2.0 / 3.0));
em.Append((cIndex+3) & 0x3, (REAL)(1.0 / 3.0));
em.Assign(0, cIndex, (REAL)(2.0 / 3.0));
em.Assign(1, (cIndex+3) & 0x3, (REAL)(1.0 / 3.0));
assert(em.GetSize() == 2);
}
assert(matrix.GetRowSize(5*cIndex + 0) == p.GetSize());
assert(matrix.GetRowSize(5*cIndex + 1) == ep.GetSize());
assert(matrix.GetRowSize(5*cIndex + 2) == em.GetSize());
}
@ -1127,36 +1151,20 @@ GregoryConverter<REAL>::computeIrregularInteriorEdgePoints(
// since Ep and Em depend on it, there should be no need to filter weights
// with value 0:
//
// The presence of overlapping contributions in the ring when a neighboring
// vertex is valence-2 interior is a nuisance and highly unlikely but still
// possible if we have a triangle near a valence-2 interior vertex -- they
// will still not be isolated at level 1.
//
// For now just use the much-less-efficient AddOrAppend() in this case
// until we can identify the overlap more specifically and have a more
// general mechanism availble to deal with it.
//
p.Append( cIndex, pWeights[0]);
ep.Append(cIndex, epWeights[0]);
em.Append(cIndex, emWeights[0]);
p.Assign( 0, cIndex, pWeights[0]);
ep.Assign(0, cIndex, epWeights[0]);
em.Assign(0, cIndex, emWeights[0]);
if (!corner.val2InRing) {
for (int i = 1; i < weightWidth; ++i) {
int pRingPoint = corner.ringPoints[i-1];
for (int i = 1; i < weightWidth; ++i) {
int pRingPoint = corner.ringPoints[i-1];
p.Append( pRingPoint, pWeights[i]);
ep.Append(pRingPoint, epWeights[i]);
em.Append(pRingPoint, emWeights[i]);
}
} else {
for (int i = 1; i < weightWidth; ++i) {
int pRingPoint = corner.ringPoints[i-1];
p.AddOrAppend( pRingPoint, pWeights[i]);
ep.AddOrAppend(pRingPoint, epWeights[i]);
em.AddOrAppend(pRingPoint, emWeights[i]);
}
p.Assign( i, pRingPoint, pWeights[i]);
ep.Assign(i, pRingPoint, epWeights[i]);
em.Assign(i, pRingPoint, emWeights[i]);
}
assert(p.GetSize() == weightWidth);
assert(ep.GetSize() == weightWidth);
assert(em.GetSize() == weightWidth);
}
@ -1193,29 +1201,36 @@ GregoryConverter<REAL>::computeIrregularBoundaryEdgePoints(
int p1 = corner.ringPoints[0];
int pN = corner.ringPoints[2*(valence-1)];
p.Append(p0, pWeights[0]);
p.Append(p1, pWeights[1]);
p.Append(pN, pWeights[N]);
p.Assign(0, p0, pWeights[0]);
p.Assign(1, p1, pWeights[1]);
p.Assign(2, pN, pWeights[N]);
assert(p.GetSize() == 3);
// If Ep is on the boundary edge, it has only two non-zero weights along
// that edge:
ep.Append(p0, epWeights[0]);
ep.Append(p1, epWeights[1]);
if (corner.faceInRing > 0) {
for (int i = 2; i < weightWidth; ++i) {
ep.Append(corner.ringPoints[i-1], epWeights[i]);
ep.Assign(0, p0, epWeights[0]);
if (corner.faceInRing == 0) {
ep.Assign(1, p1, epWeights[1]);
assert(ep.GetSize() == 2);
} else {
for (int i = 1; i < weightWidth; ++i) {
ep.Assign(i, corner.ringPoints[i-1], epWeights[i]);
}
assert(ep.GetSize() == weightWidth);
}
// If Em is on the boundary edge, it has only two non-zero weights along
// that edge:
em.Append(p0, emWeights[0]);
if (corner.faceInRing < (corner.numFaces - 1)) {
for (int i = 1; i < N; ++i) {
em.Append(corner.ringPoints[i-1], emWeights[i]);
em.Assign(0, p0, emWeights[0]);
if (corner.faceInRing == (corner.numFaces - 1)) {
em.Assign(1, pN, emWeights[N]);
assert(em.GetSize() == 2);
} else {
for (int i = 1; i <= weightWidth; ++i) {
em.Assign(i, corner.ringPoints[i-1], emWeights[i]);
}
assert(em.GetSize() == weightWidth);
}
em.Append(pN, emWeights[N]);
}
@ -1228,10 +1243,10 @@ GregoryConverter<REAL>::getIrregularFacePointSize(
CornerTopology const & adjCorner = _corners[cIndexFar];
int thisSize = corner.isSharp
? (6 - (corner.valence == 2))
: (1 + corner.ringPoints.GetSize() - corner.val2InRing);
? 6
: (1 + corner.ringPoints.GetSize());
int adjSize = (adjCorner.isRegular || adjCorner.isSharp || (adjCorner.valence == 2))
int adjSize = (adjCorner.isRegular || adjCorner.isSharp)
? 0
: (1 + adjCorner.ringPoints.GetSize() - 6);
@ -1268,9 +1283,6 @@ GregoryConverter<REAL>::computeIrregularFacePoint(
// Remember that R is to be computed about an interior edge and is
// comprised of the two pairs of points opposite the interior edge
//
// Remember also that val-2-overlap may cause two of these to be the
// same -- doesn't matter if we accumulate here will if we assign:
//
int iEdgeInterior = edgeInNearCornerRing;
int iEdgePrev = (iEdgeInterior + valence - 1) % valence;
int iEdgeNext = (iEdgeInterior + 1) % valence;
@ -1280,20 +1292,28 @@ GregoryConverter<REAL>::computeIrregularFacePoint(
rowWeights[cornerNear.ringPoints[2*iEdgeInterior + 1]] += signForSideOfEdge / 18.0f;
rowWeights[cornerNear.ringPoints[2*iEdgeNext]] += signForSideOfEdge / 9.0f;
int nWeights = 0;
for (int i = 0; i < fullRowSize; ++i) {
if (columnMask[i]) {
fNear.Append(columnMask[i] - 1, rowWeights[i]);
fNear.Assign(nWeights++, columnMask[i] - 1, rowWeights[i]);
}
}
// Complete the expected row size when val-2 interior corners induce duplicates:
if (_hasVal2InteriorCorner && (nWeights < fNear.GetSize())) {
while (nWeights < fNear.GetSize()) {
fNear.Assign(nWeights++, cIndexNear, 0.0f);
}
}
assert(fNear.GetSize() == nWeights);
}
template <typename REAL>
void
GregoryConverter<REAL>::assignRegularFacePoints(int cIndex, Matrix & matrix) const {
// Declare with 0 size for use with Append()
Point fp(matrix, 5*cIndex + 3, 0);
Point fm(matrix, 5*cIndex + 4, 0);
Point fp(matrix, 5*cIndex + 3);
Point fm(matrix, 5*cIndex + 4);
CornerTopology const & corner = _corners[cIndex];
@ -1303,18 +1323,18 @@ GregoryConverter<REAL>::assignRegularFacePoints(int cIndex, Matrix & matrix) con
// Assign regular Fp and/or Fm:
if (corner.fpIsRegular) {
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());
fp.Assign(0, cIndex, (REAL)(4.0 / 9.0));
fp.Assign(1, cPrev, (REAL)(2.0 / 9.0));
fp.Assign(2, cNext, (REAL)(2.0 / 9.0));
fp.Assign(3, cOpp, (REAL)(1.0 / 9.0));
assert(fp.GetSize() == 4);
}
if (corner.fmIsRegular) {
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());
fm.Assign(0, cIndex, (REAL)(4.0 / 9.0));
fm.Assign(1, cPrev, (REAL)(2.0 / 9.0));
fm.Assign(2, cNext, (REAL)(2.0 / 9.0));
fm.Assign(3, cOpp, (REAL)(1.0 / 9.0));
assert(fm.GetSize() == 4);
}
}
@ -1335,9 +1355,8 @@ GregoryConverter<REAL>::computeIrregularFacePoints(int cIndex,
Point ep (matrix, 5*cIndex + 1);
Point emNext(matrix, 5*cNext + 2);
// Declare with 0 size for use with Append()
Point fp(matrix, 5*cIndex + 3, 0);
Point fm(matrix, 5*cIndex + 4, 0);
Point fp(matrix, 5*cIndex + 3);
Point fm(matrix, 5*cIndex + 4);
//
// Compute the face points Fp and Fm in terms of the corner (P) and edge
@ -1795,7 +1814,6 @@ class LinearConverter {
public:
typedef REAL Weight;
typedef SparseMatrix<Weight> Matrix;
typedef SparseMatrixPoint<Weight> MatrixPoint;
public:
LinearConverter() : _sourcePatch(0) { }
LinearConverter(SourcePatch const & sourcePatch);
@ -1838,6 +1856,8 @@ LinearConverter<REAL>::Convert(Matrix & matrix) const {
matrix.Resize(4, _sourcePatch->GetNumSourcePoints(), numElements);
bool hasVal2InteriorCorner = false;
for (int cIndex = 0; cIndex < 4; ++cIndex) {
// Deal with the trivial sharp case first:
if (_sourcePatch->_corners[cIndex]._sharp) {
@ -1853,7 +1873,7 @@ LinearConverter<REAL>::Convert(Matrix & matrix) const {
if (sourceCorner._boundary) {
matrix.SetRowSize(cIndex, 3);
} else {
matrix.SetRowSize(cIndex, 1 + ringSize - sourceCorner._val2Adjacent);
matrix.SetRowSize(cIndex, 1 + ringSize);
}
Array<Index> rowIndices = matrix.SetRowColumns(cIndex);
@ -1879,18 +1899,13 @@ LinearConverter<REAL>::Convert(Matrix & matrix) const {
sourceCorner._numFaces, sourceCorner._patchFace,
&weightBuffer[0], 0, 0);
if (sourceCorner._val2Adjacent) {
MatrixPoint row(matrix, cIndex, 0);
for (int i = 0; i <= ringSize; ++i) {
row.AddOrAppend(indexBuffer[i], weightBuffer[i]);
}
} else {
std::memcpy(
&rowIndices[0], indexBuffer, (1 + ringSize) * sizeof(Index));
std::memcpy(
&rowWeights[0], weightBuffer, (1 + ringSize) * sizeof(Weight));
}
std::memcpy(&rowIndices[0], indexBuffer, (1 + ringSize) * sizeof(Index));
std::memcpy(&rowWeights[0], weightBuffer, (1 + ringSize) * sizeof(Weight));
}
hasVal2InteriorCorner |= sourceCorner._val2Interior;
}
if (hasVal2InteriorCorner) {
_removeValence2Duplicates(matrix);
}
}