subdivision.cpp

// Copyright 2024 The Manifold Authors.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//      http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include <unordered_map>

#include "impl.h"
#include "parallel.h"

template <>
struct std::hash<manifold::ivec4> {
  size_t operator()(const manifold::ivec4& p) const {
    return std::hash<int>()(p.x) ^ std::hash<int>()(p.y) ^
           std::hash<int>()(p.z) ^ std::hash<int>()(p.w);
  }
};

namespace {
using namespace manifold;

class Partition {
 public:
  // The cached partitions don't have idx - it's added to the copy returned
  // from GetPartition that contains the mapping of the input divisions into the
  // sorted divisions that are uniquely cached.
  ivec4 idx;
  ivec4 sortedDivisions;
  Vec<vec4> vertBary;
  Vec<ivec3> triVert;

  int InteriorOffset() const {
    return sortedDivisions[0] + sortedDivisions[1] + sortedDivisions[2] +
           sortedDivisions[3];
  }

  int NumInterior() const { return vertBary.size() - InteriorOffset(); }

  static Partition GetPartition(ivec4 divisions) {
    if (divisions[0] == 0) return Partition();  // skip wrong side of quad

    ivec4 sortedDiv = divisions;
    ivec4 triIdx = {0, 1, 2, 3};
    if (divisions[3] == 0) {  // triangle
      if (sortedDiv[2] > sortedDiv[1]) {
        std::swap(sortedDiv[2], sortedDiv[1]);
        std::swap(triIdx[2], triIdx[1]);
      }
      if (sortedDiv[1] > sortedDiv[0]) {
        std::swap(sortedDiv[1], sortedDiv[0]);
        std::swap(triIdx[1], triIdx[0]);
        if (sortedDiv[2] > sortedDiv[1]) {
          std::swap(sortedDiv[2], sortedDiv[1]);
          std::swap(triIdx[2], triIdx[1]);
        }
      }
    } else {  // quad
      int minIdx = 0;
      int min = divisions[minIdx];
      int next = divisions[1];
      for (const int i : {1, 2, 3}) {
        const int n = divisions[(i + 1) % 4];
        if (divisions[i] < min || (divisions[i] == min && n < next)) {
          minIdx = i;
          min = divisions[i];
          next = n;
        }
      }
      // Backwards (mirrored) quads get a separate cache key for now for
      // simplicity, so there is no reversal necessary for quads when
      // re-indexing.
      ivec4 tmp = sortedDiv;
      for (const int i : {0, 1, 2, 3}) {
        triIdx[i] = (i + minIdx) % 4;
        sortedDiv[i] = tmp[triIdx[i]];
      }
    }

    Partition partition = GetCachedPartition(sortedDiv);
    partition.idx = triIdx;

    return partition;
  }

  Vec<ivec3> Reindex(ivec4 triVerts, ivec4 edgeOffsets, bvec4 edgeFwd,
                     int interiorOffset) const {
    Vec<int> newVerts;
    newVerts.reserve(vertBary.size());
    ivec4 triIdx = idx;
    ivec4 outTri = {0, 1, 2, 3};
    if (triVerts[3] < 0 && idx[1] != Next3(idx[0])) {
      triIdx = {idx[2], idx[0], idx[1], idx[3]};
      edgeFwd = !edgeFwd;
      std::swap(outTri[0], outTri[1]);
    }
    for (const int i : {0, 1, 2, 3}) {
      if (triVerts[triIdx[i]] >= 0) newVerts.push_back(triVerts[triIdx[i]]);
    }
    for (const int i : {0, 1, 2, 3}) {
      const int n = sortedDivisions[i] - 1;
      int offset = edgeOffsets[idx[i]] + (edgeFwd[idx[i]] ? 0 : n - 1);
      for (int j = 0; j < n; ++j) {
        newVerts.push_back(offset);
        offset += edgeFwd[idx[i]] ? 1 : -1;
      }
    }
    const int offset = interiorOffset - newVerts.size();
    size_t old = newVerts.size();
    newVerts.resize_nofill(vertBary.size());
    std::iota(newVerts.begin() + old, newVerts.end(), old + offset);

    const int numTri = triVert.size();
    Vec<ivec3> newTriVert(numTri);
    for_each_n(autoPolicy(numTri), countAt(0), numTri,
               [&newTriVert, &outTri, &newVerts, this](const int tri) {
                 for (const int j : {0, 1, 2}) {
                   newTriVert[tri][outTri[j]] = newVerts[triVert[tri][j]];
                 }
               });
    return newTriVert;
  }

 private:
  static inline auto cacheLock = std::mutex();
  static inline auto cache =
      std::unordered_map<ivec4, std::unique_ptr<Partition>>();

  // This triangulation is purely topological - it depends only on the number of
  // divisions of the three sides of the triangle. This allows them to be cached
  // and reused for similar triangles. The shape of the final surface is defined
  // by the tangents and the barycentric coordinates of the new verts. For
  // triangles, the input must be sorted: n[0] >= n[1] >= n[2] > 0.
  static Partition GetCachedPartition(ivec4 n) {
    {
      auto lockGuard = std::lock_guard<std::mutex>(cacheLock);
      auto cached = cache.find(n);
      if (cached != cache.end()) {
        return *cached->second;
      }
    }
    Partition partition;
    partition.sortedDivisions = n;
    if (n[3] > 0) {  // quad
      partition.vertBary.push_back({1, 0, 0, 0});
      partition.vertBary.push_back({0, 1, 0, 0});
      partition.vertBary.push_back({0, 0, 1, 0});
      partition.vertBary.push_back({0, 0, 0, 1});
      ivec4 edgeOffsets;
      edgeOffsets[0] = 4;
      for (const int i : {0, 1, 2, 3}) {
        if (i > 0) {
          edgeOffsets[i] = edgeOffsets[i - 1] + n[i - 1] - 1;
        }
        const vec4 nextBary = partition.vertBary[(i + 1) % 4];
        for (int j = 1; j < n[i]; ++j) {
          partition.vertBary.push_back(
              la::lerp(partition.vertBary[i], nextBary, (double)j / n[i]));
        }
      }
      PartitionQuad(partition.triVert, partition.vertBary, {0, 1, 2, 3},
                    edgeOffsets, n - 1, {true, true, true, true});
    } else {  // tri
      partition.vertBary.push_back({1, 0, 0, 0});
      partition.vertBary.push_back({0, 1, 0, 0});
      partition.vertBary.push_back({0, 0, 1, 0});
      for (const int i : {0, 1, 2}) {
        const vec4 nextBary = partition.vertBary[(i + 1) % 3];
        for (int j = 1; j < n[i]; ++j) {
          partition.vertBary.push_back(
              la::lerp(partition.vertBary[i], nextBary, (double)j / n[i]));
        }
      }
      const ivec3 edgeOffsets = {3, 3 + n[0] - 1, 3 + n[0] - 1 + n[1] - 1};

      const double f = n[2] * n[2] + n[0] * n[0];
      if (n[1] == 1) {
        if (n[0] == 1) {
          partition.triVert.push_back({0, 1, 2});
        } else {
          PartitionFan(partition.triVert, {0, 1, 2}, n[0] - 1, edgeOffsets[0]);
        }
      } else if (n[1] * n[1] > f - std::sqrt(2.0) * n[0] * n[2]) {  // acute-ish
        partition.triVert.push_back({edgeOffsets[1] - 1, 1, edgeOffsets[1]});
        PartitionQuad(partition.triVert, partition.vertBary,
                      {edgeOffsets[1] - 1, edgeOffsets[1], 2, 0},
                      {-1, edgeOffsets[1] + 1, edgeOffsets[2], edgeOffsets[0]},
                      {0, n[1] - 2, n[2] - 1, n[0] - 2},
                      {true, true, true, true});
      } else {  // obtuse -> spit into two acute
        // portion of n[0] under n[2]
        const int ns =
            std::min(n[0] - 2, (int)std::round((f - n[1] * n[1]) / (2 * n[0])));
        // height from n[0]: nh <= n[2]
        const int nh =
            std::max(1., std::round(std::sqrt(n[2] * n[2] - ns * ns)));

        const int hOffset = partition.vertBary.size();
        const vec4 middleBary = partition.vertBary[edgeOffsets[0] + ns - 1];
        for (int j = 1; j < nh; ++j) {
          partition.vertBary.push_back(
              la::lerp(partition.vertBary[2], middleBary, (double)j / nh));
        }

        partition.triVert.push_back({edgeOffsets[1] - 1, 1, edgeOffsets[1]});
        PartitionQuad(
            partition.triVert, partition.vertBary,
            {edgeOffsets[1] - 1, edgeOffsets[1], 2, edgeOffsets[0] + ns - 1},
            {-1, edgeOffsets[1] + 1, hOffset, edgeOffsets[0] + ns},
            {0, n[1] - 2, nh - 1, n[0] - ns - 2}, {true, true, true, true});

        if (n[2] == 1) {
          PartitionFan(partition.triVert, {0, edgeOffsets[0] + ns - 1, 2},
                       ns - 1, edgeOffsets[0]);
        } else {
          if (ns == 1) {
            partition.triVert.push_back({hOffset, 2, edgeOffsets[2]});
            PartitionQuad(partition.triVert, partition.vertBary,
                          {hOffset, edgeOffsets[2], 0, edgeOffsets[0]},
                          {-1, edgeOffsets[2] + 1, -1, hOffset + nh - 2},
                          {0, n[2] - 2, ns - 1, nh - 2},
                          {true, true, true, false});
          } else {
            partition.triVert.push_back({hOffset - 1, 0, edgeOffsets[0]});
            PartitionQuad(
                partition.triVert, partition.vertBary,
                {hOffset - 1, edgeOffsets[0], edgeOffsets[0] + ns - 1, 2},
                {-1, edgeOffsets[0] + 1, hOffset + nh - 2, edgeOffsets[2]},
                {0, ns - 2, nh - 1, n[2] - 2}, {true, true, false, true});
          }
        }
      }
    }

    auto lockGuard = std::lock_guard<std::mutex>(cacheLock);
    cache.insert({n, std::make_unique<Partition>(partition)});
    return partition;
  }

  // Side 0 has added edges while sides 1 and 2 do not. Fan spreads from vert 2.
  static void PartitionFan(Vec<ivec3>& triVert, ivec3 cornerVerts, int added,
                           int edgeOffset) {
    int last = cornerVerts[0];
    for (int i = 0; i < added; ++i) {
      const int next = edgeOffset + i;
      triVert.push_back({last, next, cornerVerts[2]});
      last = next;
    }
    triVert.push_back({last, cornerVerts[1], cornerVerts[2]});
  }

  // Partitions are parallel to the first edge unless two consecutive edgeAdded
  // are zero, in which case a terminal triangulation is performed.
  static void PartitionQuad(Vec<ivec3>& triVert, Vec<vec4>& vertBary,
                            ivec4 cornerVerts, ivec4 edgeOffsets,
                            ivec4 edgeAdded, bvec4 edgeFwd) {
    auto GetEdgeVert = [&](int edge, int idx) {
      return edgeOffsets[edge] + (edgeFwd[edge] ? 1 : -1) * idx;
    };

    DEBUG_ASSERT(la::all(la::gequal(edgeAdded, ivec4(0))), logicErr,
                 "negative divisions!");

    int corner = -1;
    int last = 3;
    int maxEdge = -1;
    for (const int i : {0, 1, 2, 3}) {
      if (corner == -1 && edgeAdded[i] == 0 && edgeAdded[last] == 0) {
        corner = i;
      }
      if (edgeAdded[i] > 0) {
        maxEdge = maxEdge == -1 ? i : -2;
      }
      last = i;
    }
    if (corner >= 0) {  // terminate
      if (maxEdge >= 0) {
        ivec4 edge = (ivec4(0, 1, 2, 3) + maxEdge) % 4;
        const int middle = edgeAdded[maxEdge] / 2;
        triVert.push_back({cornerVerts[edge[2]], cornerVerts[edge[3]],
                           GetEdgeVert(maxEdge, middle)});
        int last = cornerVerts[edge[0]];
        for (int i = 0; i <= middle; ++i) {
          const int next = GetEdgeVert(maxEdge, i);
          triVert.push_back({cornerVerts[edge[3]], last, next});
          last = next;
        }
        last = cornerVerts[edge[1]];
        for (int i = edgeAdded[maxEdge] - 1; i >= middle; --i) {
          const int next = GetEdgeVert(maxEdge, i);
          triVert.push_back({cornerVerts[edge[2]], next, last});
          last = next;
        }
      } else {
        int sideVert = cornerVerts[0];  // initial value is unused
        for (const int j : {1, 2}) {
          const int side = (corner + j) % 4;
          if (j == 2 && edgeAdded[side] > 0) {
            triVert.push_back(
                {cornerVerts[side], GetEdgeVert(side, 0), sideVert});
          } else {
            sideVert = cornerVerts[side];
          }
          for (int i = 0; i < edgeAdded[side]; ++i) {
            const int nextVert = GetEdgeVert(side, i);
            triVert.push_back({cornerVerts[corner], sideVert, nextVert});
            sideVert = nextVert;
          }
          if (j == 2 || edgeAdded[side] == 0) {
            triVert.push_back({cornerVerts[corner], sideVert,
                               cornerVerts[(corner + j + 1) % 4]});
          }
        }
      }
      return;
    }
    // recursively partition
    const int partitions = 1 + std::min(edgeAdded[1], edgeAdded[3]);
    ivec4 newCornerVerts = {cornerVerts[1], -1, -1, cornerVerts[0]};
    ivec4 newEdgeOffsets = {edgeOffsets[1], -1,
                            GetEdgeVert(3, edgeAdded[3] + 1), edgeOffsets[0]};
    ivec4 newEdgeAdded = {0, -1, 0, edgeAdded[0]};
    bvec4 newEdgeFwd = {edgeFwd[1], true, edgeFwd[3], edgeFwd[0]};

    for (int i = 1; i < partitions; ++i) {
      const int cornerOffset1 = (edgeAdded[1] * i) / partitions;
      const int cornerOffset3 =
          edgeAdded[3] - 1 - (edgeAdded[3] * i) / partitions;
      const int nextOffset1 = GetEdgeVert(1, cornerOffset1 + 1);
      const int nextOffset3 = GetEdgeVert(3, cornerOffset3 + 1);
      const int added = std::round(la::lerp(
          (double)edgeAdded[0], (double)edgeAdded[2], (double)i / partitions));

      newCornerVerts[1] = GetEdgeVert(1, cornerOffset1);
      newCornerVerts[2] = GetEdgeVert(3, cornerOffset3);
      newEdgeAdded[0] = std::abs(nextOffset1 - newEdgeOffsets[0]) - 1;
      newEdgeAdded[1] = added;
      newEdgeAdded[2] = std::abs(nextOffset3 - newEdgeOffsets[2]) - 1;
      newEdgeOffsets[1] = vertBary.size();
      newEdgeOffsets[2] = nextOffset3;

      for (int j = 0; j < added; ++j) {
        vertBary.push_back(la::lerp(vertBary[newCornerVerts[1]],
                                    vertBary[newCornerVerts[2]],
                                    (j + 1.0) / (added + 1.0)));
      }

      PartitionQuad(triVert, vertBary, newCornerVerts, newEdgeOffsets,
                    newEdgeAdded, newEdgeFwd);

      newCornerVerts[0] = newCornerVerts[1];
      newCornerVerts[3] = newCornerVerts[2];
      newEdgeAdded[3] = newEdgeAdded[1];
      newEdgeOffsets[0] = nextOffset1;
      newEdgeOffsets[3] = newEdgeOffsets[1] + newEdgeAdded[1] - 1;
      newEdgeFwd[3] = false;
    }

    newCornerVerts[1] = cornerVerts[2];
    newCornerVerts[2] = cornerVerts[3];
    newEdgeOffsets[1] = edgeOffsets[2];
    newEdgeAdded[0] =
        edgeAdded[1] - std::abs(newEdgeOffsets[0] - edgeOffsets[1]);
    newEdgeAdded[1] = edgeAdded[2];
    newEdgeAdded[2] = std::abs(newEdgeOffsets[2] - edgeOffsets[3]) - 1;
    newEdgeOffsets[2] = edgeOffsets[3];
    newEdgeFwd[1] = edgeFwd[2];

    PartitionQuad(triVert, vertBary, newCornerVerts, newEdgeOffsets,
                  newEdgeAdded, newEdgeFwd);
  }
};
}  // namespace

namespace manifold {

/**
 * Returns the tri side index (0-2) connected to the other side of this quad if
 * this tri is part of a quad, or -1 otherwise.
 */
int Manifold::Impl::GetNeighbor(int tri) const {
  int neighbor = -1;
  for (const int i : {0, 1, 2}) {
    if (IsMarkedInsideQuad(3 * tri + i)) {
      neighbor = neighbor == -1 ? i : -2;
    }
  }
  return neighbor;
}

/**
 * For the given triangle index, returns either the three halfedge indices of
 * that triangle and halfedges[3] = -1, or if the triangle is part of a quad, it
 * returns those four indices. If the triangle is part of a quad and is not the
 * lower of the two triangle indices, it returns all -1s.
 */
ivec4 Manifold::Impl::GetHalfedges(int tri) const {
  ivec4 halfedges(-1);
  for (const int i : {0, 1, 2}) {
    halfedges[i] = 3 * tri + i;
  }
  const int neighbor = GetNeighbor(tri);
  if (neighbor >= 0) {  // quad
    const int pair = halfedge_[3 * tri + neighbor].pairedHalfedge;
    if (pair / 3 < tri) {
      return ivec4(-1);  // only process lower tri index
    }
    // The order here matters to keep small quads split the way they started, or
    // else it can create a 4-manifold edge.
    halfedges[2] = NextHalfedge(halfedges[neighbor]);
    halfedges[3] = NextHalfedge(halfedges[2]);
    halfedges[0] = NextHalfedge(pair);
    halfedges[1] = NextHalfedge(halfedges[0]);
  }
  return halfedges;
}

/**
 * Returns the BaryIndices, which gives the tri and indices (0-3), such that
 * GetHalfedges(val.tri)[val.start4] points back to this halfedge, and val.end4
 * will point to the next one. This function handles this for both triangles and
 * quads. Returns {-1, -1, -1} if the edge is the interior of a quad.
 */
Manifold::Impl::BaryIndices Manifold::Impl::GetIndices(int halfedge) const {
  int tri = halfedge / 3;
  int idx = halfedge % 3;
  const int neighbor = GetNeighbor(tri);
  if (idx == neighbor) {
    return {-1, -1, -1};
  }

  if (neighbor < 0) {  // tri
    return {tri, idx, Next3(idx)};
  } else {  // quad
    const int pair = halfedge_[3 * tri + neighbor].pairedHalfedge;
    if (pair / 3 < tri) {
      tri = pair / 3;
      idx = Next3(neighbor) == idx ? 0 : 1;
    } else {
      idx = Next3(neighbor) == idx ? 2 : 3;
    }
    return {tri, idx, (idx + 1) % 4};
  }
}

/**
 * Retained verts are part of several triangles, and it doesn't matter which one
 * the vertBary refers to. Here, whichever is last will win and it's done on the
 * CPU for simplicity for now. Using AtomicCAS on .tri should work for a GPU
 * version if desired.
 */
void Manifold::Impl::FillRetainedVerts(Vec<Barycentric>& vertBary) const {
  const int numTri = halfedge_.size() / 3;
  for (int tri = 0; tri < numTri; ++tri) {
    for (const int i : {0, 1, 2}) {
      const BaryIndices indices = GetIndices(3 * tri + i);
      if (indices.start4 < 0) continue;  // skip quad interiors
      vec4 uvw(0.0);
      uvw[indices.start4] = 1;
      vertBary[halfedge_[3 * tri + i].startVert] = {indices.tri, uvw};
    }
  }
}

/**
 * Split each edge into n pieces as defined by calling the edgeDivisions
 * function, and sub-triangulate each triangle accordingly. This function
 * doesn't run Finish(), as that is expensive and it'll need to be run after
 * the new vertices have moved, which is a likely scenario after refinement
 * (smoothing).
 */
Vec<Barycentric> Manifold::Impl::Subdivide(
    std::function<int(vec3, vec4, vec4)> edgeDivisions, bool keepInterior) {
  Vec<TmpEdge> edges = CreateTmpEdges(halfedge_);
  const int numVert = NumVert();
  const int numEdge = edges.size();
  const int numTri = NumTri();
  Vec<int> half2Edge(2 * numEdge);
  auto policy = autoPolicy(numEdge, 1e4);
  for_each_n(policy, countAt(0), numEdge,
             [&half2Edge, &edges, this](const int edge) {
               const int idx = edges[edge].halfedgeIdx;
               half2Edge[idx] = edge;
               half2Edge[halfedge_[idx].pairedHalfedge] = edge;
             });

  Vec<ivec4> faceHalfedges(numTri);
  for_each_n(policy, countAt(0), numTri, [&faceHalfedges, this](const int tri) {
    faceHalfedges[tri] = GetHalfedges(tri);
  });

  Vec<int> edgeAdded(numEdge);
  for_each_n(policy, countAt(0), numEdge,
             [&edgeAdded, &edges, edgeDivisions, this](const int i) {
               const TmpEdge edge = edges[i];
               const int hIdx = edge.halfedgeIdx;
               if (IsMarkedInsideQuad(hIdx)) {
                 edgeAdded[i] = 0;
                 return;
               }
               const vec3 vec = vertPos_[edge.first] - vertPos_[edge.second];
               const vec4 tangent0 = halfedgeTangent_.empty()
                                         ? vec4(0.0)
                                         : halfedgeTangent_[hIdx];
               const vec4 tangent1 =
                   halfedgeTangent_.empty()
                       ? vec4(0.0)
                       : halfedgeTangent_[halfedge_[hIdx].pairedHalfedge];
               edgeAdded[i] = edgeDivisions(vec, tangent0, tangent1);
             });

  if (keepInterior) {
    // Triangles where the greatest number of divisions exceeds the sum of the
    // other two sides will be triangulated as a strip, since if the sub-edges
    // were all equal length it would be degenerate. This leads to poor results
    // with RefineToTolerance, so we avoid this case by adding some extra
    // divisions to the short sides so that the triangulation has some thickness
    // and creates more interior facets.
    Vec<int> tmp(numEdge);
    for_each_n(
        policy, countAt(0), numEdge,
        [&tmp, &edgeAdded, &edges, &half2Edge, this](const int i) {
          tmp[i] = edgeAdded[i];
          const TmpEdge edge = edges[i];
          int hIdx = edge.halfedgeIdx;
          if (IsMarkedInsideQuad(hIdx)) return;

          const int thisAdded = tmp[i];
          auto Added = [&edgeAdded, &half2Edge, thisAdded, this](int hIdx) {
            int longest = 0;
            int total = 0;
            for (int _ : {0, 1, 2}) {
              const int added = edgeAdded[half2Edge[hIdx]];
              longest = la::max(longest, added);
              total += added;
              hIdx = NextHalfedge(hIdx);
              if (IsMarkedInsideQuad(hIdx)) {
                // No extra on quads
                longest = 0;
                total = 1;
                break;
              }
            }
            const int minExtra = longest * 0.2 + 1;
            const int extra = 2 * longest + minExtra - total;
            return extra > 0 ? (extra * (longest - thisAdded)) / longest : 0;
          };

          tmp[i] += la::max(Added(hIdx), Added(halfedge_[hIdx].pairedHalfedge));
        });
    edgeAdded.swap(tmp);
  }

  Vec<int> edgeOffset(numEdge);
  exclusive_scan(edgeAdded.begin(), edgeAdded.end(), edgeOffset.begin(),
                 numVert);

  Vec<Barycentric> vertBary(edgeOffset.back() + edgeAdded.back());
  const int totalEdgeAdded = vertBary.size() - numVert;
  FillRetainedVerts(vertBary);
  for_each_n(policy, countAt(0), numEdge,
             [&vertBary, &edges, &edgeAdded, &edgeOffset, this](const int i) {
               const int n = edgeAdded[i];
               const int offset = edgeOffset[i];

               const BaryIndices indices = GetIndices(edges[i].halfedgeIdx);
               if (indices.tri < 0) {
                 return;  // inside quad
               }
               const double frac = 1.0 / (n + 1);

               for (int i = 0; i < n; ++i) {
                 vec4 uvw(0.0);
                 uvw[indices.end4] = (i + 1) * frac;
                 uvw[indices.start4] = 1 - uvw[indices.end4];
                 vertBary[offset + i].uvw = uvw;
                 vertBary[offset + i].tri = indices.tri;
               }
             });

  std::vector<Partition> subTris(numTri);
  for_each_n(policy, countAt(0), numTri,
             [&subTris, &half2Edge, &edgeAdded, &faceHalfedges](int tri) {
               const ivec4 halfedges = faceHalfedges[tri];
               ivec4 divisions(0);
               for (const int i : {0, 1, 2, 3}) {
                 if (halfedges[i] >= 0) {
                   divisions[i] = edgeAdded[half2Edge[halfedges[i]]] + 1;
                 }
               }
               subTris[tri] = Partition::GetPartition(divisions);
             });

  Vec<int> triOffset(numTri);
  auto numSubTris =
      TransformIterator(subTris.begin(), [](const Partition& part) {
        return static_cast<int>(part.triVert.size());
      });
  manifold::exclusive_scan(numSubTris, numSubTris + numTri, triOffset.begin(),
                           0);

  Vec<int> interiorOffset(numTri);
  auto numInterior =
      TransformIterator(subTris.begin(), [](const Partition& part) {
        return static_cast<int>(part.NumInterior());
      });
  manifold::exclusive_scan(numInterior, numInterior + numTri,
                           interiorOffset.begin(),
                           static_cast<int>(vertBary.size()));

  Vec<ivec3> triVerts(triOffset.back() + subTris.back().triVert.size());
  vertBary.resize(interiorOffset.back() + subTris.back().NumInterior());
  Vec<TriRef> triRef(triVerts.size());
  for_each_n(
      policy, countAt(0), numTri,
      [this, &triVerts, &triRef, &vertBary, &subTris, &edgeOffset, &half2Edge,
       &triOffset, &interiorOffset, &faceHalfedges](int tri) {
        const ivec4 halfedges = faceHalfedges[tri];
        if (halfedges[0] < 0) return;
        ivec4 tri3;
        ivec4 edgeOffsets;
        bvec4 edgeFwd(false);
        for (const int i : {0, 1, 2, 3}) {
          if (halfedges[i] < 0) {
            tri3[i] = -1;
            continue;
          }
          const Halfedge& halfedge = halfedge_[halfedges[i]];
          tri3[i] = halfedge.startVert;
          edgeOffsets[i] = edgeOffset[half2Edge[halfedges[i]]];
          edgeFwd[i] = halfedge.IsForward();
        }

        Vec<ivec3> newTris = subTris[tri].Reindex(tri3, edgeOffsets, edgeFwd,
                                                  interiorOffset[tri]);
        copy(newTris.begin(), newTris.end(), triVerts.begin() + triOffset[tri]);
        auto start = triRef.begin() + triOffset[tri];
        fill(start, start + newTris.size(), meshRelation_.triRef[tri]);

        const ivec4 idx = subTris[tri].idx;
        const ivec4 vIdx = halfedges[3] >= 0 || idx[1] == Next3(idx[0])
                               ? idx
                               : ivec4(idx[2], idx[0], idx[1], idx[3]);
        ivec4 rIdx;
        for (const int i : {0, 1, 2, 3}) {
          rIdx[vIdx[i]] = i;
        }

        const auto& subBary = subTris[tri].vertBary;
        transform(subBary.begin() + subTris[tri].InteriorOffset(),
                  subBary.end(), vertBary.begin() + interiorOffset[tri],
                  [tri, rIdx](vec4 bary) {
                    return Barycentric({tri,
                                        {bary[rIdx[0]], bary[rIdx[1]],
                                         bary[rIdx[2]], bary[rIdx[3]]}});
                  });
      });
  meshRelation_.triRef = triRef;

  Vec<vec3> newVertPos(vertBary.size());
  for_each_n(policy, countAt(0), vertBary.size(),
             [&newVertPos, &vertBary, &faceHalfedges, this](const int vert) {
               const Barycentric bary = vertBary[vert];
               const ivec4 halfedges = faceHalfedges[bary.tri];
               if (halfedges[3] < 0) {
                 mat3 triPos;
                 for (const int i : {0, 1, 2}) {
                   triPos[i] = vertPos_[halfedge_[halfedges[i]].startVert];
                 }
                 newVertPos[vert] = triPos * vec3(bary.uvw);
               } else {
                 mat3x4 quadPos;
                 for (const int i : {0, 1, 2, 3}) {
                   quadPos[i] = vertPos_[halfedge_[halfedges[i]].startVert];
                 }
                 newVertPos[vert] = quadPos * bary.uvw;
               }
             });
  vertPos_ = newVertPos;

  faceNormal_.clear();

  if (numProp_ > 0) {
    const int numPropVert = NumPropVert();
    const int addedVerts = NumVert() - numVert;
    const int propOffset = numPropVert - numVert;
    // duplicate the prop verts along all new edges even though this is
    // unnecessary for edges that share the same prop verts. The duplicates will
    // be removed by CompactProps.
    Vec<double> prop(numProp_ * (numPropVert + addedVerts + totalEdgeAdded));

    // copy retained prop verts
    copy(properties_.begin(), properties_.end(), prop.begin());

    // copy interior prop verts and forward edge prop verts
    for_each_n(
        policy, countAt(0), addedVerts,
        [&prop, &vertBary, &faceHalfedges, numVert, numPropVert,
         this](const int i) {
          const int vert = numPropVert + i;
          const Barycentric bary = vertBary[numVert + i];
          const ivec4 halfedges = faceHalfedges[bary.tri];
          const int numProp = NumProp();

          for (int p = 0; p < numProp; ++p) {
            if (halfedges[3] < 0) {
              vec3 triProp;
              for (const int i : {0, 1, 2}) {
                triProp[i] =
                    properties_[halfedge_[3 * bary.tri + i].propVert * numProp +
                                p];
              }
              prop[vert * numProp + p] = la::dot(triProp, vec3(bary.uvw));
            } else {
              vec4 quadProp;
              for (const int i : {0, 1, 2, 3}) {
                quadProp[i] =
                    properties_[halfedge_[halfedges[i]].propVert * numProp + p];
              }
              prop[vert * numProp + p] = la::dot(quadProp, bary.uvw);
            }
          }
        });

    // copy backward edge prop verts, some of which will be unreferenced
    // duplicates.
    for_each_n(policy, countAt(0), numEdge,
               [this, &prop, &edges, &edgeAdded, &edgeOffset, propOffset,
                addedVerts](const int i) {
                 const int n = edgeAdded[i];
                 const int offset = edgeOffset[i] + propOffset + addedVerts;
                 const int numProp = NumProp();

                 const double frac = 1.0 / (n + 1);
                 const int halfedgeIdx =
                     halfedge_[edges[i].halfedgeIdx].pairedHalfedge;
                 const int prop0 = halfedge_[halfedgeIdx].propVert;
                 const int prop1 =
                     halfedge_[NextHalfedge(halfedgeIdx)].propVert;
                 for (int i = 0; i < n; ++i) {
                   for (int p = 0; p < numProp; ++p) {
                     prop[(offset + i) * numProp + p] = la::lerp(
                         properties_[prop0 * numProp + p],
                         properties_[prop1 * numProp + p], (i + 1) * frac);
                   }
                 }
               });

    Vec<ivec3> triProp(triVerts.size());
    for_each_n(
        policy, countAt(0), numTri,
        [this, &triProp, &subTris, &edgeOffset, &half2Edge, &triOffset,
         &interiorOffset, &faceHalfedges, propOffset,
         addedVerts](const int tri) {
          const ivec4 halfedges = faceHalfedges[tri];
          if (halfedges[0] < 0) return;

          ivec4 tri3;
          ivec4 edgeOffsets;
          bvec4 edgeFwd(true);
          for (const int i : {0, 1, 2, 3}) {
            if (halfedges[i] < 0) {
              tri3[i] = -1;
              continue;
            }
            const Halfedge& halfedge = halfedge_[halfedges[i]];
            tri3[i] = halfedge.propVert;
            edgeOffsets[i] = edgeOffset[half2Edge[halfedges[i]]];
            if (!halfedge.IsForward()) {
              if (halfedge_[halfedge.pairedHalfedge].propVert !=
                      halfedge_[NextHalfedge(halfedges[i])].propVert ||
                  halfedge_[NextHalfedge(halfedge.pairedHalfedge)].propVert !=
                      halfedge.propVert) {
                // if the edge doesn't match, point to the backward edge
                // propverts.
                edgeOffsets[i] += addedVerts;
              } else {
                edgeFwd[i] = false;
              }
            }
          }

          Vec<ivec3> newTris =
              subTris[tri].Reindex(tri3, edgeOffsets + propOffset, edgeFwd,
                                   interiorOffset[tri] + propOffset);
          copy(newTris.begin(), newTris.end(),
               triProp.begin() + triOffset[tri]);
        });

    properties_ = prop;
    CreateHalfedges(triProp, triVerts);
  } else {
    CreateHalfedges(triVerts);
  }

  return vertBary;
}

}  // namespace manifold