1 //
2 // Copyright 2018 DreamWorks Animation LLC.
3 //
4 // Licensed under the Apache License, Version 2.0 (the "Apache License")
5 // with the following modification; you may not use this file except in
6 // compliance with the Apache License and the following modification to it:
7 // Section 6. Trademarks. is deleted and replaced with:
8 //
9 // 6. Trademarks. This License does not grant permission to use the trade
10 // names, trademarks, service marks, or product names of the Licensor
11 // and its affiliates, except as required to comply with Section 4(c) of
12 // the License and to reproduce the content of the NOTICE file.
13 //
14 // You may obtain a copy of the Apache License at
15 //
16 // http://www.apache.org/licenses/LICENSE-2.0
17 //
18 // Unless required by applicable law or agreed to in writing, software
19 // distributed under the Apache License with the above modification is
20 // distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
21 // KIND, either express or implied. See the Apache License for the specific
22 // language governing permissions and limitations under the Apache License.
23 //
24
25 #include "../far/loopPatchBuilder.h"
26 #include "../vtr/stackBuffer.h"
27 #include "../sdc/loopScheme.h"
28
29 #include <cassert>
30 #include <cmath>
31 #include <cstdio>
32
33 namespace OpenSubdiv {
34 namespace OPENSUBDIV_VERSION {
35
36 using Vtr::Array;
37 using Vtr::ConstArray;
38 using Vtr::internal::StackBuffer;
39
40 namespace Far {
41
42 //
43 // A simple struct with methods to compute Loop limit points (following
44 // the pattern established for Catmull-Clark limit points)
45 //
46 // Unlike the corresponding Catmull-Clark struct, Loop limit points are
47 // computed using the limit masks provided by the Sdc Scheme for Loop.
48 //
49 template <typename REAL>
50 struct LoopLimits {
51 public:
52 //
53 // Local classes to fulfill the interface requirements of template
54 // arguments to the Sdc::Scheme methods
55 //
56 // Class fulfilling <typename VERTEX>:
57 //
58 class LimitVertex {
59 public:
LimitVertex()60 LimitVertex() : _nFaces(0), _nEdges(0) { }
LimitVertex(int faces,int edges)61 LimitVertex(int faces, int edges) : _nFaces(faces), _nEdges(edges) { }
62
Assign(int faces,int edges)63 void Assign(int faces, int edges) { _nFaces = faces; _nEdges = edges; }
64
GetNumEdges() const65 int GetNumEdges() const { return _nEdges; }
GetNumFaces() const66 int GetNumFaces() const { return _nFaces; }
67
GetSharpnessPerEdge(float sharpness[]) const68 void GetSharpnessPerEdge(float sharpness[]) const {
69 sharpness[0] = Sdc::Crease::SHARPNESS_INFINITE;
70 for (int i = 1; i < _nEdges - 1; ++i) {
71 sharpness[i] = Sdc::Crease::SHARPNESS_SMOOTH;
72 }
73 sharpness[_nEdges - 1] = Sdc::Crease::SHARPNESS_INFINITE;
74 }
75
76 private:
77 int _nFaces;
78 int _nEdges;
79 };
80
81 //
82 // Class fulfilling <typename MASK>:
83 //
84 class LimitMask {
85 public:
86 typedef REAL Weight; // Also part of the expected interface
87
88 public:
LimitMask(Weight * w)89 LimitMask(Weight* w) : _weights(w), _valence(0) { }
~LimitMask()90 ~LimitMask() { }
91
92 public: // Generic interface expected of <typename MASK>:
GetNumVertexWeights() const93 int GetNumVertexWeights() const { return 1; }
GetNumEdgeWeights() const94 int GetNumEdgeWeights() const { return _valence; }
GetNumFaceWeights() const95 int GetNumFaceWeights() const { return 0; }
96
SetNumVertexWeights(int)97 void SetNumVertexWeights(int) { }
SetNumEdgeWeights(int count)98 void SetNumEdgeWeights(int count) { _valence = count; }
SetNumFaceWeights(int)99 void SetNumFaceWeights(int) { }
100
VertexWeight(int) const101 Weight const& VertexWeight(int) const { return _weights[0]; }
EdgeWeight(int index) const102 Weight const& EdgeWeight(int index) const { return _weights[1 + index]; }
FaceWeight(int) const103 Weight const& FaceWeight(int) const { return _weights[0]; }
104
VertexWeight(int)105 Weight& VertexWeight(int) { return _weights[0]; }
EdgeWeight(int index)106 Weight& EdgeWeight( int index) { return _weights[1 + index]; }
FaceWeight(int)107 Weight& FaceWeight( int) { return _weights[0]; }
108
AreFaceWeightsForFaceCenters() const109 bool AreFaceWeightsForFaceCenters() const { return false; }
SetFaceWeightsForFaceCenters(bool)110 void SetFaceWeightsForFaceCenters(bool) { }
111
112 private:
113 Weight* _weights;
114 int _valence;
115 };
116 public:
117 typedef Sdc::Scheme<Sdc::SCHEME_LOOP> SchemeLoop;
118
119 static void ComputeInteriorPointWeights(int valence, int faceInRing,
120 REAL* pWeights, REAL* epWeights, REAL* emWeights);
121
122 static void ComputeBoundaryPointWeights(int valence, int faceInRing,
123 REAL* pWeights, REAL* epWeights, REAL* emWeights);
124 };
125
126 template <typename REAL>
127 void
ComputeInteriorPointWeights(int valence,int faceInRing,REAL * pWeights,REAL * epWeights,REAL * emWeights)128 LoopLimits<REAL>::ComputeInteriorPointWeights(int valence, int faceInRing,
129 REAL* pWeights, REAL* epWeights, REAL* emWeights) {
130
131 int ringSize = valence + 1;
132
133 LimitVertex vertex(valence, valence);
134
135 bool computeTangentPoints = epWeights && emWeights;
136 if (!computeTangentPoints) {
137 //
138 // The interior position mask is symmetric -- no need to rotate
139 // or otherwise account for orientation:
140 //
141 LimitMask pMask(pWeights);
142
143 SchemeLoop().ComputeVertexLimitMask(vertex, pMask,
144 Sdc::Crease::RULE_SMOOTH);
145 } else {
146 //
147 // The interior tangent masks will be directed along the first
148 // two edges (the second a rotation of the first). Adjust the
149 // tangent weights for a point along the tangent, then rotate
150 // according to the face within the ring:
151 //
152 int weightWidth = valence + 1;
153 StackBuffer<REAL, 32, true> tWeights(2 * weightWidth);
154 REAL * t1Weights = &tWeights[0];
155 REAL * t2Weights = t1Weights + ringSize;
156
157 LimitMask pMask(pWeights);
158 LimitMask t1Mask(t1Weights);
159 LimitMask t2Mask(t2Weights);
160
161 SchemeLoop().ComputeVertexLimitMask(vertex, pMask, t1Mask, t2Mask,
162 Sdc::Crease::RULE_SMOOTH);
163
164 //
165 // Use the subdominant eigenvalue to scale the limit tangent t1:
166 //
167 // e = (3 + cos(2*PI/valence)) / 8
168 //
169 // Combine it with a normalizing factor of (2 / valence) to account
170 // for the scale inherent in the tangent weights, and (2 / 3) to
171 // match desired placement of the cubic point in the regular case.
172 //
173 // The weights for t1 can simply be rotated around the ring to yield
174 // t2. Combine the weights for the point in a single set for t2 and
175 // then copy it into the appropriate orientation for ep and em:
176 //
177 double theta = 2.0 * M_PI / (double) valence;
178
179 REAL tanScale = (REAL) ((3.0 + 2.0 * std::cos(theta)) / (6.0 * valence));
180
181 for (int i = 0; i < ringSize; ++i) {
182 t2Weights[i] = pWeights[i] + t1Weights[i] * tanScale;
183 }
184
185 int n1 = faceInRing;
186 int n2 = valence - n1;
187
188 epWeights[0] = t2Weights[0];
189 std::memcpy(epWeights + 1, t2Weights + 1 + n2, n1 * sizeof(REAL));
190 std::memcpy(epWeights + 1 + n1, t2Weights + 1, n2 * sizeof(REAL));
191
192 n1 = (faceInRing + 1) % valence;
193 n2 = valence - n1;
194
195 emWeights[0] = t2Weights[0];
196 std::memcpy(emWeights + 1, t2Weights + 1 + n2, n1 * sizeof(REAL));
197 std::memcpy(emWeights + 1 + n1, t2Weights + 1, n2 * sizeof(REAL));
198 }
199 }
200
201 template <typename REAL>
202 void
ComputeBoundaryPointWeights(int valence,int faceInRing,REAL * pWeights,REAL * epWeights,REAL * emWeights)203 LoopLimits<REAL>::ComputeBoundaryPointWeights(int valence, int faceInRing,
204 REAL* pWeights, REAL* epWeights, REAL* emWeights) {
205
206 LimitVertex vertex(valence - 1, valence);
207
208 bool computeTangentPoints = epWeights && emWeights;
209 if (!computeTangentPoints) {
210 //
211 // The boundary position mask will be assigned non-zero weights
212 // for the vertex and its first and last edges:
213 //
214 LimitMask pMask(pWeights);
215
216 SchemeLoop().ComputeVertexLimitMask(vertex, pMask,
217 Sdc::Crease::RULE_CREASE);
218 } else {
219 //
220 // The boundary tangent masks need more explicit handling than
221 // the interior. One of the tangents will be along the boundary
222 // and the other towards the interior, but one or both edge
223 // points may be along interior edges. A boundary edge point is
224 // easy to deal with once identified, but interior edge points
225 // need a numerical rotation of the interior tangent to orient it
226 // along the desired edge.
227 //
228 int weightWidth = valence + 1;
229 StackBuffer<REAL, 32, true> tWeights(2 * weightWidth);
230 REAL * t1Weights = &tWeights[0];
231 REAL * t2Weights = t1Weights + weightWidth;
232
233 REAL t1Leading = (REAL) (1.0 / 6.0);
234 REAL t1Trailing = -(REAL) (1.0 / 6.0);
235
236 REAL t2Scale = (REAL) (1.0 / 24.0);
237
238 LimitMask pMask(pWeights);
239 LimitMask t1Mask(t1Weights);
240 LimitMask t2Mask(t2Weights);
241
242 SchemeLoop().ComputeVertexLimitMask(vertex, pMask, t1Mask, t2Mask,
243 Sdc::Crease::RULE_CREASE);
244
245 bool epOnLeadingEdge = (faceInRing == 0);
246 bool emOnTrailingEdge = (faceInRing == (valence - 1));
247
248 if (epOnLeadingEdge) {
249 std::memset(&epWeights[0], 0, weightWidth * sizeof(REAL));
250
251 epWeights[0] = (REAL) (2.0 / 3.0);
252 epWeights[1] = (REAL) (1.0 / 3.0);
253 } else {
254 int iEdgeNext = faceInRing;
255 REAL faceAngle = (REAL) (M_PI / (double)(valence - 1));
256 REAL faceAngleNext = faceAngle * (REAL)iEdgeNext;
257 REAL cosAngleNext = std::cos(faceAngleNext);
258 REAL sinAngleNext = std::sin(faceAngleNext);
259
260 for (int i = 0; i < weightWidth; ++i) {
261 epWeights[i] = t2Scale * t2Weights[i] * sinAngleNext;
262 }
263 epWeights[0] += pWeights[0];
264 epWeights[1] += pWeights[1] + t1Leading * cosAngleNext;
265 epWeights[valence] += pWeights[valence] + t1Trailing * cosAngleNext;
266 }
267
268 if (emOnTrailingEdge) {
269 std::memset(&emWeights[0], 0, weightWidth * sizeof(REAL));
270
271 emWeights[0] = (REAL) (2.0 / 3.0);
272 emWeights[valence] = (REAL) (1.0 / 3.0);
273 } else {
274 int iEdgePrev = (faceInRing + 1) % valence;
275 REAL faceAngle = (REAL) (M_PI / (double)(valence - 1));
276 REAL faceAnglePrev = faceAngle * (REAL)iEdgePrev;
277 REAL cosAnglePrev = std::cos(faceAnglePrev);
278 REAL sinAnglePrev = std::sin(faceAnglePrev);
279
280 for (int i = 0; i < weightWidth; ++i) {
281 emWeights[i] = t2Scale * t2Weights[i] * sinAnglePrev;
282 }
283 emWeights[0] += pWeights[0];
284 emWeights[1] += pWeights[1] + t1Leading * cosAnglePrev;
285 emWeights[valence] += pWeights[valence] + t1Trailing * cosAnglePrev;
286 }
287 }
288 }
289
290
291 //
292 // SparseMatrixRow
293 //
294 // This was copied from the CatmarkPatchBuilder as a starting point, so
295 // comments below relate to the state of CatmarkPatchBuilder...
296 //
297 // This is a utility class representing a row of a SparseMatrix -- which
298 // in turn corresponds to a point of a resulting patch. Instances of this
299 // class are intended to encapsulate the contributions of a point and be
300 // passed to functions as such.
301 //
302 // (Consider moving this to PatchBuilder as a protected class or maybe a
303 // public class within SparseMatrix itself, e.g. SparseMatrix<REAL>::Row.)
304 //
305 namespace {
306 template <typename REAL>
307 class SparseMatrixRow {
308 public:
SparseMatrixRow(SparseMatrix<REAL> & matrix,int row)309 SparseMatrixRow(SparseMatrix<REAL> & matrix, int row) :
310 _size(matrix.GetRowSize(row)),
311 _indices(matrix.SetRowColumns(row).begin()),
312 _weights(matrix.SetRowElements(row).begin()) { }
313
GetSize() const314 int GetSize() const { return _size; }
315
SetWeight(int i,REAL weight)316 void SetWeight(int i, REAL weight) { _weights[i] = weight; }
317
Assign(int rowEntry,Index index,REAL weight)318 void Assign(int rowEntry, Index index, REAL weight) {
319 _indices[rowEntry] = index;
320 _weights[rowEntry] = weight;
321 }
322
Copy(SparseMatrixRow<REAL> const & other)323 void Copy(SparseMatrixRow<REAL> const & other) {
324 assert(GetSize() == other.GetSize());
325 std::memcpy(_indices, other._indices, _size * sizeof(Index));
326 std::memcpy(_weights, other._weights, _size * sizeof(REAL));
327 }
328
329 public:
330 int _size;
331 Index * _indices;
332 REAL * _weights;
333 };
334 } // end namespace
335
336 //
337 // Simple utility functions for dealing with SparseMatrix:
338 //
339 namespace {
340 #ifdef __INTEL_COMPILER
341 #pragma warning (push)
342 #pragma warning disable 1572
343 #endif
344 template <typename REAL>
_isZero(REAL w)345 inline bool _isZero(REAL w) { return (w == (REAL)0.0); }
346 #ifdef __INTEL_COMPILER
347 #pragma warning (pop)
348 #endif
349
350 template <typename REAL>
351 void
_initializeFullMatrix(SparseMatrix<REAL> & M,int nRows,int nColumns)352 _initializeFullMatrix(SparseMatrix<REAL> & M, int nRows, int nColumns) {
353
354 M.Resize(nRows, nColumns, nRows * nColumns);
355
356 // Fill row 0 with index for every column:
357 M.SetRowSize(0, nColumns);
358 Array<int> row0Columns = M.SetRowColumns(0);
359 for (int i = 0; i < nColumns; ++i) {
360 row0Columns[i] = i;
361 }
362
363 // Copy row 0's indices into all other rows:
364 for (int row = 1; row < nRows; ++row) {
365 M.SetRowSize(row, nColumns);
366 Array<int> dstRowColumns = M.SetRowColumns(row);
367 std::memcpy(&dstRowColumns[0], &row0Columns[0], nColumns * sizeof(int));
368 }
369 }
370
371 template <typename REAL>
372 void
_resizeMatrix(SparseMatrix<REAL> & matrix,int numRows,int numColumns,int numElements,int const rowSizes[])373 _resizeMatrix(SparseMatrix<REAL> & matrix,
374 int numRows, int numColumns, int numElements,
375 int const rowSizes[]) {
376
377 matrix.Resize(numRows, numColumns, numElements);
378 for (int i = 0; i < numRows; ++i) {
379 matrix.SetRowSize(i, rowSizes[i]);
380 }
381 assert(matrix.GetNumElements() == numElements);
382 }
383
384 template <typename REAL>
385 void
_addSparsePointToFullRow(REAL * fullRow,SparseMatrixRow<REAL> const & p,REAL s,int * indexMask)386 _addSparsePointToFullRow(REAL * fullRow,
387 SparseMatrixRow<REAL> const & p,
388 REAL s, int * indexMask) {
389
390 for (int i = 0; i < p.GetSize(); ++i) {
391 int index = p._indices[i];
392
393 fullRow[index] += s * p._weights[i];
394
395 indexMask[index] = 1 + index;
396 }
397 }
398
399 template <typename REAL>
400 void
_combineSparsePointsInFullRow(SparseMatrixRow<REAL> & p,REAL aCoeff,SparseMatrixRow<REAL> const & a,REAL bCoeff,SparseMatrixRow<REAL> const & b,int rowSize,REAL * rowBuffer,int * maskBuffer)401 _combineSparsePointsInFullRow(SparseMatrixRow<REAL> & p,
402 REAL aCoeff, SparseMatrixRow<REAL> const & a,
403 REAL bCoeff, SparseMatrixRow<REAL> const & b,
404 int rowSize, REAL * rowBuffer, int * maskBuffer) {
405
406 std::memset(maskBuffer, 0, rowSize * sizeof(int));
407 std::memset(rowBuffer, 0, rowSize * sizeof(REAL));
408
409 _addSparsePointToFullRow(rowBuffer, a, aCoeff, maskBuffer);
410 _addSparsePointToFullRow(rowBuffer, b, bCoeff, maskBuffer);
411
412 int nWeights = 0;
413 for (int i = 0; i < rowSize; ++i) {
414 if (maskBuffer[i]) {
415 p.Assign(nWeights++, maskBuffer[i] - 1, rowBuffer[i]);
416 }
417 }
418 assert(nWeights <= p.GetSize());
419 for (int i = nWeights; i < p.GetSize(); ++i, ++nWeights) {
420 p.Assign(i, 0, 0.0f);
421 }
422 }
423
424 template <typename REAL>
425 void
_addSparseRowToFull(REAL * fullRow,SparseMatrix<REAL> const & M,int sparseRow,REAL s)426 _addSparseRowToFull(REAL * fullRow,
427 SparseMatrix<REAL> const & M, int sparseRow, REAL s) {
428
429 ConstArray<int> indices = M.GetRowColumns(sparseRow);
430 ConstArray<REAL> weights = M.GetRowElements(sparseRow);
431
432 for (int i = 0; i < indices.size(); ++i) {
433 fullRow[indices[i]] += s * weights[i];
434 }
435 }
436
437 template <typename REAL>
438 void
_combineSparseMatrixRowsInFull(SparseMatrix<REAL> & dstMatrix,int dstRowIndex,SparseMatrix<REAL> const & srcMatrix,int numSrcRows,int const srcRowIndices[],REAL const * srcRowWeights)439 _combineSparseMatrixRowsInFull(SparseMatrix<REAL> & dstMatrix, int dstRowIndex,
440 SparseMatrix<REAL> const & srcMatrix, int numSrcRows,
441 int const srcRowIndices[], REAL const * srcRowWeights) {
442
443 REAL * dstRow = &dstMatrix.SetRowElements(dstRowIndex)[0];
444
445 std::memset(dstRow, 0, dstMatrix.GetNumColumns() * sizeof(REAL));
446
447 for (int i = 0; i < numSrcRows; ++i) {
448 if (!_isZero(srcRowWeights[i])) {
449 _addSparseRowToFull(dstRow, srcMatrix, srcRowIndices[i], srcRowWeights[i]);
450 }
451 }
452 }
453
454 template <typename REAL>
455 void
_matrixPrintDensity(const char * prefix,SparseMatrix<REAL> const & M)456 _matrixPrintDensity(const char* prefix, SparseMatrix<REAL> const & M) {
457
458 int fullSize = M.GetNumRows() * M.GetNumColumns();
459 int sparseSize = M.GetNumElements();
460
461 int nonZeroSize = 0;
462 for (int i = 0; i < M.GetNumRows(); ++i) {
463 ConstArray<REAL> elements = M.GetRowElements(i);
464 for (int j = 0; j < elements.size(); ++j) {
465 nonZeroSize += (elements[j] != 0);
466 }
467 }
468 printf("%s(%dx%d = %d): elements = %d, non-zero = %d, density = %.1f\n",
469 prefix, M.GetNumRows(), M.GetNumColumns(), fullSize,
470 sparseSize, nonZeroSize, (REAL)nonZeroSize * 100.0f / (REAL)fullSize);
471 }
472
473 //
474 // The valence-2 interior case poses problems for the way patch points
475 // are computed as combinations of source points and stored as a row in
476 // a SparseMatrix. An interior vertex of valence-2 causes duplicate
477 // vertices to appear in the 1-rings of its neighboring vertices and we
478 // want the entries of a SparseMatrix row to be unique.
479 //
480 // For the most part, this does not pose a problem while the matrix (set
481 // of patch points) is being constructed, so we leave those duplicate
482 // entries in place and deal with them as a post-process here.
483 //
484 // The SourcePatch is also sensitive to the presence of such valence-2
485 // vertices for its own reasons (it needs to identifiy a unique set of
486 // source points from a set of corner rings), so a simple query of its
487 // corners indicates when this post-process is necessary. (And since
488 // this case is a rare occurrence, efficiency is not a major concern.)
489 //
490 template <typename REAL>
_removeValence2Duplicates(SparseMatrix<REAL> & M)491 void _removeValence2Duplicates(SparseMatrix<REAL> & M) {
492
493 // This will later be determined by the PatchBuilder member:
494 int regFaceSize = 4;
495
496 SparseMatrix<REAL> T;
497 T.Resize(M.GetNumRows(), M.GetNumColumns(), M.GetNumElements());
498
499 int nRows = M.GetNumRows();
500 for (int row = 0; row < nRows; ++row) {
501 int srcRowSize = M.GetRowSize(row);
502
503 int const * srcIndices = M.GetRowColumns(row).begin();
504 REAL const * srcWeights = M.GetRowElements(row).begin();
505
506 // Scan the entries to see if there are duplicates -- copy
507 // the row if not, otherwise, need to compress it:
508 bool cornerUsed[4] = { false, false, false, false };
509
510 int srcDupCount = 0;
511 for (int i = 0; i < srcRowSize; ++i) {
512 int srcIndex = srcIndices[i];
513 if (srcIndex < regFaceSize) {
514 srcDupCount += (int) cornerUsed[srcIndex];
515 cornerUsed[srcIndex] = true;
516 }
517 }
518
519 // Size this row for the destination and copy or compress:
520 T.SetRowSize(row, srcRowSize - srcDupCount);
521
522 int* dstIndices = T.SetRowColumns(row).begin();
523 REAL* dstWeights = T.SetRowElements(row).begin();
524
525 if (srcDupCount) {
526 REAL * cornerDstPtr[4] = { 0, 0, 0, 0 };
527
528 for (int i = 0; i < srcRowSize; ++i) {
529 int srcIndex = *srcIndices++;
530 REAL srcWeight = *srcWeights++;
531
532 if (srcIndex < regFaceSize) {
533 if (cornerDstPtr[srcIndex]) {
534 *cornerDstPtr[srcIndex] += srcWeight;
535 continue;
536 } else {
537 cornerDstPtr[srcIndex] = dstWeights;
538 }
539 }
540 *dstIndices++ = srcIndex;
541 *dstWeights++ = srcWeight;
542 }
543 } else {
544 std::memcpy(&dstIndices[0], &srcIndices[0], srcRowSize * sizeof(Index));
545 std::memcpy(&dstWeights[0], &srcWeights[0], srcRowSize * sizeof(REAL));
546 }
547 }
548 M.Swap(T);
549 }
550 } // end namespace for SparseMatrix utilities
551
552
553 //
554 // GregoryTriConverter
555 //
556 // The GregoryTriConverter class provides a change-of-basis matrix from source
557 // vertices in a Loop mesh to the 18 control points of a quartic Gregory triangle.
558 //
559 // The quartic triangle is first constructed as a cubic/quartic hybrid -- with
560 // cubic boundary curves and cross-boundary continuity formulated in terms of
561 // cubics. The result is then raised to a full quartic once continuity across
562 // all boundaries is achieved. In most cases 2 of the 3 boundaries will be
563 // cubic (though now represented as quartic) and only one boundary need be a
564 // true quartic to meet a regular Box-spline patch.
565 //
566 // Control points are labeled using the convention adopted for quads, with
567 // Ep and Em referring to the "plus" and "minus" edge points and similarly
568 // for the face points Fp and Fm. The additional quartic "mid-edge" points
569 // associated with each boundary are referred to as M.
570 //
571 template <typename REAL>
572 class GregoryTriConverter {
573 public:
574 typedef REAL Weight;
575 typedef SparseMatrix<Weight> Matrix;
576 typedef SparseMatrixRow<Weight> Point;
577 public:
GregoryTriConverter()578 GregoryTriConverter() : _numSourcePoints(0) { }
579 GregoryTriConverter(SourcePatch const & sourcePatch);
580 GregoryTriConverter(SourcePatch const & sourcePatch, Matrix & sparseMatrix);
581
582 void Initialize(SourcePatch const & sourcePatch);
583
IsIsolatedInteriorPatch() const584 bool IsIsolatedInteriorPatch() const { return _isIsolatedInteriorPatch; }
HasVal2InteriorCorner() const585 bool HasVal2InteriorCorner() const { return _hasVal2InteriorCorner; }
GetIsolatedInteriorCorner() const586 int GetIsolatedInteriorCorner() const { return _isolatedCorner; }
GetIsolatedInteriorValence() const587 int GetIsolatedInteriorValence() const { return _isolatedValence; }
588
589 void Convert(Matrix & sparseMatrix) const;
590
591 private:
592 //
593 // Local nested class for GregoryTriConverter to cache information for
594 // the corners of the source patch. It copies some information from the
595 // SourcePatch so that we don't have to keep it around, but it contains
596 // additional information relevant to the determination of the Gregory
597 // points -- most notably classifications of the face-points and the
598 // cosines of angles for the face corners that are used repeatedly.
599 //
600 struct CornerTopology {
601 // Basic flags copied from the SourcePatch
602 unsigned int isBoundary : 1;
603 unsigned int isSharp : 1;
604 unsigned int isDart : 1;
605 unsigned int isRegular : 1;
606 unsigned int isVal2Int : 1;
607 unsigned int isCorner : 1;
608
609 // Flags for edge- and face-points relating to adjacent corners:
610 unsigned int epOnBoundary : 1;
611 unsigned int emOnBoundary : 1;
612
613 unsigned int fpIsRegular : 1;
614 unsigned int fmIsRegular : 1;
615 unsigned int fpIsCopied : 1;
616 unsigned int fmIsCopied : 1;
617
618 // Other values stored for repeated use:
619 int valence;
620 int numFaces;
621 int faceInRing;
622
623 REAL faceAngle;
624 REAL cosFaceAngle;
625
626 // Its useful to have the ring for each corner immediately available:
627 StackBuffer<int, 30, true> ringPoints;
628 };
629
630 //
631 // Methods to resize the matrix before populating it:
632 //
633 void resizeMatrixUnisolated(Matrix & matrix) const;
634
635 void resizeMatrixIsolatedIrregular(Matrix & matrix, int irregCornerIndex,
636 int irregValence) const;
637
638 //
639 // Methods to compute the various rows of points in the matrix:
640 //
641 void assignRegularEdgePoints(int cornerIndex, Matrix & matrix) const;
642 void computeIrregularEdgePoints(int cornerIndex, Matrix & matrix,
643 Weight *weightBuffer) const;
644
645 void assignRegularFacePoints(int cornerIndex, Matrix & matrix) const;
646 void computeIrregularFacePoints(int cornerIndex, Matrix & matrix,
647 Weight *weightBuffer, int *indexBuffer) const;
648
649 void computeIrregularInteriorEdgePoints(int cornerIndex,
650 Point & P, Point & Ep, Point & Em,
651 Weight *weightBuffer) const;
652 void computeIrregularBoundaryEdgePoints(int cornerIndex,
653 Point & P, Point & Ep, Point & Em,
654 Weight *weightBuffer) const;
655
656 int getIrregularFacePointSize(int cornerNear, int cornerFar) const;
657 void computeIrregularFacePoint(
658 int cornerNear, int edgeInNearRing, int cornerFar,
659 Point const & p, Point const & eNear, Point const & eFar,
660 Point & fNear, REAL signForSideOfEdge /* -1.0 or 1.0 */,
661 Weight *rowWeights, int *columnMask) const;
662
663 void assignRegularMidEdgePoint(int edgeIndex, Matrix & matrix) const;
664 void computeIrregularMidEdgePoint(int edgeIndex, Matrix & matrix,
665 Weight *rowWeights, int *columnMask) const;
666
667 void promoteCubicEdgePointsToQuartic(Matrix & matrix,
668 Weight *rowWeights, int *columnMask) const;
669
670 private:
671 int _numSourcePoints;
672 int _maxValence;
673
674 bool _isIsolatedInteriorPatch;
675 bool _hasVal2InteriorCorner;
676 int _isolatedCorner;
677 int _isolatedValence;
678
679 CornerTopology _corners[3];
680 };
681
682
683 //
684 // GregoryTriConverter
685 //
686 // Construction and initialization/computation of the change-of-basis
687 // matrix to a Gregory patch.
688 //
689 template <typename REAL>
GregoryTriConverter(SourcePatch const & sourcePatch)690 GregoryTriConverter<REAL>::GregoryTriConverter(
691 SourcePatch const & sourcePatch) {
692
693 Initialize(sourcePatch);
694 }
695
696 template <typename REAL>
GregoryTriConverter(SourcePatch const & sourcePatch,Matrix & sparseMatrix)697 GregoryTriConverter<REAL>::GregoryTriConverter(
698 SourcePatch const & sourcePatch, Matrix & sparseMatrix) {
699
700 Initialize(sourcePatch);
701 Convert(sparseMatrix);
702 }
703
704 template <typename REAL>
705 void
Initialize(SourcePatch const & sourcePatch)706 GregoryTriConverter<REAL>::Initialize(SourcePatch const & sourcePatch) {
707
708 //
709 // Allocate and gather the 1-rings for the corner vertices and other
710 // topological information for more immediate access:
711 //
712 int width = sourcePatch.GetNumSourcePoints();
713 _numSourcePoints = width;
714 _maxValence = sourcePatch.GetMaxValence();
715
716 int boundaryCount = 0;
717 int irregularCount = 0;
718 int irregularCorner = -1;
719 int irregularValence = -1;
720 int sharpCount = 0;
721 int val2IntCount = 0;
722
723 for (int cIndex = 0; cIndex < 3; ++cIndex) {
724 SourcePatch::Corner srcCorner = sourcePatch._corners[cIndex];
725
726 CornerTopology& corner = _corners[cIndex];
727
728 corner.isBoundary = srcCorner._boundary;
729 corner.isSharp = srcCorner._sharp;
730 corner.isDart = srcCorner._dart;
731 corner.isCorner = (srcCorner._numFaces == 1);
732 corner.numFaces = srcCorner._numFaces;
733 corner.faceInRing = srcCorner._patchFace;
734 corner.isVal2Int = srcCorner._val2Interior;
735 corner.valence = corner.numFaces + corner.isBoundary;
736
737 corner.isRegular = ((corner.numFaces << corner.isBoundary) == 6)
738 && !corner.isSharp;
739 if (corner.isRegular) {
740 corner.faceAngle = (REAL)(M_PI / 3.0);
741 corner.cosFaceAngle = 0.5f;
742 } else {
743 corner.faceAngle =
744 (corner.isBoundary ? REAL(M_PI) : REAL(2.0 * M_PI))
745 / REAL(corner.numFaces);
746 corner.cosFaceAngle = std::cos(corner.faceAngle);
747 }
748
749 corner.ringPoints.SetSize(sourcePatch.GetCornerRingSize(cIndex));
750 sourcePatch.GetCornerRingPoints(cIndex, corner.ringPoints);
751
752 // Accumulate topology information to categorize the patch as a whole:
753 boundaryCount += corner.isBoundary;
754 if (!corner.isRegular) {
755 irregularCount ++;
756 irregularCorner = cIndex;
757 irregularValence = corner.valence;
758 }
759 sharpCount += corner.isSharp;
760 val2IntCount += corner.isVal2Int;
761 }
762
763 // Make a second pass to assign tags dependent on adjacent corners
764 for (int cIndex = 0; cIndex < 3; ++cIndex) {
765 CornerTopology& corner = _corners[cIndex];
766
767 int cNext = (cIndex + 1) % 3;
768 int cPrev = (cIndex + 2) % 3;
769
770 corner.epOnBoundary = false;
771 corner.emOnBoundary = false;
772
773 //
774 // Identify if the face points are regular or shared/copied from
775 // one of the pair:
776 //
777 corner.fpIsRegular = corner.isRegular && _corners[cNext].isRegular;
778 corner.fmIsRegular = corner.isRegular && _corners[cPrev].isRegular;
779
780 corner.fpIsCopied = false;
781 corner.fmIsCopied = false;
782
783 if (corner.isBoundary) {
784 corner.epOnBoundary = (corner.faceInRing == 0);
785 corner.emOnBoundary = (corner.faceInRing == (corner.numFaces - 1));
786
787 // Both face points are same when one of the two corners' edges
788 // is discontinuous -- one is then copied from the other (unless
789 // regular)
790 if (corner.numFaces > 1) {
791 if (corner.epOnBoundary) {
792 corner.fpIsRegular = corner.fmIsRegular;
793 corner.fpIsCopied = !corner.fpIsRegular;
794 }
795 if (corner.emOnBoundary) {
796 corner.fmIsRegular = corner.fpIsRegular;
797 corner.fmIsCopied = !corner.fmIsRegular;
798 }
799 } else {
800 // The case of a corner patch is always regular
801 corner.fpIsRegular = true;
802 corner.fmIsRegular = true;
803 }
804 }
805 }
806 _isIsolatedInteriorPatch = (irregularCount == 1) && (boundaryCount == 0) &&
807 (irregularValence > 2) && (sharpCount == 0);
808 if (_isIsolatedInteriorPatch) {
809 _isolatedCorner = irregularCorner;
810 _isolatedValence = irregularValence;
811 }
812 _hasVal2InteriorCorner = (val2IntCount > 0);
813 }
814
815 template <typename REAL>
816 void
Convert(Matrix & matrix) const817 GregoryTriConverter<REAL>::Convert(Matrix & matrix) const {
818
819 //
820 // Initialize the sparse matrix to accomodate the coefficients for each
821 // row/point -- identify common topological cases to treat more easily
822 // (and note that specializing the population of the matrix may also be
823 // worthwhile in such cases)
824 //
825 if (_isIsolatedInteriorPatch) {
826 resizeMatrixIsolatedIrregular(matrix, _isolatedCorner, _isolatedValence);
827 } else {
828 resizeMatrixUnisolated(matrix);
829 }
830
831 //
832 // Compute the corner and edge points P, Ep and Em first. Since face
833 // points Fp and Fm involve edge points for two adjacent corners, their
834 // computation must follow:
835 //
836 int maxRingSize = 1 + _maxValence;
837 int weightBufferSize = std::max(3 * maxRingSize, 2 * _numSourcePoints);
838
839 StackBuffer<Weight, 128, true> weightBuffer(weightBufferSize);
840 StackBuffer<int, 128, true> indexBuffer(weightBufferSize);
841
842 for (int cIndex = 0; cIndex < 3; ++cIndex) {
843 if (_corners[cIndex].isRegular) {
844 assignRegularEdgePoints(cIndex, matrix);
845 } else {
846 computeIrregularEdgePoints(cIndex, matrix, weightBuffer);
847 }
848 }
849
850 for (int cIndex = 0; cIndex < 3; ++cIndex) {
851 if (_corners[cIndex].fpIsRegular || _corners[cIndex].fmIsRegular) {
852 assignRegularFacePoints(cIndex, matrix);
853 }
854 if (!_corners[cIndex].fpIsRegular || !_corners[cIndex].fmIsRegular) {
855 computeIrregularFacePoints(cIndex, matrix, weightBuffer, indexBuffer);
856 }
857 }
858
859 for (int eIndex = 0; eIndex < 3; ++eIndex) {
860 CornerTopology const & c0 = _corners[eIndex];
861 CornerTopology const & c1 = _corners[(eIndex + 1) % 3];
862
863 bool isBoundaryEdge = c0.epOnBoundary && c1.emOnBoundary;
864 bool isDartEdge = c0.epOnBoundary != c1.emOnBoundary;
865 if (isBoundaryEdge || (c0.isRegular && c1.isRegular && !isDartEdge)) {
866 assignRegularMidEdgePoint(eIndex, matrix);
867 } else {
868 computeIrregularMidEdgePoint(eIndex, matrix, weightBuffer, indexBuffer);
869 }
870 }
871 promoteCubicEdgePointsToQuartic(matrix, weightBuffer, indexBuffer);
872
873 if (_hasVal2InteriorCorner) {
874 _removeValence2Duplicates(matrix);
875 }
876 }
877
878 template <typename REAL>
879 void
resizeMatrixIsolatedIrregular(Matrix & matrix,int cornerIndex,int cornerValence) const880 GregoryTriConverter<REAL>::resizeMatrixIsolatedIrregular(
881 Matrix & matrix, int cornerIndex, int cornerValence) const {
882
883 int irregRingSize = 1 + cornerValence;
884
885 int irregCorner = cornerIndex;
886 int irregPlus = (cornerIndex + 1) % 3;
887 int irregMinus = (cornerIndex + 2) % 3;
888
889 int rowSizes[18];
890 int * rowSizePtr = 0;
891
892 rowSizePtr = rowSizes + irregCorner * 5;
893 *rowSizePtr++ = irregRingSize;
894 *rowSizePtr++ = irregRingSize;
895 *rowSizePtr++ = irregRingSize;
896 *rowSizePtr++ = 3 + irregRingSize;
897 *rowSizePtr++ = 3 + irregRingSize;
898
899 rowSizePtr = rowSizes + irregPlus * 5;
900 *rowSizePtr++ = 7;
901 *rowSizePtr++ = 7;
902 *rowSizePtr++ = 7;
903 *rowSizePtr++ = 5;
904 *rowSizePtr++ = 3 + irregRingSize;
905
906 rowSizePtr = rowSizes + irregMinus * 5;
907 *rowSizePtr++ = 7;
908 *rowSizePtr++ = 7;
909 *rowSizePtr++ = 7;
910 *rowSizePtr++ = 3 + irregRingSize;
911 *rowSizePtr++ = 5;
912
913 // The 3 quartic mid-edge points are not grouped with corners:
914 rowSizes[15 + irregCorner] = 3 + irregRingSize;
915 rowSizes[15 + irregPlus] = 4;
916 rowSizes[15 + irregMinus] = 3 + irregRingSize;
917
918 int numElements = 9*irregRingSize + 74;
919
920 _resizeMatrix(matrix, 18, _numSourcePoints, numElements, rowSizes);
921 }
922
923 template <typename REAL>
924 void
resizeMatrixUnisolated(Matrix & matrix) const925 GregoryTriConverter<REAL>::resizeMatrixUnisolated(Matrix & matrix) const {
926
927 int rowSizes[18];
928
929 int numElements = 0;
930
931 for (int cIndex = 0; cIndex < 3; ++cIndex) {
932 int * rowSize = rowSizes + cIndex*5;
933
934 CornerTopology const & corner = _corners[cIndex];
935
936 // First, the corner and pair of edge points:
937 if (corner.isRegular) {
938 if (! corner.isBoundary) {
939 rowSize[0] = 7;
940 rowSize[1] = 7;
941 rowSize[2] = 7;
942 } else {
943 rowSize[0] = 3;
944 rowSize[1] = corner.epOnBoundary ? 3 : 5;
945 rowSize[2] = corner.emOnBoundary ? 3 : 5;
946 }
947 } else {
948 if (corner.isSharp) {
949 rowSize[0] = 1;
950 rowSize[1] = 2;
951 rowSize[2] = 2;
952 } else if (! corner.isBoundary) {
953 int ringSize = 1 + corner.valence;
954 rowSize[0] = ringSize;
955 rowSize[1] = ringSize;
956 rowSize[2] = ringSize;
957 } else if (corner.numFaces > 1) {
958 int ringSize = 1 + corner.valence;
959 rowSize[0] = 3;
960 rowSize[1] = corner.epOnBoundary ? 3 : ringSize;
961 rowSize[2] = corner.emOnBoundary ? 3 : ringSize;
962 } else {
963 rowSize[0] = 3;
964 rowSize[1] = 3;
965 rowSize[2] = 3;
966 }
967 }
968 numElements += rowSize[0] + rowSize[1] + rowSize[2];
969
970 // Second, the pair of face points:
971 rowSize[3] = 5 - corner.epOnBoundary - corner.emOnBoundary;
972 rowSize[4] = 5 - corner.epOnBoundary - corner.emOnBoundary;
973 if (!corner.fpIsRegular || !corner.fmIsRegular) {
974 int cNext = (cIndex + 1) % 3;
975 int cPrev = (cIndex + 2) % 3;
976 if (!corner.fpIsRegular) {
977 rowSize[3] = getIrregularFacePointSize(cIndex,
978 corner.fpIsCopied ? cPrev : cNext);
979 }
980 if (!corner.fmIsRegular) {
981 rowSize[4] = getIrregularFacePointSize(cIndex,
982 corner.fmIsCopied ? cNext : cPrev);
983 }
984 }
985 numElements += rowSize[3] + rowSize[4];
986
987 // Third, the quartic mid-edge boundary point (edge following corner):
988 int cNext = (cIndex + 1) % 3;
989 CornerTopology const & cornerNext = _corners[cNext];
990
991 int & midEdgeSize = rowSizes[15 + cIndex];
992
993 if (corner.epOnBoundary && cornerNext.emOnBoundary) {
994 midEdgeSize = 2;
995 } else if (corner.isRegular && cornerNext.isRegular &&
996 (corner.epOnBoundary == cornerNext.emOnBoundary)) {
997 midEdgeSize = 4;
998 } else {
999 midEdgeSize = getIrregularFacePointSize(cIndex, cNext);
1000 }
1001 numElements += midEdgeSize;
1002 }
1003 _resizeMatrix(matrix, 18, _numSourcePoints, numElements, rowSizes);
1004 }
1005
1006 template <typename REAL>
1007 void
assignRegularEdgePoints(int cIndex,Matrix & matrix) const1008 GregoryTriConverter<REAL>::assignRegularEdgePoints(int cIndex, Matrix & matrix) const {
1009
1010 Point p (matrix, 5*cIndex + 0);
1011 Point ep(matrix, 5*cIndex + 1);
1012 Point em(matrix, 5*cIndex + 2);
1013
1014 CornerTopology const & corner = _corners[cIndex];
1015
1016 int const * cRing = corner.ringPoints;
1017
1018 if (! corner.isBoundary) {
1019 REAL const pScale = (REAL) (1.0 / 12.0);
1020
1021 p.Assign(0, cIndex, 0.5f);
1022 p.Assign(1, cRing[0], pScale);
1023 p.Assign(2, cRing[1], pScale);
1024 p.Assign(3, cRing[2], pScale);
1025 p.Assign(4, cRing[3], pScale);
1026 p.Assign(5, cRing[4], pScale);
1027 p.Assign(6, cRing[5], pScale);
1028 assert(p.GetSize() == 7);
1029
1030 // Identify the edges along Ep and Em and those opposite them:
1031 REAL const eWeights[6] = { 7.0f, 5.0f, 1.0f, -1.0f, 1.0f, 5.0f };
1032 REAL const eScale = (REAL) (1.0 / 36.0);
1033
1034 int iEdgeEp = corner.faceInRing;
1035 int iEdgeEm = (corner.faceInRing + 1) % 6;
1036
1037 ep.Assign(0, cIndex, 0.5f);
1038 ep.Assign(1, cRing[iEdgeEp], eScale * eWeights[0]);
1039 ep.Assign(2, cRing[(iEdgeEp + 1) % 6], eScale * eWeights[1]);
1040 ep.Assign(3, cRing[(iEdgeEp + 2) % 6], eScale * eWeights[2]);
1041 ep.Assign(4, cRing[(iEdgeEp + 3) % 6], eScale * eWeights[3]);
1042 ep.Assign(5, cRing[(iEdgeEp + 4) % 6], eScale * eWeights[4]);
1043 ep.Assign(6, cRing[(iEdgeEp + 5) % 6], eScale * eWeights[5]);
1044 assert(ep.GetSize() == 7);
1045
1046 em.Assign(0, cIndex, 0.5f);
1047 em.Assign(1, cRing[iEdgeEm], eScale * eWeights[0]);
1048 em.Assign(2, cRing[(iEdgeEm + 1) % 6], eScale * eWeights[1]);
1049 em.Assign(3, cRing[(iEdgeEm + 2) % 6], eScale * eWeights[2]);
1050 em.Assign(4, cRing[(iEdgeEm + 3) % 6], eScale * eWeights[3]);
1051 em.Assign(5, cRing[(iEdgeEm + 4) % 6], eScale * eWeights[4]);
1052 em.Assign(6, cRing[(iEdgeEm + 5) % 6], eScale * eWeights[5]);
1053 assert(em.GetSize() == 7);
1054 } else {
1055 REAL oneThird = (REAL) (1.0 / 3.0);
1056 REAL twoThirds = (REAL) (2.0 / 3.0);
1057 REAL oneSixth = (REAL) (1.0 / 6.0);
1058
1059 p.Assign(0, cIndex, twoThirds);
1060 p.Assign(1, cRing[0], oneSixth);
1061 p.Assign(2, cRing[3], oneSixth);
1062 assert(p.GetSize() == 3);
1063
1064 //
1065 // We have three triangles here, and the two edge points may be along two
1066 // of four edges -- two of which are interior and require weights adjusted
1067 // from above to account for phantom points (yielding {1/2, 1/6, 1/6, 1/6})
1068 //
1069 if (corner.epOnBoundary) {
1070 ep.Assign(0, cIndex, twoThirds);
1071 ep.Assign(1, cRing[0], oneThird);
1072 ep.Assign(2, cRing[3], 0.0f);
1073 assert(ep.GetSize() == 3);
1074 } else {
1075 ep.Assign(0, cIndex, 0.5f);
1076 ep.Assign(1, cRing[1], oneSixth);
1077 ep.Assign(2, cRing[2], oneSixth);
1078 ep.Assign(3, cRing[corner.emOnBoundary ? 3 : 0], oneSixth);
1079 ep.Assign(4, cRing[corner.emOnBoundary ? 0 : 3], 0.0f);
1080 assert(ep.GetSize() == 5);
1081 }
1082
1083 if (corner.emOnBoundary) {
1084 em.Assign(0, cIndex, twoThirds);
1085 em.Assign(1, cRing[3], oneThird);
1086 em.Assign(2, cRing[0], 0.0f);
1087 assert(em.GetSize() == 3);
1088 } else {
1089 em.Assign(0, cIndex, 0.5f);
1090 em.Assign(1, cRing[1], oneSixth);
1091 em.Assign(2, cRing[2], oneSixth);
1092 em.Assign(3, cRing[corner.epOnBoundary ? 0 : 3], oneSixth);
1093 em.Assign(4, cRing[corner.epOnBoundary ? 3 : 0], 0.0f);
1094 assert(em.GetSize() == 5);
1095 }
1096 }
1097 }
1098
1099 template <typename REAL>
1100 void
computeIrregularEdgePoints(int cIndex,Matrix & matrix,Weight * weightBuffer) const1101 GregoryTriConverter<REAL>::computeIrregularEdgePoints(int cIndex,
1102 Matrix & matrix, Weight *weightBuffer) const {
1103
1104 Point p (matrix, 5*cIndex + 0);
1105 Point ep(matrix, 5*cIndex + 1);
1106 Point em(matrix, 5*cIndex + 2);
1107
1108 //
1109 // The corner and edge points P, Ep and Em are completely determined
1110 // by the 1-ring of vertices around (and including) the corner vertex.
1111 // We combine full sets of coefficients for the vertex and its 1-ring.
1112 //
1113 CornerTopology const & corner = _corners[cIndex];
1114
1115 if (corner.isSharp) {
1116 //
1117 // The sharp case -- both interior and boundary...
1118 //
1119 p.Assign(0, cIndex, 1.0f);
1120 assert(p.GetSize() == 1);
1121
1122 // Approximating these for now, pending future investigation...
1123 ep.Assign(0, cIndex, (REAL) (2.0 / 3.0));
1124 ep.Assign(1, (cIndex+1) % 3, (REAL) (1.0 / 3.0));
1125 assert(ep.GetSize() == 2);
1126
1127 em.Assign(0, cIndex, (REAL) (2.0 / 3.0));
1128 em.Assign(1, (cIndex+2) % 3, (REAL) (1.0 / 3.0));
1129 assert(em.GetSize() == 2);
1130 } else if (! corner.isBoundary) {
1131 //
1132 // The irregular interior case:
1133 //
1134 computeIrregularInteriorEdgePoints(cIndex, p, ep, em, weightBuffer);
1135 } else if (corner.numFaces > 1) {
1136 //
1137 // The irregular boundary case:
1138 //
1139 computeIrregularBoundaryEdgePoints(cIndex, p, ep, em, weightBuffer);
1140 } else {
1141 //
1142 // The irregular/smooth corner case:
1143 //
1144 p.Assign(0, cIndex, (REAL) (4.0 / 6.0));
1145 p.Assign(1, (cIndex+1) % 3, (REAL) (1.0 / 6.0));
1146 p.Assign(2, (cIndex+2) % 3, (REAL) (1.0 / 6.0));
1147 assert(p.GetSize() == 3);
1148
1149 ep.Assign(0, cIndex, (REAL) (2.0 / 3.0));
1150 ep.Assign(1, (cIndex+1) % 3, (REAL) (1.0 / 3.0));
1151 ep.Assign(2, (cIndex+2) % 3, 0.0f);
1152 assert(ep.GetSize() == 3);
1153
1154 em.Assign(0, cIndex, (REAL) (2.0 / 3.0));
1155 em.Assign(1, (cIndex+2) % 3, (REAL) (1.0 / 3.0));
1156 em.Assign(2, (cIndex+1) % 3, 0.0f);
1157 assert(em.GetSize() == 3);
1158 }
1159 }
1160
1161
1162 template <typename REAL>
1163 void
computeIrregularInteriorEdgePoints(int cIndex,Point & p,Point & ep,Point & em,Weight * ringWeights) const1164 GregoryTriConverter<REAL>::computeIrregularInteriorEdgePoints(
1165 int cIndex,
1166 Point& p, Point& ep, Point& em,
1167 Weight *ringWeights) const {
1168
1169 CornerTopology const & corner = _corners[cIndex];
1170
1171 int valence = corner.valence;
1172 int weightWidth = 1 + valence;
1173
1174 Weight* pWeights = &ringWeights[0];
1175 Weight* epWeights = pWeights + weightWidth;
1176 Weight* emWeights = pWeights + weightWidth * 2;
1177
1178 //
1179 // The interior (smooth) case -- invoke the public static method that
1180 // computes pre-allocated ring weights for P, Ep and Em:
1181 //
1182 LoopLimits<REAL>::ComputeInteriorPointWeights(valence, corner.faceInRing,
1183 pWeights, epWeights, emWeights);
1184
1185 //
1186 // Transer the weights for the ring into the stencil form of the required
1187 // Point type. The limit mask for position involves all ring weights, and
1188 // since Ep and Em depend on it, there should be no need to filter weights
1189 // with value 0:
1190 //
1191 p.Assign( 0, cIndex, pWeights[0]);
1192 ep.Assign(0, cIndex, epWeights[0]);
1193 em.Assign(0, cIndex, emWeights[0]);
1194
1195 for (int i = 1; i < weightWidth; ++i) {
1196 int pRingPoint = corner.ringPoints[i-1];
1197
1198 p.Assign( i, pRingPoint, pWeights[i]);
1199 ep.Assign(i, pRingPoint, epWeights[i]);
1200 em.Assign(i, pRingPoint, emWeights[i]);
1201 }
1202 assert(p.GetSize() == weightWidth);
1203 assert(ep.GetSize() == weightWidth);
1204 assert(em.GetSize() == weightWidth);
1205 }
1206
1207
1208 template <typename REAL>
1209 void
computeIrregularBoundaryEdgePoints(int cIndex,Point & p,Point & ep,Point & em,Weight * ringWeights) const1210 GregoryTriConverter<REAL>::computeIrregularBoundaryEdgePoints(
1211 int cIndex,
1212 Point& p, Point& ep, Point& em,
1213 Weight *ringWeights) const {
1214
1215 CornerTopology const & corner = _corners[cIndex];
1216
1217 int valence = corner.valence;
1218 int weightWidth = 1 + corner.valence;
1219
1220 Weight* pWeights = &ringWeights[0];
1221 Weight* epWeights = pWeights + weightWidth;
1222 Weight* emWeights = pWeights + weightWidth * 2;
1223
1224 //
1225 // The boundary (smooth) case -- invoke the public static method that
1226 // computes pre-allocated ring weights for P, Ep and Em:
1227 //
1228 LoopLimits<REAL>::ComputeBoundaryPointWeights(valence, corner.faceInRing,
1229 pWeights, epWeights, emWeights);
1230
1231 //
1232 // Transfer ring weights into points -- exploiting cases where they
1233 // are known to be non-zero only along the two boundary edges:
1234 //
1235 int N = weightWidth - 1;
1236
1237 int p0 = cIndex;
1238 int p1 = corner.ringPoints[0];
1239 int pN = corner.ringPoints[valence-1];
1240
1241 p.Assign(0, p0, pWeights[0]);
1242 p.Assign(1, p1, pWeights[1]);
1243 p.Assign(2, pN, pWeights[N]);
1244 assert(p.GetSize() == 3);
1245
1246 // If Ep is on the boundary edge, it has only two non-zero weights along
1247 // that edge:
1248 ep.Assign(0, p0, epWeights[0]);
1249 if (corner.epOnBoundary) {
1250 ep.Assign(1, p1, epWeights[1]);
1251 ep.Assign(2, pN, 0.0f);
1252 assert(ep.GetSize() == 3);
1253 } else {
1254 for (int i = 1; i < weightWidth; ++i) {
1255 ep.Assign(i, corner.ringPoints[i-1], epWeights[i]);
1256 }
1257 assert(ep.GetSize() == weightWidth);
1258 }
1259
1260 // If Em is on the boundary edge, it has only two non-zero weights along
1261 // that edge:
1262 em.Assign(0, p0, emWeights[0]);
1263 if (corner.emOnBoundary) {
1264 em.Assign(1, pN, emWeights[N]);
1265 em.Assign(2, p1, 0.0f);
1266 assert(em.GetSize() == 3);
1267 } else {
1268 for (int i = 1; i < weightWidth; ++i) {
1269 em.Assign(i, corner.ringPoints[i-1], emWeights[i]);
1270 }
1271 assert(em.GetSize() == weightWidth);
1272 }
1273 }
1274
1275
1276 template <typename REAL>
1277 int
getIrregularFacePointSize(int cIndexNear,int cIndexFar) const1278 GregoryTriConverter<REAL>::getIrregularFacePointSize(
1279 int cIndexNear, int cIndexFar) const {
1280
1281 CornerTopology const & nearCorner = _corners[cIndexNear];
1282 CornerTopology const & farCorner = _corners[cIndexFar];
1283
1284 if (nearCorner.isSharp && farCorner.isSharp) return 2;
1285
1286 int nearSize = nearCorner.ringPoints.GetSize() - 3;
1287 int farSize = farCorner.ringPoints.GetSize() - 3;
1288
1289 return 4 + (((nearSize > 0) && !nearCorner.isSharp) ? nearSize : 0)
1290 + (((farSize > 0) && !farCorner.isSharp) ? farSize : 0);
1291 }
1292
1293 template <typename REAL>
1294 void
computeIrregularFacePoint(int cIndexNear,int edgeInNearCornerRing,int cIndexFar,Point const & p,Point const & eNear,Point const & eFar,Point & fNear,REAL signForSideOfEdge,Weight * rowWeights,int * columnMask) const1295 GregoryTriConverter<REAL>::computeIrregularFacePoint(
1296 int cIndexNear, int edgeInNearCornerRing, int cIndexFar,
1297 Point const & p, Point const & eNear, Point const & eFar, Point & fNear,
1298 REAL signForSideOfEdge, Weight *rowWeights, int *columnMask) const {
1299
1300 CornerTopology const & cornerNear = _corners[cIndexNear];
1301 CornerTopology const & cornerFar = _corners[cIndexFar];
1302
1303 int valence = cornerNear.valence;
1304
1305 Weight cosNear = cornerNear.cosFaceAngle;
1306 Weight cosFar = cornerFar.cosFaceAngle;
1307
1308 //
1309 // From Loop, Schaefer et al, a face point F is computed as:
1310 //
1311 // F = (1/d) * (c0 * P0 + (d - 2*c0 - c1) * E0 + 2*c1 * E1 + R)
1312 //
1313 // where d = 3 for quads and d = 4 for triangles, and R is:
1314 //
1315 // R = 1/3 (Mm + Mp) + 2/3 (Cm + Cp)
1316 //
1317 // where Mm and Mp are the midpoints of the two adjacent edges
1318 // and Cm and Cp are the midpoints of the two adjacent faces.
1319 //
1320 Weight pCoeff = cosFar / 4.0f;
1321 Weight eNearCoeff = (4.0f - 2.0f * cosNear - cosFar) / 4.0f;
1322 Weight eFarCoeff = 2.0f * cosNear / 4.0f;
1323
1324 int fullRowSize = _numSourcePoints;
1325 std::memset(&columnMask[0], 0, fullRowSize * sizeof(int));
1326 std::memset(&rowWeights[0], 0, fullRowSize * sizeof(Weight));
1327
1328 _addSparsePointToFullRow(rowWeights, p, pCoeff, columnMask);
1329 _addSparsePointToFullRow(rowWeights, eNear, eNearCoeff, columnMask);
1330 _addSparsePointToFullRow(rowWeights, eFar, eFarCoeff, columnMask);
1331
1332 // Remember that R is to be computed about an interior edge and is
1333 // comprised of the two pairs of points opposite the interior edge
1334 //
1335 int iEdgeInterior = edgeInNearCornerRing;
1336 int iEdgePrev = (iEdgeInterior + valence - 1) % valence;
1337 int iEdgeNext = (iEdgeInterior + 1) % valence;
1338
1339 Weight rScale = (REAL) (0.25 * (7.0 / 18.0));
1340
1341 rowWeights[cornerNear.ringPoints[iEdgePrev]] += -signForSideOfEdge * rScale;
1342 rowWeights[cornerNear.ringPoints[iEdgeNext]] += signForSideOfEdge * rScale;
1343
1344 int nWeights = 0;
1345 for (int i = 0; i < fullRowSize; ++i) {
1346 if (columnMask[i]) {
1347 fNear.Assign(nWeights++, columnMask[i] - 1, rowWeights[i]);
1348 }
1349 }
1350
1351 // Complete the expected row size when val-2 corners induce duplicates:
1352 if (_hasVal2InteriorCorner && (nWeights < fNear.GetSize())) {
1353 while (nWeights < fNear.GetSize()) {
1354 fNear.Assign(nWeights++, cIndexNear, 0.0f);
1355 }
1356 }
1357 assert(fNear.GetSize() == nWeights);
1358 }
1359
1360 template <typename REAL>
1361 void
assignRegularFacePoints(int cIndex,Matrix & matrix) const1362 GregoryTriConverter<REAL>::assignRegularFacePoints(int cIndex, Matrix & matrix) const {
1363
1364 CornerTopology const & corner = _corners[cIndex];
1365
1366 int cNext = (cIndex+1) % 3;
1367 int cPrev = (cIndex+2) % 3;
1368
1369 int const * cRing = corner.ringPoints;
1370
1371 //
1372 // Regular face-points are computed the same for both face-points of a
1373 // a corner (fp and fm), so iterate through both and make appropriate
1374 // assignments when tagged as regular:
1375 //
1376 for (int fIsFm = 0; fIsFm < 2; ++fIsFm) {
1377 bool fIsRegular = fIsFm ? corner.fmIsRegular : corner.fpIsRegular;
1378 if (!fIsRegular) continue;
1379
1380 Point f(matrix, 5*cIndex + 3 + fIsFm);
1381
1382 if (corner.isCorner) {
1383 f.Assign(0, cIndex, 0.5f);
1384 f.Assign(1, cNext, 0.25f);
1385 f.Assign(2, cPrev, 0.25f);
1386 assert(f.GetSize() == 3);
1387 } else if (corner.epOnBoundary) {
1388 // Face is the first/leading face of the boundary ring:
1389 f.Assign(0, cIndex, (REAL) (11.0 / 24.0));
1390 f.Assign(1, cRing[0], (REAL) ( 7.0 / 24.0));
1391 f.Assign(2, cRing[1], (REAL) ( 5.0 / 24.0));
1392 f.Assign(3, cRing[2], (REAL) ( 1.0 / 24.0));
1393 assert(f.GetSize() == 4);
1394 } else if (corner.emOnBoundary) {
1395 // Face is the last/trailing face of the boundary ring:
1396 f.Assign(0, cIndex, (REAL) (11.0 / 24.0));
1397 f.Assign(1, cRing[3], (REAL) ( 7.0 / 24.0));
1398 f.Assign(2, cRing[2], (REAL) ( 5.0 / 24.0));
1399 f.Assign(3, cRing[1], (REAL) ( 1.0 / 24.0));
1400 assert(f.GetSize() == 4);
1401 } else {
1402 // Face is interior or the middle face of the boundary:
1403 int eNext = corner.isBoundary ? 0 : ((corner.faceInRing + 5) % 6);
1404 int ePrev = corner.isBoundary ? 3 : ((corner.faceInRing + 2) % 6);
1405
1406 f.Assign(0, cIndex, (REAL) (10.0 / 24.0));
1407 f.Assign(1, cPrev, 0.25f);
1408 f.Assign(2, cNext, 0.25f);
1409 f.Assign(3, cRing[ePrev], (REAL) ( 1.0 / 24.0));
1410 f.Assign(4, cRing[eNext], (REAL) ( 1.0 / 24.0));
1411 assert(f.GetSize() == 5);
1412 }
1413 }
1414 }
1415
1416 template <typename REAL>
1417 void
computeIrregularFacePoints(int cIndex,Matrix & matrix,Weight * rowWeights,int * columnMask) const1418 GregoryTriConverter<REAL>::computeIrregularFacePoints(int cIndex,
1419 Matrix & matrix, Weight *rowWeights, int *columnMask) const {
1420
1421 // Identify neighboring corners:
1422 CornerTopology const & corner = _corners[cIndex];
1423
1424 int cNext = (cIndex+1) % 3;
1425 int cPrev = (cIndex+2) % 3;
1426
1427 Point epPrev(matrix, 5*cPrev + 1);
1428 Point em (matrix, 5*cIndex + 2);
1429 Point p (matrix, 5*cIndex + 0);
1430 Point ep (matrix, 5*cIndex + 1);
1431 Point emNext(matrix, 5*cNext + 2);
1432
1433 Point fp(matrix, 5*cIndex + 3);
1434 Point fm(matrix, 5*cIndex + 4);
1435
1436 //
1437 // Compute the face points Fp and Fm in terms of the corner (P) and edge
1438 // points (Ep and Em) previously computed. The caller provides a buffer
1439 // of the appropriate size (twice the width of the matrix) to use for
1440 // combining weights, along with an integer buffer used to identify
1441 // non-zero weights and preserve the sparsity of the combinations (note
1442 // they use index + 1 to detect index 0 when cleared with 0 entries).
1443 //
1444 if (!corner.fpIsRegular && !corner.fpIsCopied) {
1445 int iEdgeP = corner.faceInRing;
1446 computeIrregularFacePoint(cIndex, iEdgeP, cNext,
1447 p, ep, emNext, fp, 1.0, rowWeights, columnMask);
1448 }
1449 if (!corner.fmIsRegular && !corner.fmIsCopied) {
1450 int iEdgeM = (corner.faceInRing + 1) % corner.valence;
1451 computeIrregularFacePoint(cIndex, iEdgeM, cPrev,
1452 p, em, epPrev, fm, -1.0, rowWeights, columnMask);
1453 }
1454
1455 // Copy Fp or Fm now that any shared values were computed above:
1456 if (corner.fpIsCopied) {
1457 fp.Copy(fm);
1458 }
1459 if (corner.fmIsCopied) {
1460 fm.Copy(fp);
1461 }
1462
1463 if (!corner.fpIsRegular) assert(matrix.GetRowSize(5*cIndex + 3) == fp.GetSize());
1464 if (!corner.fmIsRegular) assert(matrix.GetRowSize(5*cIndex + 4) == fm.GetSize());
1465 }
1466
1467 template <typename REAL>
1468 void
assignRegularMidEdgePoint(int edgeIndex,Matrix & matrix) const1469 GregoryTriConverter<REAL>::assignRegularMidEdgePoint(int edgeIndex,
1470 Matrix & matrix) const {
1471
1472 Point M(matrix, 15 + edgeIndex);
1473
1474 CornerTopology const & corner = _corners[edgeIndex];
1475 if (corner.epOnBoundary) {
1476 // Trivial boundary edge case -- midway between two corners
1477
1478 M.Assign(0, edgeIndex, 0.5f);
1479 M.Assign(1, (edgeIndex + 1) % 3, 0.5f);
1480 assert(M.GetSize() == 2);
1481 } else {
1482 // Regular case -- two corners and two vertices opposite the edge
1483
1484 int oppositeInRing = corner.isBoundary ?
1485 (corner.faceInRing - 1) : ((corner.faceInRing + 5) % 6);
1486 int oppositeVertex = corner.ringPoints[oppositeInRing];
1487
1488 M.Assign(0, edgeIndex, (REAL) (1.0 / 3.0));
1489 M.Assign(1, (edgeIndex + 1) % 3, (REAL) (1.0 / 3.0));
1490 M.Assign(2, (edgeIndex + 2) % 3, (REAL) (1.0 / 6.0));
1491 M.Assign(3, oppositeVertex, (REAL) (1.0 / 6.0));
1492 assert(M.GetSize() == 4);
1493 }
1494 }
1495
1496 template <typename REAL>
1497 void
computeIrregularMidEdgePoint(int edgeIndex,Matrix & matrix,Weight * rowWeights,int * columnMask) const1498 GregoryTriConverter<REAL>::computeIrregularMidEdgePoint(int edgeIndex, Matrix & matrix,
1499 Weight * rowWeights, int * columnMask) const {
1500 //
1501 // General case -- interpolate midway between cubic edge points E0 and E1:
1502 //
1503 int cIndex0 = edgeIndex;
1504 int cIndex1 = (edgeIndex + 1) % 3;
1505
1506 Point E0p(matrix, 5 * (cIndex0) + 1);
1507 Point E1m(matrix, 5 * (cIndex1) + 2);
1508
1509 Point M(matrix, 15 + edgeIndex);
1510
1511 _combineSparsePointsInFullRow(M, (REAL)0.5f, E0p, (REAL)0.5f, E1m,
1512 _numSourcePoints, rowWeights, columnMask);
1513 }
1514
1515 template <typename REAL>
1516 void
promoteCubicEdgePointsToQuartic(Matrix & matrix,Weight * rowWeights,int * columnMask) const1517 GregoryTriConverter<REAL>::promoteCubicEdgePointsToQuartic(Matrix & matrix,
1518 Weight * rowWeights, int * columnMask) const {
1519 //
1520 // Re-assign all regular edge-point weights with quartic coefficients,
1521 // so only perform general combinations for the irregular case.
1522 //
1523 REAL const onBoundaryWeights[3] = { 16, 7, 1 };
1524 REAL const regBoundaryWeights[5] = { 13, 3, 3, 4, 1 };
1525 REAL const regInteriorWeights[7] = { 12, 4, 3, 1, 0, 1, 3 };
1526
1527 REAL const oneOver24 = (REAL) (1.0 / 24.0);
1528
1529 for (int cIndex = 0; cIndex < 3; ++cIndex) {
1530 CornerTopology const & corner = _corners[cIndex];
1531
1532 //
1533 // Ordering of weight values for symmetric ep and em is the same, so
1534 // we can re-assign in a loop of 2 for {ep, em}
1535 //
1536 Point P(matrix, 5 * cIndex);
1537
1538 for (int ePair = 0; ePair < 2; ++ePair) {
1539 Point E(matrix, 5 * cIndex + 1 + ePair);
1540
1541 REAL const * weightsToReassign = 0;
1542
1543 bool eOnBoundary = ePair ? corner.emOnBoundary : corner.epOnBoundary;
1544 if (eOnBoundary && !corner.isSharp) {
1545 assert(E.GetSize() == 3);
1546 weightsToReassign = onBoundaryWeights;
1547 } else if (corner.isRegular) {
1548 if (corner.isBoundary) {
1549 assert(E.GetSize() == 5);
1550 weightsToReassign = regBoundaryWeights;
1551 } else {
1552 assert(E.GetSize() == 7);
1553 weightsToReassign = regInteriorWeights;
1554 }
1555 }
1556 if (weightsToReassign) {
1557 for (int i = 0; i < E.GetSize(); ++i) {
1558 E.SetWeight(i, weightsToReassign[i] * oneOver24);
1559 }
1560 } else {
1561 _combineSparsePointsInFullRow(E, (REAL)0.25f, P, (REAL)0.75f, E,
1562 _numSourcePoints, rowWeights, columnMask);
1563 }
1564 }
1565 }
1566 }
1567
1568 namespace {
1569 template <typename REAL>
1570 void
_printPoint(SparseMatrixRow<REAL> const & p,bool printIndices=true,bool printWeights=true)1571 _printPoint(SparseMatrixRow<REAL> const & p, bool printIndices = true,
1572 bool printWeights = true) {
1573 printf(" Point size = %d:\n", p._size);
1574 if (printIndices) {
1575 printf(" Indices: ");
1576 for (int j = 0; j < p._size; ++j) {
1577 printf("%6d ", p._indices[j]);
1578 }
1579 printf("\n");
1580 }
1581 if (printWeights) {
1582 printf(" Weights: ");
1583 for (int j = 0; j < p._size; ++j) {
1584 printf("%6.3f ", (REAL)p._weights[j]);
1585 }
1586 printf("\n");
1587 }
1588 }
1589
1590 template <typename REAL>
1591 void
_printMatrix(SparseMatrix<REAL> const & matrix,bool printIndices=true,bool printWeights=true)1592 _printMatrix(SparseMatrix<REAL> const & matrix, bool printIndices = true,
1593 bool printWeights = true) {
1594
1595 printf("Matrix %d x %d, %d elements:\n",
1596 matrix.GetNumRows(), matrix.GetNumColumns(), matrix.GetNumElements());
1597
1598 for (int i = 0; i < matrix.GetNumRows(); ++i) {
1599 int rowSize = matrix.GetRowSize(i);
1600 printf(" Row %d (size = %d):\n", i, rowSize);
1601
1602 if (printIndices) {
1603 ConstArray<int> indices = matrix.GetRowColumns(i);
1604 printf(" Indices: ");
1605 for (int j = 0; j < rowSize; ++j) {
1606 printf("%6d ", indices[j]);
1607 }
1608 printf("\n");
1609 }
1610 if (printWeights) {
1611 ConstArray<REAL> weights = matrix.GetRowElements(i);
1612 printf(" Weights: ");
1613 for (int j = 0; j < rowSize; ++j) {
1614 printf("%6.3f ", (REAL)weights[j]);
1615 }
1616 printf("\n");
1617 }
1618 }
1619 }
1620
1621 void
_printSourcePatch(SourcePatch const & patch,bool printCornerInfo=true,bool printRingPoints=true)1622 _printSourcePatch(SourcePatch const & patch, bool printCornerInfo = true,
1623 bool printRingPoints = true) {
1624
1625 printf("SoucePatch - %d corners, %d points:\n",
1626 patch._numCorners, patch._numSourcePoints);
1627
1628 if (printCornerInfo) {
1629 printf(" Corner info:\n");
1630 for (int i = 0; i < patch._numCorners; ++i) {
1631 printf("%6d: boundary = %d, sharp = %d, numFaces = %d, in-ring = %d, ringSize = %d\n", i,
1632 patch._corners[i]._boundary, patch._corners[i]._sharp, patch._corners[i]._numFaces,
1633 patch._corners[i]._patchFace, patch._ringSizes[i]);
1634 }
1635 }
1636 if (printRingPoints) {
1637 StackBuffer<Index,64,true> ringPoints;
1638 printf(" Ring points:\n");
1639 for (int i = 0; i < patch._numCorners; ++i) {
1640 int ringSize = patch._ringSizes[i];
1641
1642 ringPoints.SetSize(ringSize);
1643 patch.GetCornerRingPoints(i, ringPoints);
1644
1645 printf("%6d: ", i);
1646 for (int j = 0; j < ringSize; ++j) {
1647 printf("%d ", ringPoints[j]);
1648 }
1649 printf("\n");
1650 }
1651 }
1652 }
1653 }
1654
1655
1656 //
1657 // Not using the same "Converter" pattern used above and in the Catmark scheme,
1658 // since the two remaining conversions are fairly trivial.
1659 //
1660 template <typename REAL>
1661 void
convertToLinear(SourcePatch const & sourcePatch,SparseMatrix<REAL> & matrix)1662 convertToLinear(SourcePatch const & sourcePatch, SparseMatrix<REAL> & matrix) {
1663
1664 typedef REAL Weight;
1665
1666 StackBuffer<Index,64,true> indexBuffer(1 + sourcePatch.GetMaxRingSize());
1667 StackBuffer<Weight,64,true> weightBuffer(1 + sourcePatch.GetMaxRingSize());
1668
1669 int numElements = sourcePatch.GetCornerRingSize(0) +
1670 sourcePatch.GetCornerRingSize(1) +
1671 sourcePatch.GetCornerRingSize(2);
1672
1673 matrix.Resize(3, sourcePatch.GetNumSourcePoints(), numElements);
1674
1675 bool hasVal2InteriorCorner = false;
1676
1677 for (int cIndex = 0; cIndex < 3; ++cIndex) {
1678 SourcePatch::Corner const & sourceCorner = sourcePatch._corners[cIndex];
1679
1680 int ringSize = sourcePatch.GetCornerRingSize(cIndex);
1681 if (sourceCorner._sharp) {
1682 matrix.SetRowSize(cIndex, 1);
1683 } else if (sourceCorner._boundary) {
1684 matrix.SetRowSize(cIndex, 3);
1685 } else {
1686 matrix.SetRowSize(cIndex, 1 + ringSize);
1687 }
1688
1689 Array<Index> rowIndices = matrix.SetRowColumns(cIndex);
1690 Array<Weight> rowWeights = matrix.SetRowElements(cIndex);
1691
1692 indexBuffer[0] = cIndex;
1693 sourcePatch.GetCornerRingPoints(cIndex, &indexBuffer[1]);
1694
1695 if (sourceCorner._sharp) {
1696 rowIndices[0] = cIndex;
1697 rowWeights[0] = 1.0f;
1698 } else if (sourceCorner._boundary) {
1699 LoopLimits<REAL>::ComputeBoundaryPointWeights(
1700 1 + sourceCorner._numFaces, sourceCorner._patchFace,
1701 &weightBuffer[0], 0, 0);
1702
1703 rowIndices[0] = indexBuffer[0];
1704 rowIndices[1] = indexBuffer[1];
1705 rowIndices[2] = indexBuffer[ringSize];
1706
1707 rowWeights[0] = weightBuffer[0];
1708 rowWeights[1] = weightBuffer[1];
1709 rowWeights[2] = weightBuffer[ringSize];
1710 } else {
1711 LoopLimits<REAL>::ComputeInteriorPointWeights(
1712 sourceCorner._numFaces, sourceCorner._patchFace,
1713 &weightBuffer[0], 0, 0);
1714
1715 std::memcpy(&rowIndices[0], indexBuffer, rowIndices.size() * sizeof(Index));
1716 std::memcpy(&rowWeights[0], weightBuffer, rowWeights.size() * sizeof(Weight));
1717 }
1718 hasVal2InteriorCorner |= sourceCorner._val2Interior;
1719 }
1720 if (hasVal2InteriorCorner) {
1721 _removeValence2Duplicates(matrix);
1722 }
1723 }
1724
1725 template <typename REAL>
1726 void
convertToGregory(SourcePatch const & sourcePatch,SparseMatrix<REAL> & matrix)1727 convertToGregory(SourcePatch const & sourcePatch, SparseMatrix<REAL> & matrix) {
1728
1729 GregoryTriConverter<REAL> gregoryConverter(sourcePatch);
1730 gregoryConverter.Convert(matrix);
1731 }
1732
1733 template <typename REAL>
1734 void
convertToLoop(SourcePatch const & sourcePatch,SparseMatrix<REAL> & matrix)1735 convertToLoop(SourcePatch const & sourcePatch, SparseMatrix<REAL> & matrix) {
1736
1737 //
1738 // Unlike quads, there are not enough degrees of freedom in the regular patch
1739 // to enforce interpolation of the limit point and tangent at the corner while
1740 // preserving the two adjoining boundary curves. Since we end up destroying
1741 // neighboring conintuity in doing so, we use a fully constructed Gregory
1742 // patch here for the isolated corner case as well as the general case.
1743 //
1744 // Unfortunately, the regular patch here -- a quartic Box-spline triangle --
1745 // is not as flexible as the BSpline patches for Catmark. Unlike BSplines
1746 // and Bezier patches, the Box-splines do not span the full space of possible
1747 // shapes (only 12 control points in a space spanned by 15 polynomials for
1748 // the quartic case). So it is possible to construct shapes with a Gregory
1749 // or Bezier triangle that cannot be represented by the Box-spline.
1750 //
1751 // The solution fits a Box-spline patch to the constructed Gregory patch with
1752 // a set of constraints. With quartic boundary curves, 12 constraints on the
1753 // boundary curve make this tightly constrained. Such a set of constraints
1754 // is rank deficient (11 instead of 12) so an additional constraint on the
1755 // midpoint of the patch is included and a conversion matrix is constructed
1756 // from the pseudo-inverse of the 13 constraints.
1757 //
1758 // For the full 12x15 conversion matrix from 15-point quartic Bezier patch
1759 // back to a Box spline patch, the matrix rows and columns are ordered
1760 // according to control point orientations used elsewhere. Correllation of
1761 // points between the Bezier and Gregory points is as follows:
1762 //
1763 // Q0 Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11 Q12 Q13 Q14
1764 // G0 G1 G15 G7 G5 G2 G3,4 G8,9 G6 G17 G13,14 G16 G11 G12 G10
1765 //
1766 // As with conversion from Gregory to BSpline for Catmark, one of the face
1767 // points is chosen as a Bezier point in the conversion rather than combining
1768 // the pair (which would avoid slight asymmetric artefacts of the choice).
1769 // And given the solution now depends primarily on the boundary, its not
1770 // necessary to construct a full Gregory patch with enforced continuity.
1771 //
1772 REAL const gregoryToLoopMatrix[12][15] = {
1773 { 8.214411f, 7.571190f, -7.690082f, 2.237840f, -1.118922f,-16.428828f, 0.666666f, 0.666666f,
1774 2.237835f, 6.309870f, 0.666666f, -1.690100f, -0.428812f, -0.428805f, 0.214407f },
1775 { -0.304687f, 0.609374f, 6.752593f, 0.609374f, -0.304687f, 0.609378f, -3.333333f, -3.333333f,
1776 0.609378f, -1.247389f, -3.333333f, -1.247389f, 3.276037f, 3.276037f, -1.638020f },
1777 { -1.118922f, 2.237840f, -7.690082f, 7.571190f, 8.214411f, 2.237835f, 0.666666f, 0.666666f,
1778 -16.428828f, -1.690100f, 0.666666f, 6.309870f, -0.428805f, -0.428812f, 0.214407f },
1779 { 8.214411f,-16.428828f, 6.309870f, -0.428812f, 0.214407f, 7.571190f, 0.666666f, 0.666666f,
1780 -0.428805f, -7.690082f, 0.666666f, -1.690100f, 2.237840f, 2.237835f, -1.118922f },
1781 { -0.813368f, 1.626735f, -0.773435f, -1.039929f, 0.519965f, 1.626735f, 0.666666f, 0.666666f,
1782 -1.039930f, -0.773435f, 0.666666f, 1.226558f, -1.039929f, -1.039930f, 0.519965f },
1783 { 0.519965f, -1.039929f, -0.773435f, 1.626735f, -0.813368f, -1.039930f, 0.666666f, 0.666666f,
1784 1.626735f, 1.226558f, 0.666666f, -0.773435f, -1.039930f, -1.039929f, 0.519965f },
1785 { 0.214407f, -0.428812f, 6.309870f,-16.428828f, 8.214411f, -0.428805f, 0.666666f, 0.666666f,
1786 7.571190f, -1.690100f, 0.666666f, -7.690082f, 2.237835f, 2.237840f, -1.118922f },
1787 { -0.304687f, 0.609378f, -1.247389f, 3.276037f, -1.638020f, 0.609374f, -3.333333f, -3.333333f,
1788 3.276037f, 6.752593f, -3.333333f, -1.247389f, 0.609374f, 0.609378f, -0.304687f },
1789 { 0.519965f, -1.039930f, 1.226558f, -1.039930f, 0.519965f, -1.039929f, 0.666666f, 0.666666f,
1790 -1.039929f, -0.773435f, 0.666666f, -0.773435f, 1.626735f, 1.626735f, -0.813368f },
1791 { -1.638020f, 3.276037f, -1.247389f, 0.609378f, -0.304687f, 3.276037f, -3.333333f, -3.333333f,
1792 0.609374f, -1.247389f, -3.333333f, 6.752593f, 0.609378f, 0.609374f, -0.304687f },
1793 { -1.118922f, 2.237835f, -1.690100f, -0.428805f, 0.214407f, 2.237840f, 0.666666f, 0.666666f,
1794 -0.428812f, -7.690082f, 0.666666f, 6.309870f, 7.571190f,-16.428828f, 8.214411f },
1795 { 0.214407f, -0.428805f, -1.690100f, 2.237835f, -1.118922f, -0.428812f, 0.666666f, 0.666666f,
1796 2.237840f, 6.309870f, 0.666666f, -7.690082f,-16.428828f, 7.571190f, 8.214411f }
1797 };
1798 int const gRowIndices[15] = { 0,1,15,7,5, 2,4,8,6, 17,14,16, 11,12, 10 };
1799
1800 SparseMatrix<REAL> G;
1801 convertToGregory<REAL>(sourcePatch, G);
1802
1803 _initializeFullMatrix(matrix, 12, G.GetNumColumns());
1804
1805 for (int i = 0; i < 12; ++i) {
1806 REAL const * gRowWeights = gregoryToLoopMatrix[i];
1807 _combineSparseMatrixRowsInFull(matrix, i, G, 15, gRowIndices, gRowWeights);
1808 }
1809 }
1810
1811
1812 //
1813 // Internal utilities more relevant to the LoopPatchBuilder:
1814 //
1815 namespace {
1816 //
1817 // The patch type associated with each basis for Loop -- quickly
1818 // indexed from an array. The patch type here is essentially the
1819 // triangle form of each basis.
1820 //
1821 PatchDescriptor::Type patchTypeFromBasisArray[] = {
1822 PatchDescriptor::NON_PATCH, // undefined
1823 PatchDescriptor::LOOP, // regular
1824 PatchDescriptor::GREGORY_TRIANGLE, // Gregory
1825 PatchDescriptor::TRIANGLES, // linear
1826 PatchDescriptor::NON_PATCH }; // Bezier -- for future use
1827 };
1828
LoopPatchBuilder(TopologyRefiner const & refiner,Options const & options)1829 LoopPatchBuilder::LoopPatchBuilder(
1830 TopologyRefiner const& refiner, Options const& options) :
1831 PatchBuilder(refiner, options) {
1832
1833 _regPatchType = patchTypeFromBasisArray[_options.regBasisType];
1834 _irregPatchType = (_options.irregBasisType == BASIS_UNSPECIFIED)
1835 ? _regPatchType
1836 : patchTypeFromBasisArray[_options.irregBasisType];
1837
1838 _nativePatchType = PatchDescriptor::LOOP;
1839 _linearPatchType = PatchDescriptor::TRIANGLES;
1840 }
1841
~LoopPatchBuilder()1842 LoopPatchBuilder::~LoopPatchBuilder() {
1843 }
1844
1845 PatchDescriptor::Type
patchTypeFromBasis(BasisType basis) const1846 LoopPatchBuilder::patchTypeFromBasis(BasisType basis) const {
1847
1848 return patchTypeFromBasisArray[(int)basis];
1849 }
1850
1851 template <typename REAL>
1852 int
convertSourcePatch(SourcePatch const & sourcePatch,PatchDescriptor::Type patchType,SparseMatrix<REAL> & matrix) const1853 LoopPatchBuilder::convertSourcePatch(SourcePatch const & sourcePatch,
1854 PatchDescriptor::Type patchType,
1855 SparseMatrix<REAL> & matrix) const {
1856
1857 assert(_schemeType == Sdc::SCHEME_LOOP);
1858
1859 if (patchType == PatchDescriptor::LOOP) {
1860 convertToLoop<REAL>(sourcePatch, matrix);
1861 } else if (patchType == PatchDescriptor::TRIANGLES) {
1862 convertToLinear<REAL>(sourcePatch, matrix);
1863 } else if (patchType == PatchDescriptor::GREGORY_TRIANGLE) {
1864 convertToGregory<REAL>(sourcePatch, matrix);
1865 } else {
1866 assert("Unknown or unsupported patch type" == 0);
1867 }
1868 return matrix.GetNumRows();
1869 }
1870
1871 int
convertToPatchType(SourcePatch const & sourcePatch,PatchDescriptor::Type patchType,SparseMatrix<float> & matrix) const1872 LoopPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
1873 PatchDescriptor::Type patchType,
1874 SparseMatrix<float> & matrix) const {
1875 return convertSourcePatch(sourcePatch, patchType, matrix);
1876 }
1877 int
convertToPatchType(SourcePatch const & sourcePatch,PatchDescriptor::Type patchType,SparseMatrix<double> & matrix) const1878 LoopPatchBuilder::convertToPatchType(SourcePatch const & sourcePatch,
1879 PatchDescriptor::Type patchType,
1880 SparseMatrix<double> & matrix) const {
1881 return convertSourcePatch(sourcePatch, patchType, matrix);
1882 }
1883
1884
1885 } // end namespace Far
1886
1887 } // end namespace OPENSUBDIV_VERSION
1888 } // end namespace OpenSubdiv
1889