text: //
  //   Copyright 2014 DreamWorks Animation LLC.
  //
  //   Licensed under the Apache License, Version 2.0 (the "Apache License")
  //   with the following modification; you may not use this file except in
  //   compliance with the Apache License and the following modification to it:
  //   Section 6. Trademarks. is deleted and replaced with:
  //
  //   6. Trademarks. This License does not grant permission to use the trade
  //      names, trademarks, service marks, or product names of the Licensor
  //      and its affiliates, except as required to comply with Section 4(c) of
  //      the License and to reproduce the content of the NOTICE file.
  //
  //   You may obtain a copy of the Apache License at
  //
  //       http://www.apache.org/licenses/LICENSE-2.0
  //
  //   Unless required by applicable law or agreed to in writing, software
  //   distributed under the Apache License with the above modification is
  //   distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
  //   KIND, either express or implied. See the Apache License for the specific
  //   language governing permissions and limitations under the Apache License.
  //
  #ifndef OPENSUBDIV3_SDC_CATMARK_SCHEME_H
  #define OPENSUBDIV3_SDC_CATMARK_SCHEME_H
  
  #include "../version.h"
  
  #include "../sdc/scheme.h"
  
  #include <cassert>
  #include <cmath>
  
  namespace OpenSubdiv {
  namespace OPENSUBDIV_VERSION {
  
  namespace Sdc {
  
  //
  //  Specializations for Scheme<SCHEME_CATMARK>:
  //
  
  //
  //  Catmark traits:
  //
  template <>
  inline Split Scheme<SCHEME_CATMARK>::GetTopologicalSplitType() { return SPLIT_TO_QUADS; }
  
  template <>
  inline int Scheme<SCHEME_CATMARK>::GetRegularFaceSize() { return 4; }
  
  template <>
  inline int Scheme<SCHEME_CATMARK>::GetRegularVertexValence() { return 4; }
  
  template <>
  inline int Scheme<SCHEME_CATMARK>::GetLocalNeighborhoodSize() { return 1; }
  
  
  //
  //  Masks for edge-vertices:  the hard Crease mask does not need to be specialized
  //  (simply the midpoint), so all that is left is the Smooth case:
  //
  //  The Smooth mask is complicated by the need to support the "triangle subdivision"
  //  option, which applies different weighting in the presence of triangles.  It is
  //  up for debate as to whether this is useful or not -- we may be able to deprecate
  //  this option.
  //
  template <>
  template <typename EDGE, typename MASK>
  inline void
  Scheme<SCHEME_CATMARK>::assignSmoothMaskForEdge(EDGE const& edge, MASK& mask) const {
  
      typedef typename MASK::Weight Weight;
  
      int faceCount = edge.GetNumFaces();
  
      mask.SetNumVertexWeights(2);
      mask.SetNumEdgeWeights(0);
      mask.SetNumFaceWeights(faceCount);
      mask.SetFaceWeightsForFaceCenters(true);
  
      //
      //  Determine if we need to inspect incident faces and apply alternate weighting for
      //  triangles -- and if so, determine which of the two are triangles.
      //
      bool face0IsTri = false;
      bool face1IsTri = false;
      bool useTriangleOption = (_options.GetTriangleSubdivision() == Options::TRI_SUB_SMOOTH);
      if (useTriangleOption) {
          if (faceCount == 2) {
              //
              //  Ideally we want to avoid this inspection when we have already subdivided at
              //  least once -- need something in the Edge interface to help avoid this, e.g.
              //  an IsRegular() query, the subdivision level...
              //
              int vertsPerFace[2];
              edge.GetNumVerticesPerFace(vertsPerFace);
  
              face0IsTri = (vertsPerFace[0] == 3);
              face1IsTri = (vertsPerFace[1] == 3);
              useTriangleOption = face0IsTri || face1IsTri;
          } else {
              useTriangleOption = false;
          }
      }
  
      if (! useTriangleOption) {
          mask.VertexWeight(0) = 0.25f;
          mask.VertexWeight(1) = 0.25f;
  
          if (faceCount == 2) {
              mask.FaceWeight(0) = 0.25f;
              mask.FaceWeight(1) = 0.25f;
          } else {
              Weight fWeight = 0.5f / (Weight)faceCount;
              for (int i = 0; i < faceCount; ++i) {
                  mask.FaceWeight(i) = fWeight;
              }
          }
      } else {
          //
          //  This mimics the implementation in Hbr in terms of order of operations.
          //
          const Weight CATMARK_SMOOTH_TRI_EDGE_WEIGHT = 0.470f;
  
          Weight f0Weight = face0IsTri ? CATMARK_SMOOTH_TRI_EDGE_WEIGHT : 0.25f;
          Weight f1Weight = face1IsTri ? CATMARK_SMOOTH_TRI_EDGE_WEIGHT : 0.25f;
  
          Weight fWeight = 0.5f * (f0Weight + f1Weight);
          Weight vWeight = 0.5f * (1.0f - 2.0f * fWeight);
  
          mask.VertexWeight(0) = vWeight;
          mask.VertexWeight(1) = vWeight;
  
          mask.FaceWeight(0) = fWeight;
          mask.FaceWeight(1) = fWeight;
      }
  }
  
  
  //
  //  Masks for vertex-vertices:  the hard Corner mask does not need to be specialized
  //  (simply the vertex itself), leaving the Crease and Smooth cases (Dart is smooth):
  //
  template <>
  template <typename VERTEX, typename MASK>
  inline void
  Scheme<SCHEME_CATMARK>::assignCreaseMaskForVertex(VERTEX const& vertex, MASK& mask,
                                                    int const creaseEnds[2]) const {
      typedef typename MASK::Weight Weight;
  
      int valence = vertex.GetNumEdges();
  
      mask.SetNumVertexWeights(1);
      mask.SetNumEdgeWeights(valence);
      mask.SetNumFaceWeights(0);
      mask.SetFaceWeightsForFaceCenters(false);
  
      Weight vWeight = 0.75f;
      Weight eWeight = 0.125f;
  
      mask.VertexWeight(0) = vWeight;
      for (int i = 0; i < valence; ++i) {
          mask.EdgeWeight(i) = 0.0f;
      }
      mask.EdgeWeight(creaseEnds[0]) = eWeight;
      mask.EdgeWeight(creaseEnds[1]) = eWeight;
  }
  
  template <>
  template <typename VERTEX, typename MASK>
  inline void
  Scheme<SCHEME_CATMARK>::assignSmoothMaskForVertex(VERTEX const& vertex, MASK& mask) const {
  
      typedef typename MASK::Weight Weight;
  
      //
      //  A Smooth vertex must be manifold and interior -- manifold boundary vertices will be
      //  Creases and non-manifold vertices of any kind will be Corners or Creases.  If smooth
      //  rules for non-manifold vertices are ever defined, this will need adjusting:
      //
      assert(vertex.GetNumFaces() == vertex.GetNumEdges());
  
      int valence = vertex.GetNumFaces();
  
      mask.SetNumVertexWeights(1);
      mask.SetNumEdgeWeights(valence);
      mask.SetNumFaceWeights(valence);
      mask.SetFaceWeightsForFaceCenters(true);
  
      Weight vWeight = (Weight)(valence - 2) / (Weight)valence;
      Weight fWeight = 1.0f / (Weight)(valence * valence);
      Weight eWeight = fWeight;
  
      mask.VertexWeight(0) = vWeight;
      for (int i = 0; i < valence; ++i) {
          mask.EdgeWeight(i) = eWeight;
          mask.FaceWeight(i) = fWeight;
      }
  }
  
  //
  //  Limit masks for position:
  //
  template <>
  template <typename VERTEX, typename MASK>
  inline void
  Scheme<SCHEME_CATMARK>::assignCornerLimitMask(VERTEX const& /* vertex */, MASK& posMask) const {
  
      posMask.SetNumVertexWeights(1);
      posMask.SetNumEdgeWeights(0);
      posMask.SetNumFaceWeights(0);
      posMask.SetFaceWeightsForFaceCenters(false);
  
      posMask.VertexWeight(0) = 1.0f;
  }
  
  template <>
  template <typename VERTEX, typename MASK>
  inline void
  Scheme<SCHEME_CATMARK>::assignCreaseLimitMask(VERTEX const& vertex, MASK& posMask,
                                                int const creaseEnds[2]) const {
  
      typedef typename MASK::Weight Weight;
  
      int valence = vertex.GetNumEdges();
  
      posMask.SetNumVertexWeights(1);
      posMask.SetNumEdgeWeights(valence);
      posMask.SetNumFaceWeights(0);
      posMask.SetFaceWeightsForFaceCenters(false);
  
      Weight vWeight = 2.0f / 3.0f;
      Weight eWeight = 1.0f / 6.0f;
  
      posMask.VertexWeight(0) = vWeight;
      for (int i = 0; i < valence; ++i) {
          posMask.EdgeWeight(i) = 0.0f;
      }
      posMask.EdgeWeight(creaseEnds[0]) = eWeight;
      posMask.EdgeWeight(creaseEnds[1]) = eWeight;
  }
  
  template <>
  template <typename VERTEX, typename MASK>
  inline void
  Scheme<SCHEME_CATMARK>::assignSmoothLimitMask(VERTEX const& vertex, MASK& posMask) const {
  
      typedef typename MASK::Weight Weight;
  
      int valence = vertex.GetNumFaces();
      if (valence == 2) {
          assignCornerLimitMask(vertex, posMask);
          return;
      }
  
      posMask.SetNumVertexWeights(1);
      posMask.SetNumEdgeWeights(valence);
      posMask.SetNumFaceWeights(valence);
      posMask.SetFaceWeightsForFaceCenters(false);
  
      //  Specialize for the regular case:
      if (valence == 4) {
          Weight fWeight = 1.0f / 36.0f;
          Weight eWeight = 1.0f /  9.0f;
          Weight vWeight = 4.0f /  9.0f;
  
          posMask.VertexWeight(0) = vWeight;
  
          posMask.EdgeWeight(0) = eWeight;
          posMask.EdgeWeight(1) = eWeight;
          posMask.EdgeWeight(2) = eWeight;
          posMask.EdgeWeight(3) = eWeight;
  
          posMask.FaceWeight(0) = fWeight;
          posMask.FaceWeight(1) = fWeight;
          posMask.FaceWeight(2) = fWeight;
          posMask.FaceWeight(3) = fWeight;
      } else {
          Weight fWeight = 1.0f / (Weight)(valence * (valence + 5.0f));
          Weight eWeight = 4.0f * fWeight;
          Weight vWeight = (Weight)(1.0f - valence * (eWeight + fWeight));
  
          posMask.VertexWeight(0) = vWeight;
          for (int i = 0; i < valence; ++i) {
              posMask.EdgeWeight(i) = eWeight;
              posMask.FaceWeight(i) = fWeight;
          }
      }
  }
  
  //
  //  Limit masks for tangents -- these are stubs for now, or have a temporary
  //  implementation
  //
  template <>
  template <typename VERTEX, typename MASK>
  inline void
  Scheme<SCHEME_CATMARK>::assignCornerLimitTangentMasks(VERTEX const& vertex,
          MASK& tan1Mask, MASK& tan2Mask) const {
  
      int valence = vertex.GetNumEdges();
  
      tan1Mask.SetNumVertexWeights(1);
      tan1Mask.SetNumEdgeWeights(valence);
      tan1Mask.SetNumFaceWeights(0);
      tan1Mask.SetFaceWeightsForFaceCenters(false);
  
      tan2Mask.SetNumVertexWeights(1);
      tan2Mask.SetNumEdgeWeights(valence);
      tan2Mask.SetNumFaceWeights(0);
      tan2Mask.SetFaceWeightsForFaceCenters(false);
  
      //  Should be at least 2 edges -- be sure to clear weights for any more:
      tan1Mask.VertexWeight(0) = -1.0f;
      tan1Mask.EdgeWeight(0)   =  1.0f;
      tan1Mask.EdgeWeight(1)   =  0.0f;
  
      tan2Mask.VertexWeight(0) = -1.0f;
      tan2Mask.EdgeWeight(0)   =  0.0f;
      tan2Mask.EdgeWeight(1)   =  1.0f;
  
      for (int i = 2; i < valence; ++i) {
          tan1Mask.EdgeWeight(i) = 0.0f;
          tan2Mask.EdgeWeight(i) = 0.0f;
      }
  }
  
  template <>
  template <typename VERTEX, typename MASK>
  inline void
  Scheme<SCHEME_CATMARK>::assignCreaseLimitTangentMasks(VERTEX const& vertex,
          MASK& tan1Mask, MASK& tan2Mask, int const creaseEnds[2]) const {
  
      typedef typename MASK::Weight Weight;
  
      //
      //  First, the tangent along the crease:
      //      The first crease edge is considered the "leading" edge of the span
      //  of surface for which we are evaluating tangents and the second edge the
      //  "trailing edge".  By convention, the tangent along the crease is oriented
      //  in the direction of the leading edge.
      //
      int numEdges = vertex.GetNumEdges();
      int numFaces = vertex.GetNumFaces();
  
      tan1Mask.SetNumVertexWeights(1);
      tan1Mask.SetNumEdgeWeights(numEdges);
      tan1Mask.SetNumFaceWeights(numFaces);
      tan1Mask.SetFaceWeightsForFaceCenters(false);
  
      tan1Mask.VertexWeight(0) = 0.0f;
      for (int i = 0; i < numEdges; ++i) {
          tan1Mask.EdgeWeight(i) = 0.0f;
      }
      for (int i = 0; i < numFaces; ++i) {
          tan1Mask.FaceWeight(i) = 0.0f;
      }
  
      tan1Mask.EdgeWeight(creaseEnds[0]) =  0.5f;
      tan1Mask.EdgeWeight(creaseEnds[1]) = -0.5f;
  
      //
      //  Second, the tangent across the interior faces:
      //      Note this is ambigous for an interior vertex.  We currently return
      //  the tangent for the surface in the counter-clockwise span between the
      //  leading and trailing edges that form the crease.  Given the expected
      //  computation of a surface normal as Tan1 X Tan2, this tangent should be
      //  oriented "inward" from the crease/boundary -- across the surface rather
      //  than outward and away from it.
      //
      tan2Mask.SetNumVertexWeights(1);
      tan2Mask.SetNumEdgeWeights(numEdges);
      tan2Mask.SetNumFaceWeights(numFaces);
      tan2Mask.SetFaceWeightsForFaceCenters(false);
  
      //  Prepend weights of 0 preceding the crease:
      for (int i = 0; i < creaseEnds[0]; ++i) {
          tan2Mask.EdgeWeight(i) = 0.0f;
          tan2Mask.FaceWeight(i) = 0.0f;
      }
  
      //  Assign weights to crease edge and interior points:
      int interiorEdgeCount = creaseEnds[1] - creaseEnds[0] - 1;
      if (interiorEdgeCount == 1) {
          //  The regular case -- uniform B-spline cross-tangent:
  
          tan2Mask.VertexWeight(0) = -4.0f / 6.0f;
  
          tan2Mask.EdgeWeight(creaseEnds[0])     = -1.0f / 6.0f;
          tan2Mask.EdgeWeight(creaseEnds[0] + 1) =  4.0f / 6.0f;
          tan2Mask.EdgeWeight(creaseEnds[1])     = -1.0f / 6.0f;
  
          tan2Mask.FaceWeight(creaseEnds[0])     = 1.0f / 6.0f;
          tan2Mask.FaceWeight(creaseEnds[0] + 1) = 1.0f / 6.0f;
      } else if (interiorEdgeCount > 1) {
          //  The irregular case -- formulae from Biermann et al:
  
          double k     = (double) (interiorEdgeCount + 1);
          double theta = M_PI / k;
  
          double cosTheta = std::cos(theta);
          double sinTheta = std::sin(theta);
  
          //  Loop/Schaefer use a different divisor here (3*k + cos(theta)):
          double commonDenom = 1.0f / (k * (3.0f + cosTheta));
          double R = (cosTheta + 1.0f) / sinTheta;
  
          double vertexWeight = 4.0f * R * (cosTheta - 1.0f);
          double creaseWeight = -R * (1.0f + 2.0f * cosTheta);
  
          tan2Mask.VertexWeight(0) = (Weight) (vertexWeight * commonDenom);
  
          tan2Mask.EdgeWeight(creaseEnds[0]) = (Weight) (creaseWeight * commonDenom);
          tan2Mask.EdgeWeight(creaseEnds[1]) = (Weight) (creaseWeight * commonDenom);
  
          tan2Mask.FaceWeight(creaseEnds[0]) = (Weight) (sinTheta * commonDenom);
  
          double sinThetaI      = 0.0f;
          double sinThetaIplus1 = sinTheta;
          for (int i = 1; i < k; ++i) {
              sinThetaI      = sinThetaIplus1;
              sinThetaIplus1 = std::sin((i+1)*theta);
  
              tan2Mask.EdgeWeight(creaseEnds[0] + i) = (Weight) ((4.0f * sinThetaI) * commonDenom);
              tan2Mask.FaceWeight(creaseEnds[0] + i) = (Weight) ((sinThetaI + sinThetaIplus1) * commonDenom);
          }
      } else {
          //  Special case for a single face -- simple average of boundary edges:
  
          tan2Mask.VertexWeight(0) = -6.0f;
  
          tan2Mask.EdgeWeight(creaseEnds[0]) = 3.0f;
          tan2Mask.EdgeWeight(creaseEnds[1]) = 3.0f;
  
          tan2Mask.FaceWeight(creaseEnds[0]) = 0.0f;
      }
  
      //  Append weights of 0 following the crease:
      for (int i = creaseEnds[1]; i < numFaces; ++i) {
          tan2Mask.FaceWeight(i) = 0.0f;
      }
      for (int i = creaseEnds[1] + 1; i < numEdges; ++i) {
          tan2Mask.EdgeWeight(i) = 0.0f;
      }
  }
  
  template <>
  template <typename VERTEX, typename MASK>
  inline void
  Scheme<SCHEME_CATMARK>::assignSmoothLimitTangentMasks(VERTEX const& vertex,
          MASK& tan1Mask, MASK& tan2Mask) const {
  
      typedef typename MASK::Weight Weight;
  
      int valence = vertex.GetNumFaces();
      if (valence == 2) {
          assignCornerLimitTangentMasks(vertex, tan1Mask, tan2Mask);
          return;
      }
  
      //  Compute tan1 initially -- tan2 is simply a rotation:
      tan1Mask.SetNumVertexWeights(1);
      tan1Mask.SetNumEdgeWeights(valence);
      tan1Mask.SetNumFaceWeights(valence);
      tan1Mask.SetFaceWeightsForFaceCenters(false);
  
      tan1Mask.VertexWeight(0) = 0.0f;
  
      if (valence == 4) {
          tan1Mask.EdgeWeight(0) =  4.0f;
          tan1Mask.EdgeWeight(1) =  0.0f;
          tan1Mask.EdgeWeight(2) = -4.0f;
          tan1Mask.EdgeWeight(3) =  0.0f;
  
          tan1Mask.FaceWeight(0) =  1.0f;
          tan1Mask.FaceWeight(1) = -1.0f;
          tan1Mask.FaceWeight(2) = -1.0f;
          tan1Mask.FaceWeight(3) =  1.0f;
      } else {
          double theta = 2.0f * M_PI / (double)valence;
  
          double cosTheta     = std::cos(theta);
          double cosHalfTheta = std::cos(theta * 0.5f);
  
          double lambda = (5.0f / 16.0f) + (1.0f / 16.0f) *
                  (cosTheta + cosHalfTheta * std::sqrt(2.0f * (9.0f + cosTheta)));
  
          double edgeWeightScale = 4.0f;
          double faceWeightScale = 1.0f / (4.0f * lambda - 1.0f);
  
          for (int i = 0; i < valence; ++i) {
              double cosThetaI      = std::cos(  i  * theta);
              double cosThetaIplus1 = std::cos((i+1)* theta);
  
              tan1Mask.EdgeWeight(i) = (Weight) (edgeWeightScale * cosThetaI);
              tan1Mask.FaceWeight(i) = (Weight) (faceWeightScale * (cosThetaI + cosThetaIplus1));
          }
      }
  
      //  Now rotate/copy tan1 weights to tan2:
      tan2Mask.SetNumVertexWeights(1);
      tan2Mask.SetNumEdgeWeights(valence);
      tan2Mask.SetNumFaceWeights(valence);
      tan2Mask.SetFaceWeightsForFaceCenters(false);
  
      tan2Mask.VertexWeight(0) = 0.0f;
      if (valence == 4) {
          tan2Mask.EdgeWeight(0) =  0.0f;
          tan2Mask.EdgeWeight(1) =  4.0f;
          tan2Mask.EdgeWeight(2) =  0.0f;
          tan2Mask.EdgeWeight(3) = -4.0f;
  
          tan2Mask.FaceWeight(0) =  1.0f;
          tan2Mask.FaceWeight(1) =  1.0f;
          tan2Mask.FaceWeight(2) = -1.0f;
          tan2Mask.FaceWeight(3) = -1.0f;
      } else {
          tan2Mask.EdgeWeight(0) = tan1Mask.EdgeWeight(valence-1);
          tan2Mask.FaceWeight(0) = tan1Mask.FaceWeight(valence-1);
          for (int i = 1; i < valence; ++i) {
              tan2Mask.EdgeWeight(i) = tan1Mask.EdgeWeight(i-1);
              tan2Mask.FaceWeight(i) = tan1Mask.FaceWeight(i-1);
          }
      }
  }
  
  } // end namespace sdc
  
  } // end namespace OPENSUBDIV_VERSION
  using namespace OPENSUBDIV_VERSION;
  } // end namespace OpenSubdiv
  
  #endif /* OPENSUBDIV3_SDC_CATMARK_SCHEME_H */
  