boolean_result.cpp

// Copyright 2021 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 <algorithm>
#include <array>
#include <map>

#include "boolean3.h"
#include "parallel.h"
#include "utils.h"

#if (MANIFOLD_PAR == 1) && __has_include(<tbb/concurrent_map.h>)
#define TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS 1
#include <tbb/concurrent_map.h>
#include <tbb/parallel_for.h>

template <typename K, typename V>
using concurrent_map = tbb::concurrent_map<K, V>;
#else
template <typename K, typename V>
// not really concurrent when tbb is disabled
using concurrent_map = std::map<K, V>;
#endif

using namespace manifold;

template <>
struct std::hash<std::pair<int, int>> {
  size_t operator()(const std::pair<int, int>& p) const {
    return std::hash<int>()(p.first) ^ std::hash<int>()(p.second);
  }
};

namespace {

constexpr int kParallelThreshold = 128;

struct AbsSum {
  int operator()(int a, int b) const { return abs(a) + abs(b); }
};

struct DuplicateVerts {
  VecView<vec3> vertPosR;
  VecView<const int> inclusion;
  VecView<const int> vertR;
  VecView<const vec3> vertPosP;

  void operator()(const int vert) {
    const int n = std::abs(inclusion[vert]);
    for (int i = 0; i < n; ++i) {
      vertPosR[vertR[vert] + i] = vertPosP[vert];
    }
  }
};

template <bool atomic>
struct CountVerts {
  VecView<Halfedge> halfedges;
  VecView<int> count;
  VecView<const int> inclusion;

  void operator()(size_t i) {
    if (atomic)
      AtomicAdd(count[i / 3], std::abs(inclusion[halfedges[i].startVert]));
    else
      count[i / 3] += std::abs(inclusion[halfedges[i].startVert]);
  }
};

template <const bool inverted, const bool atomic>
struct CountNewVerts {
  VecView<int> countP;
  VecView<int> countQ;
  VecView<const int> i12;
  const Vec<std::array<int, 2>>& pq;
  VecView<const Halfedge> halfedges;

  void operator()(const int idx) {
    int edgeP = pq[idx][inverted ? 1 : 0];
    int faceQ = pq[idx][inverted ? 0 : 1];
    int inclusion = std::abs(i12[idx]);

    if (atomic) {
      AtomicAdd(countQ[faceQ], inclusion);
      const Halfedge half = halfedges[edgeP];
      AtomicAdd(countP[edgeP / 3], inclusion);
      AtomicAdd(countP[half.pairedHalfedge / 3], inclusion);
    } else {
      countQ[faceQ] += inclusion;
      const Halfedge half = halfedges[edgeP];
      countP[edgeP / 3] += inclusion;
      countP[half.pairedHalfedge / 3] += inclusion;
    }
  }
};

std::tuple<Vec<int>, Vec<int>> SizeOutput(
    Manifold::Impl& outR, const Manifold::Impl& inP, const Manifold::Impl& inQ,
    const Vec<int>& i03, const Vec<int>& i30, const Vec<int>& i12,
    const Vec<int>& i21, const Vec<std::array<int, 2>>& p1q2,
    const Vec<std::array<int, 2>>& p2q1, bool invertQ) {
  ZoneScoped;
  Vec<int> sidesPerFacePQ(inP.NumTri() + inQ.NumTri(), 0);
  // note: numFaceR <= facePQ2R.size() = sidesPerFacePQ.size() + 1

  auto sidesPerFaceP = sidesPerFacePQ.view(0, inP.NumTri());
  auto sidesPerFaceQ = sidesPerFacePQ.view(inP.NumTri(), inQ.NumTri());

  if (inP.halfedge_.size() >= 1e5) {
    for_each(ExecutionPolicy::Par, countAt(0_uz), countAt(inP.halfedge_.size()),
             CountVerts<true>({inP.halfedge_, sidesPerFaceP, i03}));
    for_each(ExecutionPolicy::Par, countAt(0_uz), countAt(inQ.halfedge_.size()),
             CountVerts<true>({inQ.halfedge_, sidesPerFaceQ, i30}));
  } else {
    for_each(ExecutionPolicy::Seq, countAt(0_uz), countAt(inP.halfedge_.size()),
             CountVerts<false>({inP.halfedge_, sidesPerFaceP, i03}));
    for_each(ExecutionPolicy::Seq, countAt(0_uz), countAt(inQ.halfedge_.size()),
             CountVerts<false>({inQ.halfedge_, sidesPerFaceQ, i30}));
  }

  if (i12.size() >= 1e5) {
    for_each_n(ExecutionPolicy::Par, countAt(0), i12.size(),
               CountNewVerts<false, true>(
                   {sidesPerFaceP, sidesPerFaceQ, i12, p1q2, inP.halfedge_}));
    for_each_n(ExecutionPolicy::Par, countAt(0), i21.size(),
               CountNewVerts<true, true>(
                   {sidesPerFaceQ, sidesPerFaceP, i21, p2q1, inQ.halfedge_}));
  } else {
    for_each_n(ExecutionPolicy::Seq, countAt(0), i12.size(),
               CountNewVerts<false, false>(
                   {sidesPerFaceP, sidesPerFaceQ, i12, p1q2, inP.halfedge_}));
    for_each_n(ExecutionPolicy::Seq, countAt(0), i21.size(),
               CountNewVerts<true, false>(
                   {sidesPerFaceQ, sidesPerFaceP, i21, p2q1, inQ.halfedge_}));
  }

  Vec<int> facePQ2R(inP.NumTri() + inQ.NumTri() + 1, 0);
  auto keepFace = TransformIterator(sidesPerFacePQ.begin(),
                                    [](int x) { return x > 0 ? 1 : 0; });

  inclusive_scan(keepFace, keepFace + sidesPerFacePQ.size(),
                 facePQ2R.begin() + 1);
  int numFaceR = facePQ2R.back();
  facePQ2R.resize(inP.NumTri() + inQ.NumTri());

  outR.faceNormal_.resize_nofill(numFaceR);

  Vec<size_t> tmpBuffer(outR.faceNormal_.size());
  auto faceIds = TransformIterator(countAt(0_uz), [&sidesPerFacePQ](size_t i) {
    if (sidesPerFacePQ[i] > 0) return i;
    return std::numeric_limits<size_t>::max();
  });

  auto next =
      copy_if(faceIds, faceIds + inP.faceNormal_.size(), tmpBuffer.begin(),
              [](size_t v) { return v != std::numeric_limits<size_t>::max(); });

  gather(tmpBuffer.begin(), next, inP.faceNormal_.begin(),
         outR.faceNormal_.begin());

  auto faceIdsQ =
      TransformIterator(countAt(0_uz), [&sidesPerFacePQ, &inP](size_t i) {
        if (sidesPerFacePQ[i + inP.faceNormal_.size()] > 0) return i;
        return std::numeric_limits<size_t>::max();
      });
  auto end =
      copy_if(faceIdsQ, faceIdsQ + inQ.faceNormal_.size(), next,
              [](size_t v) { return v != std::numeric_limits<size_t>::max(); });

  if (invertQ) {
    gather(next, end,
           TransformIterator(inQ.faceNormal_.begin(), Negate<vec3>()),
           outR.faceNormal_.begin() + std::distance(tmpBuffer.begin(), next));
  } else {
    gather(next, end, inQ.faceNormal_.begin(),
           outR.faceNormal_.begin() + std::distance(tmpBuffer.begin(), next));
  }

  auto newEnd = remove(sidesPerFacePQ.begin(), sidesPerFacePQ.end(), 0);
  Vec<int> faceEdge(newEnd - sidesPerFacePQ.begin() + 1, 0);
  inclusive_scan(sidesPerFacePQ.begin(), newEnd, faceEdge.begin() + 1);
  outR.halfedge_.resize(faceEdge.back());

  return std::make_tuple(faceEdge, facePQ2R);
}

struct EdgePos {
  double edgePos;
  int vert;
  int collisionId;
  bool isStart;
};

// thread sanitizer doesn't really know how to check when there are too many
// mutex
#if defined(__has_feature)
#if __has_feature(thread_sanitizer)
__attribute__((no_sanitize("thread")))
#endif
#endif
void AddNewEdgeVerts(
    // we need concurrent_map because we will be adding things concurrently
    concurrent_map<int, std::vector<EdgePos>> &edgesP,
    concurrent_map<std::pair<int, int>, std::vector<EdgePos>> &edgesNew,
    const Vec<std::array<int, 2>> &p1q2, const Vec<int> &i12, const Vec<int> &v12R,
    const Vec<Halfedge> &halfedgeP, bool forward, size_t offset) {
  ZoneScoped;
  // For each edge of P that intersects a face of Q (p1q2), add this vertex to
  // P's corresponding edge vector and to the two new edges, which are
  // intersections between the face of Q and the two faces of P attached to the
  // edge. The direction and duplicity are given by i12, while v12R remaps to
  // the output vert index. When forward is false, all is reversed.
  auto process = [&](std::function<void(size_t)> lock,
                     std::function<void(size_t)> unlock, size_t i) {
    const int edgeP = p1q2[i][forward ? 0 : 1];
    const int faceQ = p1q2[i][forward ? 1 : 0];
    const int vert = v12R[i];
    const int inclusion = i12[i];

    Halfedge halfedge = halfedgeP[edgeP];
    std::pair<int, int> keyRight = {halfedge.pairedHalfedge / 3, faceQ};
    if (!forward) std::swap(keyRight.first, keyRight.second);

    std::pair<int, int> keyLeft = {edgeP / 3, faceQ};
    if (!forward) std::swap(keyLeft.first, keyLeft.second);

    bool direction = inclusion < 0;
    std::hash<std::pair<int, int>> pairHasher;
    std::array<std::tuple<bool, size_t, std::vector<EdgePos>*>, 3> edges = {
        std::make_tuple(direction, std::hash<int>{}(edgeP), &edgesP[edgeP]),
        std::make_tuple(direction ^ !forward,  // revert if not forward
                        pairHasher(keyRight), &edgesNew[keyRight]),
        std::make_tuple(direction ^ forward,  // revert if forward
                        pairHasher(keyLeft), &edgesNew[keyLeft])};
    for (const auto& tuple : edges) {
      lock(std::get<1>(tuple));
      for (int j = 0; j < std::abs(inclusion); ++j)
        std::get<2>(tuple)->push_back(
            {0.0, vert + j, static_cast<int>(i + offset), std::get<0>(tuple)});
      unlock(std::get<1>(tuple));
      direction = !direction;
    }
  };
#if (MANIFOLD_PAR == 1) && __has_include(<tbb/concurrent_map.h>)
  // parallelize operations, requires concurrent_map so we can only enable this
  // with tbb
  if (p1q2.size() > kParallelThreshold) {
    // ideally we should have 1 mutex per key, but kParallelThreshold is enough
    // to avoid contention for most of the cases
    std::array<std::mutex, kParallelThreshold> mutexes;
    static tbb::affinity_partitioner ap;
    auto processFun = std::bind(
        process, [&](size_t hash) { mutexes[hash % mutexes.size()].lock(); },
        [&](size_t hash) { mutexes[hash % mutexes.size()].unlock(); },
        std::placeholders::_1);
    tbb::parallel_for(
        tbb::blocked_range<size_t>(0_uz, p1q2.size(), 32),
        [&](const tbb::blocked_range<size_t>& range) {
          for (size_t i = range.begin(); i != range.end(); i++) processFun(i);
        },
        ap);
    return;
  }
#endif
  auto processFun =
      std::bind(process, [](size_t) {}, [](size_t) {}, std::placeholders::_1);
  for (size_t i = 0; i < p1q2.size(); ++i) processFun(i);
}

std::vector<Halfedge> PairUp(std::vector<EdgePos>& edgePos) {
  // Pair start vertices with end vertices to form edges. The choice of pairing
  // is arbitrary for the manifoldness guarantee, but must be ordered to be
  // geometrically valid. If the order does not go start-end-start-end... then
  // the input and output are not geometrically valid and this algorithm becomes
  // a heuristic.
  DEBUG_ASSERT(edgePos.size() % 2 == 0, topologyErr,
               "Non-manifold edge! Not an even number of points.");
  size_t nEdges = edgePos.size() / 2;
  auto middle = std::partition(edgePos.begin(), edgePos.end(),
                               [](EdgePos x) { return x.isStart; });
  DEBUG_ASSERT(static_cast<size_t>(middle - edgePos.begin()) == nEdges,
               topologyErr, "Non-manifold edge!");
  auto cmp = [](EdgePos a, EdgePos b) {
    return a.edgePos < b.edgePos ||
           // we also sort by collisionId to make things deterministic
           (a.edgePos == b.edgePos && a.collisionId < b.collisionId);
  };
  std::stable_sort(edgePos.begin(), middle, cmp);
  std::stable_sort(middle, edgePos.end(), cmp);
  std::vector<Halfedge> edges;
  for (size_t i = 0; i < nEdges; ++i)
    edges.push_back({edgePos[i].vert, edgePos[i + nEdges].vert, -1});
  return edges;
}

void AppendPartialEdges(Manifold::Impl& outR, Vec<char>& wholeHalfedgeP,
                        Vec<int>& facePtrR,
                        concurrent_map<int, std::vector<EdgePos>>& edgesP,
                        Vec<TriRef>& halfedgeRef, const Manifold::Impl& inP,
                        const Vec<int>& i03, const Vec<int>& vP2R,
                        const Vec<int>::IterC faceP2R, bool forward) {
  ZoneScoped;
  // Each edge in the map is partially retained; for each of these, look up
  // their original verts and include them based on their winding number (i03),
  // while remapping them to the output using vP2R. Use the verts position
  // projected along the edge vector to pair them up, then distribute these
  // edges to their faces.
  Vec<Halfedge>& halfedgeR = outR.halfedge_;
  const Vec<vec3>& vertPosP = inP.vertPos_;
  const Vec<Halfedge>& halfedgeP = inP.halfedge_;

  for (auto& value : edgesP) {
    const int edgeP = value.first;
    std::vector<EdgePos>& edgePosP = value.second;

    const Halfedge& halfedge = halfedgeP[edgeP];
    wholeHalfedgeP[edgeP] = false;
    wholeHalfedgeP[halfedge.pairedHalfedge] = false;

    const int vStart = halfedge.startVert;
    const int vEnd = halfedge.endVert;
    const vec3 edgeVec = vertPosP[vEnd] - vertPosP[vStart];
    // Fill in the edge positions of the old points.
    for (EdgePos& edge : edgePosP) {
      edge.edgePos = la::dot(outR.vertPos_[edge.vert], edgeVec);
    }

    int inclusion = i03[vStart];
    EdgePos edgePos = {la::dot(outR.vertPos_[vP2R[vStart]], edgeVec),
                       vP2R[vStart], std::numeric_limits<int>::max(),
                       inclusion > 0};
    for (int j = 0; j < std::abs(inclusion); ++j) {
      edgePosP.push_back(edgePos);
      ++edgePos.vert;
    }

    inclusion = i03[vEnd];
    edgePos = {la::dot(outR.vertPos_[vP2R[vEnd]], edgeVec), vP2R[vEnd],
               std::numeric_limits<int>::max(), inclusion < 0};
    for (int j = 0; j < std::abs(inclusion); ++j) {
      edgePosP.push_back(edgePos);
      ++edgePos.vert;
    }

    // sort edges into start/end pairs along length
    std::vector<Halfedge> edges = PairUp(edgePosP);

    // add halfedges to result
    const int faceLeftP = edgeP / 3;
    const int faceLeft = faceP2R[faceLeftP];
    const int faceRightP = halfedge.pairedHalfedge / 3;
    const int faceRight = faceP2R[faceRightP];
    // Negative inclusion means the halfedges are reversed, which means our
    // reference is now to the endVert instead of the startVert, which is one
    // position advanced CCW. This is only valid if this is a retained vert; it
    // will be ignored later if the vert is new.
    const TriRef forwardRef = {forward ? 0 : 1, -1, faceLeftP, -1};
    const TriRef backwardRef = {forward ? 0 : 1, -1, faceRightP, -1};

    for (Halfedge e : edges) {
      const int forwardEdge = facePtrR[faceLeft]++;
      const int backwardEdge = facePtrR[faceRight]++;

      e.pairedHalfedge = backwardEdge;
      halfedgeR[forwardEdge] = e;
      halfedgeRef[forwardEdge] = forwardRef;

      std::swap(e.startVert, e.endVert);
      e.pairedHalfedge = forwardEdge;
      halfedgeR[backwardEdge] = e;
      halfedgeRef[backwardEdge] = backwardRef;
    }
  }
}

void AppendNewEdges(
    Manifold::Impl& outR, Vec<int>& facePtrR,
    concurrent_map<std::pair<int, int>, std::vector<EdgePos>>& edgesNew,
    Vec<TriRef>& halfedgeRef, const Vec<int>& facePQ2R, const int numFaceP) {
  ZoneScoped;
  // Pair up each edge's verts and distribute to faces based on indices in key.
  Vec<Halfedge>& halfedgeR = outR.halfedge_;
  Vec<vec3>& vertPosR = outR.vertPos_;

  for (auto& value : edgesNew) {
    const int faceP = value.first.first;
    const int faceQ = value.first.second;
    std::vector<EdgePos>& edgePos = value.second;

    Box bbox;
    for (auto edge : edgePos) {
      bbox.Union(vertPosR[edge.vert]);
    }
    const vec3 size = bbox.Size();
    // Order the points along their longest dimension.
    const int i = (size.x > size.y && size.x > size.z) ? 0
                  : size.y > size.z                    ? 1
                                                       : 2;
    for (auto& edge : edgePos) {
      edge.edgePos = vertPosR[edge.vert][i];
    }

    // sort edges into start/end pairs along length.
    std::vector<Halfedge> edges = PairUp(edgePos);

    // add halfedges to result
    const int faceLeft = facePQ2R[faceP];
    const int faceRight = facePQ2R[numFaceP + faceQ];
    const TriRef forwardRef = {0, -1, faceP, -1};
    const TriRef backwardRef = {1, -1, faceQ, -1};
    for (Halfedge e : edges) {
      const int forwardEdge = facePtrR[faceLeft]++;
      const int backwardEdge = facePtrR[faceRight]++;

      e.pairedHalfedge = backwardEdge;
      halfedgeR[forwardEdge] = e;
      halfedgeRef[forwardEdge] = forwardRef;

      std::swap(e.startVert, e.endVert);
      e.pairedHalfedge = forwardEdge;
      halfedgeR[backwardEdge] = e;
      halfedgeRef[backwardEdge] = backwardRef;
    }
  }
}

struct DuplicateHalfedges {
  VecView<Halfedge> halfedgesR;
  VecView<TriRef> halfedgeRef;
  VecView<int> facePtr;
  VecView<const char> wholeHalfedgeP;
  VecView<const Halfedge> halfedgesP;
  VecView<const int> i03;
  VecView<const int> vP2R;
  VecView<const int> faceP2R;
  const bool forward;

  void operator()(const int idx) {
    if (!wholeHalfedgeP[idx]) return;
    Halfedge halfedge = halfedgesP[idx];
    if (!halfedge.IsForward()) return;

    const int inclusion = i03[halfedge.startVert];
    if (inclusion == 0) return;
    if (inclusion < 0) {  // reverse
      std::swap(halfedge.startVert, halfedge.endVert);
    }
    halfedge.startVert = vP2R[halfedge.startVert];
    halfedge.endVert = vP2R[halfedge.endVert];
    const int faceLeftP = idx / 3;
    const int newFace = faceP2R[faceLeftP];
    const int faceRightP = halfedge.pairedHalfedge / 3;
    const int faceRight = faceP2R[faceRightP];
    // Negative inclusion means the halfedges are reversed, which means our
    // reference is now to the endVert instead of the startVert, which is one
    // position advanced CCW.
    const TriRef forwardRef = {forward ? 0 : 1, -1, faceLeftP, -1};
    const TriRef backwardRef = {forward ? 0 : 1, -1, faceRightP, -1};

    for (int i = 0; i < std::abs(inclusion); ++i) {
      int forwardEdge = AtomicAdd(facePtr[newFace], 1);
      int backwardEdge = AtomicAdd(facePtr[faceRight], 1);
      halfedge.pairedHalfedge = backwardEdge;

      halfedgesR[forwardEdge] = halfedge;
      halfedgesR[backwardEdge] = {halfedge.endVert, halfedge.startVert,
                                  forwardEdge};
      halfedgeRef[forwardEdge] = forwardRef;
      halfedgeRef[backwardEdge] = backwardRef;

      ++halfedge.startVert;
      ++halfedge.endVert;
    }
  }
};

void AppendWholeEdges(Manifold::Impl& outR, Vec<int>& facePtrR,
                      Vec<TriRef>& halfedgeRef, const Manifold::Impl& inP,
                      const Vec<char> wholeHalfedgeP, const Vec<int>& i03,
                      const Vec<int>& vP2R, VecView<const int> faceP2R,
                      bool forward) {
  ZoneScoped;
  for_each_n(
      autoPolicy(inP.halfedge_.size()), countAt(0), inP.halfedge_.size(),
      DuplicateHalfedges({outR.halfedge_, halfedgeRef, facePtrR, wholeHalfedgeP,
                          inP.halfedge_, i03, vP2R, faceP2R, forward}));
}

struct MapTriRef {
  VecView<const TriRef> triRefP;
  VecView<const TriRef> triRefQ;
  const int offsetQ;

  void operator()(TriRef& triRef) {
    const int tri = triRef.faceID;
    const bool PQ = triRef.meshID == 0;
    triRef = PQ ? triRefP[tri] : triRefQ[tri];
    if (!PQ) triRef.meshID += offsetQ;
  }
};

void UpdateReference(Manifold::Impl& outR, const Manifold::Impl& inP,
                     const Manifold::Impl& inQ, bool invertQ) {
  const int offsetQ = Manifold::Impl::meshIDCounter_;
  for_each_n(
      autoPolicy(outR.NumTri(), 1e5), outR.meshRelation_.triRef.begin(),
      outR.NumTri(),
      MapTriRef({inP.meshRelation_.triRef, inQ.meshRelation_.triRef, offsetQ}));

  for (const auto& pair : inP.meshRelation_.meshIDtransform) {
    outR.meshRelation_.meshIDtransform[pair.first] = pair.second;
  }
  for (const auto& pair : inQ.meshRelation_.meshIDtransform) {
    outR.meshRelation_.meshIDtransform[pair.first + offsetQ] = pair.second;
    outR.meshRelation_.meshIDtransform[pair.first + offsetQ].backSide ^=
        invertQ;
  }
}

struct Barycentric {
  VecView<vec3> uvw;
  VecView<const TriRef> ref;
  VecView<const vec3> vertPosP;
  VecView<const vec3> vertPosQ;
  VecView<const vec3> vertPosR;
  VecView<const Halfedge> halfedgeP;
  VecView<const Halfedge> halfedgeQ;
  VecView<const Halfedge> halfedgeR;
  const double epsilon;

  void operator()(const int tri) {
    const TriRef refPQ = ref[tri];
    if (halfedgeR[3 * tri].startVert < 0) return;

    const int triPQ = refPQ.faceID;
    const bool PQ = refPQ.meshID == 0;
    const auto& vertPos = PQ ? vertPosP : vertPosQ;
    const auto& halfedge = PQ ? halfedgeP : halfedgeQ;

    mat3 triPos;
    for (const int j : {0, 1, 2})
      triPos[j] = vertPos[halfedge[3 * triPQ + j].startVert];

    for (const int i : {0, 1, 2}) {
      const int vert = halfedgeR[3 * tri + i].startVert;
      uvw[3 * tri + i] = GetBarycentric(vertPosR[vert], triPos, epsilon);
    }
  }
};

void CreateProperties(Manifold::Impl& outR, const Manifold::Impl& inP,
                      const Manifold::Impl& inQ) {
  ZoneScoped;
  const int numPropP = inP.NumProp();
  const int numPropQ = inQ.NumProp();
  const int numProp = std::max(numPropP, numPropQ);
  outR.numProp_ = numProp;
  if (numProp == 0) return;

  const int numTri = outR.NumTri();
  Vec<vec3> bary(outR.halfedge_.size());
  for_each_n(autoPolicy(numTri, 1e4), countAt(0), numTri,
             Barycentric({bary, outR.meshRelation_.triRef, inP.vertPos_,
                          inQ.vertPos_, outR.vertPos_, inP.halfedge_,
                          inQ.halfedge_, outR.halfedge_, outR.epsilon_}));

  using Entry = std::pair<ivec3, int>;
  int idMissProp = outR.NumVert();
  std::vector<std::vector<Entry>> propIdx(outR.NumVert() + 1);
  std::vector<int> propMissIdx[2];
  propMissIdx[0].resize(inQ.NumPropVert(), -1);
  propMissIdx[1].resize(inP.NumPropVert(), -1);

  outR.properties_.reserve(outR.NumVert() * numProp);
  int idx = 0;

  for (int tri = 0; tri < numTri; ++tri) {
    // Skip collapsed triangles
    if (outR.halfedge_[3 * tri].startVert < 0) continue;

    const TriRef ref = outR.meshRelation_.triRef[tri];
    const bool PQ = ref.meshID == 0;
    const int oldNumProp = PQ ? numPropP : numPropQ;
    const auto& properties = PQ ? inP.properties_ : inQ.properties_;
    const auto& halfedge = PQ ? inP.halfedge_ : inQ.halfedge_;

    for (const int i : {0, 1, 2}) {
      const int vert = outR.halfedge_[3 * tri + i].startVert;
      const vec3& uvw = bary[3 * tri + i];

      ivec4 key(PQ, idMissProp, -1, -1);
      if (oldNumProp > 0) {
        int edge = -2;
        for (const int j : {0, 1, 2}) {
          if (uvw[j] == 1) {
            // On a retained vert, the propVert must also match
            key[2] = halfedge[3 * ref.faceID + j].propVert;
            edge = -1;
            break;
          }
          if (uvw[j] == 0) edge = j;
        }
        if (edge >= 0) {
          // On an edge, both propVerts must match
          const int p0 = halfedge[3 * ref.faceID + Next3(edge)].propVert;
          const int p1 = halfedge[3 * ref.faceID + Prev3(edge)].propVert;
          key[1] = vert;
          key[2] = std::min(p0, p1);
          key[3] = std::max(p0, p1);
        } else if (edge == -2) {
          key[1] = vert;
        }
      }

      if (key.y == idMissProp && key.z >= 0) {
        // only key.x/key.z matters
        auto& entry = propMissIdx[key.x][key.z];
        if (entry >= 0) {
          outR.halfedge_[3 * tri + i].propVert = entry;
          continue;
        }
        entry = idx;
      } else {
        auto& bin = propIdx[key.y];
        bool bFound = false;
        for (const auto& b : bin) {
          if (b.first == ivec3(key.x, key.z, key.w)) {
            bFound = true;
            outR.halfedge_[3 * tri + i].propVert = b.second;
            break;
          }
        }
        if (bFound) continue;
        bin.push_back(std::make_pair(ivec3(key.x, key.z, key.w), idx));
      }

      outR.halfedge_[3 * tri + i].propVert = idx++;
      for (int p = 0; p < numProp; ++p) {
        if (p < oldNumProp) {
          vec3 oldProps;
          for (const int j : {0, 1, 2})
            oldProps[j] =
                properties[oldNumProp * halfedge[3 * ref.faceID + j].propVert +
                           p];
          outR.properties_.push_back(la::dot(uvw, oldProps));
        } else {
          outR.properties_.push_back(0);
        }
      }
    }
  }
}

void ReorderHalfedges(VecView<Halfedge>& halfedges) {
  ZoneScoped;
  // halfedges in the same face are added in non-deterministic order, so we have
  // to reorder them for determinism

  // step 1: reorder within the same face, such that the halfedge with the
  // smallest starting vertex is placed first
  for_each(autoPolicy(halfedges.size() / 3), countAt(0),
           countAt(halfedges.size() / 3), [&halfedges](size_t tri) {
             std::array<Halfedge, 3> face = {halfedges[tri * 3],
                                             halfedges[tri * 3 + 1],
                                             halfedges[tri * 3 + 2]};
             int index = 0;
             for (int i : {1, 2})
               if (face[i].startVert < face[index].startVert) index = i;
             for (int i : {0, 1, 2})
               halfedges[tri * 3 + i] = face[(index + i) % 3];
           });
  // step 2: fix paired halfedge
  for_each(autoPolicy(halfedges.size() / 3), countAt(0),
           countAt(halfedges.size() / 3), [&halfedges](size_t tri) {
             for (int i : {0, 1, 2}) {
               Halfedge& curr = halfedges[tri * 3 + i];
               int oppositeFace = curr.pairedHalfedge / 3;
               int index = -1;
               for (int j : {0, 1, 2})
                 if (curr.startVert == halfedges[oppositeFace * 3 + j].endVert)
                   index = j;
               curr.pairedHalfedge = oppositeFace * 3 + index;
             }
           });
}

}  // namespace

namespace manifold {

Manifold::Impl Boolean3::Result(OpType op) const {
#ifdef MANIFOLD_DEBUG
  Timer assemble;
  assemble.Start();
#endif

  DEBUG_ASSERT((expandP_ > 0) == (op == OpType::Add), logicErr,
               "Result op type not compatible with constructor op type.");
  const int c1 = op == OpType::Intersect ? 0 : 1;
  const int c2 = op == OpType::Add ? 1 : 0;
  const int c3 = op == OpType::Intersect ? 1 : -1;

  if (inP_.status_ != Manifold::Error::NoError) {
    auto impl = Manifold::Impl();
    impl.status_ = inP_.status_;
    return impl;
  }
  if (inQ_.status_ != Manifold::Error::NoError) {
    auto impl = Manifold::Impl();
    impl.status_ = inQ_.status_;
    return impl;
  }

  if (inP_.IsEmpty()) {
    if (!inQ_.IsEmpty() && op == OpType::Add) {
      return inQ_;
    }
    return Manifold::Impl();
  } else if (inQ_.IsEmpty()) {
    if (op == OpType::Intersect) {
      return Manifold::Impl();
    }
    return inP_;
  }

  if (!valid) {
    auto impl = Manifold::Impl();
    impl.status_ = Manifold::Error::ResultTooLarge;
    return impl;
  }

  const bool invertQ = op == OpType::Subtract;

  // Convert winding numbers to inclusion values based on operation type.
  Vec<int> i12(x12_.size());
  Vec<int> i21(x21_.size());
  Vec<int> i03(w03_.size());
  Vec<int> i30(w30_.size());

  transform(x12_.begin(), x12_.end(), i12.begin(),
            [c3](int v) { return c3 * v; });
  transform(x21_.begin(), x21_.end(), i21.begin(),
            [c3](int v) { return c3 * v; });
  transform(w03_.begin(), w03_.end(), i03.begin(),
            [c1, c3](int v) { return c1 + c3 * v; });
  transform(w30_.begin(), w30_.end(), i30.begin(),
            [c2, c3](int v) { return c2 + c3 * v; });

  Vec<int> vP2R(inP_.NumVert());
  exclusive_scan(i03.begin(), i03.end(), vP2R.begin(), 0, AbsSum());
  int numVertR = AbsSum()(vP2R.back(), i03.back());
  const int nPv = numVertR;

  Vec<int> vQ2R(inQ_.NumVert());
  exclusive_scan(i30.begin(), i30.end(), vQ2R.begin(), numVertR, AbsSum());
  numVertR = AbsSum()(vQ2R.back(), i30.back());
  const int nQv = numVertR - nPv;

  Vec<int> v12R(v12_.size());
  if (v12_.size() > 0) {
    exclusive_scan(i12.begin(), i12.end(), v12R.begin(), numVertR, AbsSum());
    numVertR = AbsSum()(v12R.back(), i12.back());
  }
  const int n12 = numVertR - nPv - nQv;

  Vec<int> v21R(v21_.size());
  if (v21_.size() > 0) {
    exclusive_scan(i21.begin(), i21.end(), v21R.begin(), numVertR, AbsSum());
    numVertR = AbsSum()(v21R.back(), i21.back());
  }
  const int n21 = numVertR - nPv - nQv - n12;

  // Create the output Manifold
  Manifold::Impl outR;

  if (numVertR == 0) return outR;

  outR.epsilon_ = std::max(inP_.epsilon_, inQ_.epsilon_);
  outR.tolerance_ = std::max(inP_.tolerance_, inQ_.tolerance_);

  outR.vertPos_.resize_nofill(numVertR);
  // Add vertices, duplicating for inclusion numbers not in [-1, 1].
  // Retained vertices from P and Q:
  for_each_n(autoPolicy(inP_.NumVert(), 1e4), countAt(0), inP_.NumVert(),
             DuplicateVerts({outR.vertPos_, i03, vP2R, inP_.vertPos_}));
  for_each_n(autoPolicy(inQ_.NumVert(), 1e4), countAt(0), inQ_.NumVert(),
             DuplicateVerts({outR.vertPos_, i30, vQ2R, inQ_.vertPos_}));
  // New vertices created from intersections:
  for_each_n(autoPolicy(i12.size(), 1e4), countAt(0), i12.size(),
             DuplicateVerts({outR.vertPos_, i12, v12R, v12_}));
  for_each_n(autoPolicy(i21.size(), 1e4), countAt(0), i21.size(),
             DuplicateVerts({outR.vertPos_, i21, v21R, v21_}));

  PRINT(nPv << " verts from inP");
  PRINT(nQv << " verts from inQ");
  PRINT(n12 << " new verts from edgesP -> facesQ");
  PRINT(n21 << " new verts from facesP -> edgesQ");

  // Build up new polygonal faces from triangle intersections. At this point the
  // calculation switches from parallel to serial.

  // Level 3

  // This key is the forward halfedge index of P or Q. Only includes intersected
  // edges.
  concurrent_map<int, std::vector<EdgePos>> edgesP, edgesQ;
  // This key is the face index of <P, Q>
  concurrent_map<std::pair<int, int>, std::vector<EdgePos>> edgesNew;

  AddNewEdgeVerts(edgesP, edgesNew, p1q2_, i12, v12R, inP_.halfedge_, true, 0);
  AddNewEdgeVerts(edgesQ, edgesNew, p2q1_, i21, v21R, inQ_.halfedge_, false,
                  p1q2_.size());

  v12R.clear();
  v21R.clear();

  // Level 4
  Vec<int> faceEdge;
  Vec<int> facePQ2R;
  std::tie(faceEdge, facePQ2R) =
      SizeOutput(outR, inP_, inQ_, i03, i30, i12, i21, p1q2_, p2q1_, invertQ);

  i12.clear();
  i21.clear();

  // This gets incremented for each halfedge that's added to a face so that the
  // next one knows where to slot in.
  Vec<int> facePtrR = faceEdge;
  // Intersected halfedges are marked false.
  Vec<char> wholeHalfedgeP(inP_.halfedge_.size(), true);
  Vec<char> wholeHalfedgeQ(inQ_.halfedge_.size(), true);
  // The halfedgeRef contains the data that will become triRef once the faces
  // are triangulated.
  Vec<TriRef> halfedgeRef(2 * outR.NumEdge());

  AppendPartialEdges(outR, wholeHalfedgeP, facePtrR, edgesP, halfedgeRef, inP_,
                     i03, vP2R, facePQ2R.begin(), true);
  AppendPartialEdges(outR, wholeHalfedgeQ, facePtrR, edgesQ, halfedgeRef, inQ_,
                     i30, vQ2R, facePQ2R.begin() + inP_.NumTri(), false);

  edgesP.clear();
  edgesQ.clear();

  AppendNewEdges(outR, facePtrR, edgesNew, halfedgeRef, facePQ2R,
                 inP_.NumTri());

  edgesNew.clear();

  AppendWholeEdges(outR, facePtrR, halfedgeRef, inP_, wholeHalfedgeP, i03, vP2R,
                   facePQ2R.cview(0, inP_.NumTri()), true);
  AppendWholeEdges(outR, facePtrR, halfedgeRef, inQ_, wholeHalfedgeQ, i30, vQ2R,
                   facePQ2R.cview(inP_.NumTri(), inQ_.NumTri()), false);

  wholeHalfedgeP.clear();
  wholeHalfedgeQ.clear();
  facePtrR.clear();
  facePQ2R.clear();
  i03.clear();
  i30.clear();
  vP2R.clear();
  vQ2R.clear();

#ifdef MANIFOLD_DEBUG
  assemble.Stop();
  Timer triangulate;
  triangulate.Start();
#endif

  // Level 6

  if (ManifoldParams().intermediateChecks)
    DEBUG_ASSERT(outR.IsManifold(), logicErr, "polygon mesh is not manifold!");

  outR.Face2Tri(faceEdge, halfedgeRef);
  halfedgeRef.clear();
  faceEdge.clear();

  ReorderHalfedges(outR.halfedge_);

#ifdef MANIFOLD_DEBUG
  triangulate.Stop();
  Timer simplify;
  simplify.Start();
#endif

  if (ManifoldParams().intermediateChecks)
    DEBUG_ASSERT(outR.IsManifold(), logicErr,
                 "triangulated mesh is not manifold!");

  CreateProperties(outR, inP_, inQ_);

  UpdateReference(outR, inP_, inQ_, invertQ);

  outR.SimplifyTopology(nPv + nQv);
  outR.RemoveUnreferencedVerts();

  if (ManifoldParams().intermediateChecks)
    DEBUG_ASSERT(outR.Is2Manifold(), logicErr,
                 "simplified mesh is not 2-manifold!");

#ifdef MANIFOLD_DEBUG
  simplify.Stop();
  Timer sort;
  sort.Start();
#endif

  outR.Finish();
  outR.IncrementMeshIDs();

#ifdef MANIFOLD_DEBUG
  sort.Stop();
  if (ManifoldParams().verbose) {
    assemble.Print("Assembly");
    triangulate.Print("Triangulation");
    simplify.Print("Simplification");
    sort.Print("Sorting");
    std::cout << outR.NumVert() << " verts and " << outR.NumTri() << " tris"
              << std::endl;
  }
#endif

  return outR;
}

}  // namespace manifold