diff --git a/thirdparty/README.md b/thirdparty/README.md
index dc7eee4a0d67..0110e6daf203 100644
--- a/thirdparty/README.md
+++ b/thirdparty/README.md
@@ -567,12 +567,12 @@ in the MSVC debugger.
 ## manifold
 
 - Upstream: https://github.com/elalish/manifold
-- Version: master (36035428bc32302a9d7c9ee1ecc833fb8394a2a3, 2024)
+- Version: 3.0.1 (98b8142519d35c13e0e25cfa9fd6e3a271403be6, 2024)
 - License: Apache 2.0
 
 File extracted from upstream source:
 
-- `src/`
+- `src/` and `include/`, except from `CMakeLists.txt`, `cross_section.cpp` and `meshIO.{cpp,h}`
 - `AUTHORS`, `LICENSE`
 
 
diff --git a/thirdparty/manifold/include/manifold/manifold.h b/thirdparty/manifold/include/manifold/manifold.h
index a4f65e9c4488..9c59d7367af8 100644
--- a/thirdparty/manifold/include/manifold/manifold.h
+++ b/thirdparty/manifold/include/manifold/manifold.h
@@ -16,6 +16,10 @@
 #include <functional>
 #include <memory>
 
+#ifdef MANIFOLD_EXPORT
+#include <iostream>
+#endif
+
 #include "manifold/common.h"
 #include "manifold/vec_view.h"
 
@@ -376,6 +380,10 @@ class Manifold {
 
   struct Impl;
 
+#ifdef MANIFOLD_EXPORT
+  static Manifold ImportMeshGL64(std::istream& stream);
+#endif
+
  private:
   Manifold(std::shared_ptr<CsgNode> pNode_);
   Manifold(std::shared_ptr<Impl> pImpl_);
diff --git a/thirdparty/manifold/src/boolean_result.cpp b/thirdparty/manifold/src/boolean_result.cpp
index 44faf73e97b8..2921dde90689 100644
--- a/thirdparty/manifold/src/boolean_result.cpp
+++ b/thirdparty/manifold/src/boolean_result.cpp
@@ -857,6 +857,7 @@ Manifold::Impl Boolean3::Result(OpType op) const {
   UpdateReference(outR, inP_, inQ_, invertQ);
 
   outR.SimplifyTopology();
+  outR.RemoveUnreferencedVerts();
 
   if (ManifoldParams().intermediateChecks)
     DEBUG_ASSERT(outR.Is2Manifold(), logicErr,
diff --git a/thirdparty/manifold/src/cross_section/cross_section.cpp b/thirdparty/manifold/src/cross_section/cross_section.cpp
index ac45fc5233ec..bacdfd10e6b8 100644
--- a/thirdparty/manifold/src/cross_section/cross_section.cpp
+++ b/thirdparty/manifold/src/cross_section/cross_section.cpp
@@ -740,8 +740,8 @@ double CrossSection::Area() const { return C2::Area(GetPaths()->paths_); }
 /**
  * Return the number of vertices in the CrossSection.
  */
-int CrossSection::NumVert() const {
-  int n = 0;
+size_t CrossSection::NumVert() const {
+  size_t n = 0;
   auto paths = GetPaths()->paths_;
   for (auto p : paths) {
     n += p.size();
@@ -753,7 +753,7 @@ int CrossSection::NumVert() const {
  * Return the number of contours (both outer and inner paths) in the
  * CrossSection.
  */
-int CrossSection::NumContour() const { return GetPaths()->paths_.size(); }
+size_t CrossSection::NumContour() const { return GetPaths()->paths_.size(); }
 
 /**
  * Does the CrossSection contain any contours?
diff --git a/thirdparty/manifold/src/csg_tree.cpp b/thirdparty/manifold/src/csg_tree.cpp
index a05fce3aa6b6..e126485ecef2 100644
--- a/thirdparty/manifold/src/csg_tree.cpp
+++ b/thirdparty/manifold/src/csg_tree.cpp
@@ -112,6 +112,50 @@ std::shared_ptr<CsgLeafNode> ImplToLeaf(Manifold::Impl &&impl) {
   return std::make_shared<CsgLeafNode>(std::make_shared<Manifold::Impl>(impl));
 }
 
+std::shared_ptr<CsgLeafNode> SimpleBoolean(const Manifold::Impl &a,
+                                           const Manifold::Impl &b, OpType op) {
+#ifdef MANIFOLD_DEBUG
+  auto dump = [&]() {
+    dump_lock.lock();
+    std::cout << "LHS self-intersecting: " << a.IsSelfIntersecting()
+              << std::endl;
+    std::cout << "RHS self-intersecting: " << b.IsSelfIntersecting()
+              << std::endl;
+    if (ManifoldParams().verbose) {
+      if (op == OpType::Add)
+        std::cout << "Add";
+      else if (op == OpType::Intersect)
+        std::cout << "Intersect";
+      else
+        std::cout << "Subtract";
+      std::cout << std::endl;
+      std::cout << a;
+      std::cout << b;
+    }
+    dump_lock.unlock();
+  };
+  try {
+    Boolean3 boolean(a, b, op);
+    auto impl = boolean.Result(op);
+    if (ManifoldParams().intermediateChecks && impl.IsSelfIntersecting()) {
+      dump_lock.lock();
+      std::cout << "self intersections detected" << std::endl;
+      dump_lock.unlock();
+      throw logicErr("self intersection detected");
+    }
+    return ImplToLeaf(std::move(impl));
+  } catch (logicErr &err) {
+    dump();
+    throw err;
+  } catch (geometryErr &err) {
+    dump();
+    throw err;
+  }
+#else
+  return ImplToLeaf(Boolean3(a, b, op).Result(op));
+#endif
+}
+
 /**
  * Efficient union of a set of pairwise disjoint meshes.
  */
@@ -306,10 +350,9 @@ std::shared_ptr<CsgLeafNode> BatchBoolean(
   // common cases
   if (results.size() == 0) return std::make_shared<CsgLeafNode>();
   if (results.size() == 1) return results.front();
-  if (results.size() == 2) {
-    Boolean3 boolean(*results[0]->GetImpl(), *results[1]->GetImpl(), operation);
-    return ImplToLeaf(boolean.Result(operation));
-  }
+  if (results.size() == 2)
+    return SimpleBoolean(*results[0]->GetImpl(), *results[1]->GetImpl(),
+                         operation);
 #if (MANIFOLD_PAR == 1) && __has_include(<tbb/tbb.h>)
   tbb::task_group group;
   tbb::concurrent_priority_queue<std::shared_ptr<CsgLeafNode>, MeshCompare>
@@ -327,8 +370,7 @@ std::shared_ptr<CsgLeafNode> BatchBoolean(
         continue;
       }
       group.run([&, a, b]() {
-        Boolean3 boolean(*a->GetImpl(), *b->GetImpl(), operation);
-        queue.emplace(ImplToLeaf(boolean.Result(operation)));
+        queue.emplace(SimpleBoolean(*a->GetImpl(), *b->GetImpl(), operation));
         return group.run(process);
       });
     }
@@ -351,11 +393,8 @@ std::shared_ptr<CsgLeafNode> BatchBoolean(
     auto b = std::move(results.back());
     results.pop_back();
     // boolean operation
-    Boolean3 boolean(*a->GetImpl(), *b->GetImpl(), operation);
-    auto result = ImplToLeaf(boolean.Result(operation));
-    if (results.size() == 0) {
-      return result;
-    }
+    auto result = SimpleBoolean(*a->GetImpl(), *b->GetImpl(), operation);
+    if (results.size() == 0) return result;
     results.push_back(result);
     std::push_heap(results.begin(), results.end(), cmpFn);
   }
@@ -429,9 +468,33 @@ CsgOpNode::CsgOpNode() {}
 
 CsgOpNode::CsgOpNode(const std::vector<std::shared_ptr<CsgNode>> &children,
                      OpType op)
-    : impl_(Impl{}), op_(op) {
-  auto impl = impl_.GetGuard();
-  impl->children_ = children;
+    : impl_(children), op_(op) {}
+
+CsgOpNode::~CsgOpNode() {
+  if (impl_.UseCount() == 1) {
+    auto impl = impl_.GetGuard();
+    std::vector<std::shared_ptr<CsgOpNode>> toProcess;
+    auto handleChildren =
+        [&toProcess](std::vector<std::shared_ptr<CsgNode>> &children) {
+          while (!children.empty()) {
+            // move out so shrinking the vector will not trigger recursive drop
+            auto movedChild = std::move(children.back());
+            children.pop_back();
+            if (movedChild->GetNodeType() != CsgNodeType::Leaf)
+              toProcess.push_back(
+                  std::static_pointer_cast<CsgOpNode>(std::move(movedChild)));
+          }
+        };
+    handleChildren(*impl);
+    while (!toProcess.empty()) {
+      auto child = std::move(toProcess.back());
+      toProcess.pop_back();
+      if (impl_.UseCount() == 1) {
+        auto childImpl = child->impl_.GetGuard();
+        handleChildren(*childImpl);
+      }
+    }
+  }
 }
 
 std::shared_ptr<CsgNode> CsgOpNode::Boolean(
@@ -452,21 +515,25 @@ std::shared_ptr<CsgNode> CsgOpNode::Transform(const mat3x4 &m) const {
 }
 
 struct CsgStackFrame {
+  using Nodes = std::vector<std::shared_ptr<CsgLeafNode>>;
+
   bool finalize;
   OpType parent_op;
   mat3x4 transform;
-  std::vector<std::shared_ptr<CsgLeafNode>> *destination;
+  Nodes *positive_dest;
+  Nodes *negative_dest;
   std::shared_ptr<const CsgOpNode> op_node;
-  std::vector<std::shared_ptr<CsgLeafNode>> positive_children;
-  std::vector<std::shared_ptr<CsgLeafNode>> negative_children;
+  Nodes positive_children;
+  Nodes negative_children;
 
   CsgStackFrame(bool finalize, OpType parent_op, mat3x4 transform,
-                std::vector<std::shared_ptr<CsgLeafNode>> *parent,
+                Nodes *positive_dest, Nodes *negative_dest,
                 std::shared_ptr<const CsgOpNode> op_node)
       : finalize(finalize),
         parent_op(parent_op),
         transform(transform),
-        destination(parent),
+        positive_dest(positive_dest),
+        negative_dest(negative_dest),
         op_node(op_node) {}
 };
 
@@ -481,10 +548,10 @@ std::shared_ptr<CsgLeafNode> CsgOpNode::ToLeafNode() const {
   // when we remove it from the stack, but it is a bit more complicated and
   // there is no measurable overhead from using `shared_ptr` here...
   std::vector<std::shared_ptr<CsgStackFrame>> stack;
-  // initial node, destination is a nullptr because we don't need to put the
+  // initial node, positive_dest is a nullptr because we don't need to put the
   // result anywhere else (except in the cache_).
   stack.push_back(std::make_shared<CsgStackFrame>(
-      false, op_, la::identity, nullptr,
+      false, op_, la::identity, nullptr, nullptr,
       std::static_pointer_cast<const CsgOpNode>(shared_from_this())));
 
   // Instead of actually using recursion in the algorithm, we use an explicit
@@ -494,7 +561,7 @@ std::shared_ptr<CsgLeafNode> CsgOpNode::ToLeafNode() const {
   // Before performing boolean operations, we should make sure that all children
   // are `CsgLeafNodes`, i.e. are actual meshes that can be operated on. Hence,
   // we do it in two steps:
-  // 1. Populate `children` (`left_children` and `right_children`, see below)
+  // 1. Populate `children` (`positive_children` and `negative_children`)
   //    If the child is a `CsgOpNode`, we either collapse it or compute its
   //    boolean operation result.
   // 2. Performs boolean after populating the `children` set.
@@ -510,21 +577,23 @@ std::shared_ptr<CsgLeafNode> CsgOpNode::ToLeafNode() const {
   // reused. And in the special case where the children set only contains one
   // element, we don't need any operation, so we can collapse that as well.
   // Instead of moving `b` and `c` into the parent, and running this collapsing
-  // check until a fixed point, we remember the `destination` where we should
-  // put the `CsgLeafNode` into. Normally, the `destination` pointer point to
+  // check until a fixed point, we remember the `positive_dest` where we should
+  // put the `CsgLeafNode` into. Normally, the `positive_dest` pointer point to
   // the parent `children` set. However, when a child is being collapsed, we
-  // keep using the old `destination` pointer for the grandchildren. Hence,
+  // keep using the old `positive_dest` pointer for the grandchildren. Hence,
   // removing a node by collapsing takes O(1) time. We also need to store the
   // parent operation type for checking if the node is eligible for collapsing,
   // and transform matrix because we need to re-apply the transformation to the
   // children.
   //
-  // `Subtract` is handled differently from `Add` and `Intersect`. It is treated
-  // as two `Add` nodes, `positive_children` and `negative_children`, that
-  // should be subtracted later. This allows collapsing children `Add` nodes.
-  // For normal `Add` and `Intersect`, we only use `positive_children`.
+  // `Subtract` is handled differently from `Add` and `Intersect`.
+  // For the first operand, it is treated as normal subtract. Negative children
+  // in this operand is propagated to the parent, which is equivalent to
+  // collapsing `(a - b) - c` into `a - (b + c)`.
+  // For the remaining operands, they are treated as a nested `Add` node,
+  // collapsing `a - (b + (c + d))` into `a - (b + c + d)`.
   //
-  // `impl->children_` should always contain either the raw set of children or
+  // `impl` should always contain either the raw set of children or
   // the NOT transformed result, while `cache_` should contain the transformed
   // result. This is because `impl` can be shared between `CsgOpNode` that
   // differ in `transform_`, so we want it to be able to share the result.
@@ -532,38 +601,39 @@ std::shared_ptr<CsgLeafNode> CsgOpNode::ToLeafNode() const {
   // Recursive version (pseudocode only):
   //
   // void f(CsgOpNode node, OpType parent_op, mat3x4 transform,
-  //        std::vector<CsgLeafNode> *destination) {
+  //        Nodes *positive_dest, Nodes *negative_dest) {
   //   auto impl = node->impl_.GetGuard();
   //   // can collapse when we have the same operation as the parent and is
   //   // unique, or when we have only one children.
-  //   const bool canCollapse = (node->op_ == parent_op && IsUnique(node)) ||
-  //                            impl->children_.size() == 1;
+  //   const OpType op = node->op_;
+  //   const bool canCollapse = (op == parent_op && IsUnique(node)) ||
+  //                            impl->size() == 1;
   //   const mat3x4 transform2 = canCollapse ? transform * node->transform_
   //                                         : la::identity;
-  //   std::vector<CsgLeafNode> positive_children, negative_children;
-  //   // for subtract, we pretend the operation is Add for our children.
-  //   auto op = node->op_ == OpType::Subtract ? OpType::Add : node->op_;
-  //   for (size_t i = 0; i < impl->children_.size(); i++) {
-  //     auto child = impl->children_[i];
-  //     // negative when it is the remaining operands for Subtract
-  //     auto dest = node->op_ == OpType::Subtract && i != 0 ?
-  //       negative_children : positive_children;
-  //     if (canCollapse) dest = destination;
+  //   Nodes positive_children, negative_children;
+  //   Nodes* pos_dest = canCollapse ? positive_dest : &positive_children;
+  //   Nodes* neg_dest = canCollapse ? negative_dest : &negative_children;
+  //   for (size_t i = 0; i < impl->size(); i++) {
+  //     auto child = (*impl)[i];
+  //     const bool negative = op == OpType::Subtract && i != 0;
+  //     Nodes *dest1 = negative ? neg_dest : pos_dest;
+  //     Nodes *dest2 = (op == OpType::Subtract && i == 0) ?
+  //                      neg_dest : nullptr;
   //     if (child->GetNodeType() == CsgNodeType::Leaf)
-  //       dest.push_back(child);
+  //       dest1.push_back(child);
   //     else
-  //       f(child, op, transform2, dest);
+  //       f(child, op, transform2, dest1, dest2);
   //   }
   //   if (canCollapse) return;
   //   if (node->op_ == OpType::Add)
-  //     impl->children_ = {BatchUnion(positive_children)};
+  //     *impl = {BatchUnion(positive_children)};
   //   else if (node->op_ == OpType::Intersect)
-  //     impl->children_ = {BatchBoolean(Intersect, positive_children)};
+  //     *impl = {BatchBoolean(Intersect, positive_children)};
   //   else // subtract
-  //     impl->children_ = { BatchUnion(positive_children) -
+  //     *impl = { BatchUnion(positive_children) -
   //                         BatchUnion(negative_children)};
   //   // node local transform
-  //   node->cache_ = impl->children_[0].Transform(node.transform);
+  //   node->cache_ = (*impl)[0].Transform(node.transform);
   //   // collapsed node transforms
   //   if (destination)
   //     destination->push_back(node->cache_->Transform(transform));
@@ -574,56 +644,56 @@ std::shared_ptr<CsgLeafNode> CsgOpNode::ToLeafNode() const {
     if (frame->finalize) {
       switch (frame->op_node->op_) {
         case OpType::Add:
-          impl->children_ = {BatchUnion(frame->positive_children)};
+          *impl = {BatchUnion(frame->positive_children)};
           break;
         case OpType::Intersect: {
-          impl->children_ = {
-              BatchBoolean(OpType::Intersect, frame->positive_children)};
+          *impl = {BatchBoolean(OpType::Intersect, frame->positive_children)};
           break;
         };
         case OpType::Subtract:
           if (frame->positive_children.empty()) {
             // nothing to subtract from, so the result is empty.
-            impl->children_ = {std::make_shared<CsgLeafNode>()};
+            *impl = {std::make_shared<CsgLeafNode>()};
           } else {
             auto positive = BatchUnion(frame->positive_children);
             if (frame->negative_children.empty()) {
               // nothing to subtract, result equal to the LHS.
-              impl->children_ = {frame->positive_children[0]};
+              *impl = {frame->positive_children[0]};
             } else {
-              Boolean3 boolean(*positive->GetImpl(),
-                               *BatchUnion(frame->negative_children)->GetImpl(),
-                               OpType::Subtract);
-              impl->children_ = {ImplToLeaf(boolean.Result(OpType::Subtract))};
+              auto negative = BatchUnion(frame->negative_children);
+              *impl = {SimpleBoolean(*positive->GetImpl(), *negative->GetImpl(),
+                                     OpType::Subtract)};
             }
           }
           break;
       }
       frame->op_node->cache_ = std::static_pointer_cast<CsgLeafNode>(
-          impl->children_[0]->Transform(frame->op_node->transform_));
-      if (frame->destination != nullptr)
-        frame->destination->push_back(std::static_pointer_cast<CsgLeafNode>(
+          (*impl)[0]->Transform(frame->op_node->transform_));
+      if (frame->positive_dest != nullptr)
+        frame->positive_dest->push_back(std::static_pointer_cast<CsgLeafNode>(
             frame->op_node->cache_->Transform(frame->transform)));
       stack.pop_back();
     } else {
-      auto add_children = [&stack](std::shared_ptr<CsgNode> &node, OpType op,
-                                   mat3x4 transform, auto *destination) {
-        if (node->GetNodeType() == CsgNodeType::Leaf)
-          destination->push_back(std::static_pointer_cast<CsgLeafNode>(
-              node->Transform(transform)));
-        else
-          stack.push_back(std::make_shared<CsgStackFrame>(
-              false, op, transform, destination,
-              std::static_pointer_cast<const CsgOpNode>(node)));
-      };
+      auto add_children =
+          [&stack](std::shared_ptr<CsgNode> &node, OpType op, mat3x4 transform,
+                   CsgStackFrame::Nodes *dest1, CsgStackFrame::Nodes *dest2) {
+            if (node->GetNodeType() == CsgNodeType::Leaf)
+              dest1->push_back(std::static_pointer_cast<CsgLeafNode>(
+                  node->Transform(transform)));
+            else
+              stack.push_back(std::make_shared<CsgStackFrame>(
+                  false, op, transform, dest1, dest2,
+                  std::static_pointer_cast<const CsgOpNode>(node)));
+          };
       // op_node use_count == 2 because it is both inside one CsgOpNode
       // and in our stack.
       // if there is only one child, we can also collapse.
-      const bool canCollapse = frame->destination != nullptr &&
-                               ((frame->op_node->op_ == frame->parent_op &&
-                                 frame->op_node.use_count() <= 2 &&
-                                 frame->op_node->impl_.UseCount() == 1) ||
-                                impl->children_.size() == 1);
+      const OpType op = frame->op_node->op_;
+      const bool canCollapse =
+          frame->positive_dest != nullptr &&
+          ((op == frame->parent_op && frame->op_node.use_count() <= 2 &&
+            frame->op_node->impl_.UseCount() == 1) ||
+           impl->size() == 1);
       if (canCollapse)
         stack.pop_back();
       else
@@ -632,14 +702,17 @@ std::shared_ptr<CsgLeafNode> CsgOpNode::ToLeafNode() const {
       const mat3x4 transform =
           canCollapse ? (frame->transform * Mat4(frame->op_node->transform_))
                       : la::identity;
-      OpType op = frame->op_node->op_ == OpType::Subtract ? OpType::Add
-                                                          : frame->op_node->op_;
-      for (size_t i = 0; i < impl->children_.size(); i++) {
-        auto dest = canCollapse ? frame->destination
-                    : (frame->op_node->op_ == OpType::Subtract && i != 0)
-                        ? &frame->negative_children
-                        : &frame->positive_children;
-        add_children(impl->children_[i], op, transform, dest);
+      CsgStackFrame::Nodes *pos_dest =
+          canCollapse ? frame->positive_dest : &frame->positive_children;
+      CsgStackFrame::Nodes *neg_dest =
+          canCollapse ? frame->negative_dest : &frame->negative_children;
+      for (size_t i = 0; i < impl->size(); i++) {
+        const bool negative = op == OpType::Subtract && i != 0;
+        CsgStackFrame::Nodes *dest1 = negative ? neg_dest : pos_dest;
+        CsgStackFrame::Nodes *dest2 =
+            (op == OpType::Subtract && i == 0) ? neg_dest : nullptr;
+        add_children((*impl)[i], negative ? OpType::Add : op, transform, dest1,
+                     dest2);
       }
     }
   }
diff --git a/thirdparty/manifold/src/csg_tree.h b/thirdparty/manifold/src/csg_tree.h
index 5c24dc5a2050..fb664df00c08 100644
--- a/thirdparty/manifold/src/csg_tree.h
+++ b/thirdparty/manifold/src/csg_tree.h
@@ -74,12 +74,11 @@ class CsgOpNode final : public CsgNode {
 
   CsgNodeType GetNodeType() const override;
 
+  ~CsgOpNode();
+
  private:
-  struct Impl {
-    std::vector<std::shared_ptr<CsgNode>> children_;
-    bool forcedToLeafNodes_ = false;
-  };
-  mutable ConcurrentSharedPtr<Impl> impl_ = ConcurrentSharedPtr<Impl>(Impl{});
+  mutable ConcurrentSharedPtr<std::vector<std::shared_ptr<CsgNode>>> impl_ =
+      ConcurrentSharedPtr<std::vector<std::shared_ptr<CsgNode>>>({});
   OpType op_;
   mat3x4 transform_ = la::identity;
   // the following fields are for lazy evaluation, so they are mutable
diff --git a/thirdparty/manifold/src/impl.cpp b/thirdparty/manifold/src/impl.cpp
index 59cc0293d6da..afa45dff3020 100644
--- a/thirdparty/manifold/src/impl.cpp
+++ b/thirdparty/manifold/src/impl.cpp
@@ -17,12 +17,19 @@
 #include <algorithm>
 #include <atomic>
 #include <map>
+#include <optional>
 
 #include "./hashtable.h"
 #include "./mesh_fixes.h"
 #include "./parallel.h"
 #include "./svd.h"
 
+#ifdef MANIFOLD_EXPORT
+#include <string.h>
+
+#include <iostream>
+#endif
+
 namespace {
 using namespace manifold;
 
@@ -269,38 +276,21 @@ Manifold::Impl::Impl(Shape shape, const mat3x4 m) {
 
 void Manifold::Impl::RemoveUnreferencedVerts() {
   ZoneScoped;
-  Vec<int> vertOld2New(NumVert(), 0);
-  auto policy = autoPolicy(NumVert(), 1e5);
-  for_each(policy, halfedge_.cbegin(), halfedge_.cend(),
-           [&vertOld2New](Halfedge h) {
-             reinterpret_cast<std::atomic<int>*>(&vertOld2New[h.startVert])
-                 ->store(1, std::memory_order_relaxed);
-           });
-
-  const Vec<vec3> oldVertPos = vertPos_;
-
-  Vec<size_t> tmpBuffer(oldVertPos.size());
-  auto vertIdIter = TransformIterator(countAt(0_uz), [&vertOld2New](size_t i) {
-    if (vertOld2New[i] > 0) return i;
-    return std::numeric_limits<size_t>::max();
+  const int numVert = NumVert();
+  Vec<int> keep(numVert, 0);
+  auto policy = autoPolicy(numVert, 1e5);
+  for_each(policy, halfedge_.cbegin(), halfedge_.cend(), [&keep](Halfedge h) {
+    if (h.startVert >= 0) {
+      reinterpret_cast<std::atomic<int>*>(&keep[h.startVert])
+          ->store(1, std::memory_order_relaxed);
+    }
   });
 
-  auto next =
-      copy_if(vertIdIter, vertIdIter + tmpBuffer.size(), tmpBuffer.begin(),
-              [](size_t v) { return v != std::numeric_limits<size_t>::max(); });
-  if (next == tmpBuffer.end()) return;
-
-  gather(tmpBuffer.begin(), next, oldVertPos.begin(), vertPos_.begin());
-
-  vertPos_.resize(std::distance(tmpBuffer.begin(), next));
-
-  exclusive_scan(vertOld2New.begin(), vertOld2New.end(), vertOld2New.begin());
-
-  for_each(policy, halfedge_.begin(), halfedge_.end(),
-           [&vertOld2New](Halfedge& h) {
-             h.startVert = vertOld2New[h.startVert];
-             h.endVert = vertOld2New[h.endVert];
-           });
+  for_each_n(policy, countAt(0), numVert, [&keep, this](int v) {
+    if (keep[v] == 0) {
+      vertPos_[v] = vec3(NAN);
+    }
+  });
 }
 
 void Manifold::Impl::InitializeOriginal(bool keepFaceID) {
@@ -429,7 +419,8 @@ void Manifold::Impl::CreateHalfedges(const Vec<ivec3>& triVerts) {
     return edge[a] < edge[b];
   });
 
-  // Mark opposed triangles for removal
+  // Mark opposed triangles for removal - this may strand unreferenced verts
+  // which are removed later by RemoveUnreferencedVerts() and Finish().
   const int numEdge = numHalfedge / 2;
   for (int i = 0; i < numEdge; ++i) {
     const int pair0 = ids[i];
@@ -441,6 +432,8 @@ void Manifold::Impl::CreateHalfedges(const Vec<ivec3>& triVerts) {
       if (h0.startVert != h1.endVert || h0.endVert != h1.startVert) break;
       if (halfedge_[NextHalfedge(pair0)].endVert ==
           halfedge_[NextHalfedge(pair1)].endVert) {
+        h0 = {-1, -1, -1};
+        h1 = {-1, -1, -1};
         // Reorder so that remaining edges pair up
         if (k != i + numEdge) std::swap(ids[i + numEdge], ids[k]);
         break;
@@ -456,12 +449,11 @@ void Manifold::Impl::CreateHalfedges(const Vec<ivec3>& triVerts) {
   for_each_n(policy, countAt(0), numEdge, [this, &ids, numEdge](int i) {
     const int pair0 = ids[i];
     const int pair1 = ids[i + numEdge];
-    halfedge_[pair0].pairedHalfedge = pair1;
-    halfedge_[pair1].pairedHalfedge = pair0;
+    if (halfedge_[pair0].startVert >= 0) {
+      halfedge_[pair0].pairedHalfedge = pair1;
+      halfedge_[pair1].pairedHalfedge = pair0;
+    }
   });
-
-  // When opposed triangles are removed, they may strand unreferenced verts.
-  RemoveUnreferencedVerts();
 }
 
 /**
@@ -686,4 +678,90 @@ SparseIndices Manifold::Impl::VertexCollisionsZ(VecView<const vec3> vertsIn,
     return collider_.Collisions<false, false>(vertsIn);
 }
 
+#ifdef MANIFOLD_DEBUG
+std::ostream& operator<<(std::ostream& stream, const Manifold::Impl& impl) {
+  stream << std::setprecision(17);  // for double precision
+  stream << "# ======= begin mesh ======" << std::endl;
+  stream << "# tolerance = " << impl.tolerance_ << std::endl;
+  stream << "# epsilon = " << impl.epsilon_ << std::endl;
+  // TODO: Mesh relation, vertex normal and face normal
+  for (const vec3& v : impl.vertPos_)
+    stream << "v " << v.x << " " << v.y << " " << v.z << std::endl;
+  std::vector<ivec3> triangles;
+  triangles.reserve(impl.halfedge_.size() / 3);
+  for (size_t i = 0; i < impl.halfedge_.size(); i += 3)
+    triangles.emplace_back(impl.halfedge_[i].startVert + 1,
+                           impl.halfedge_[i + 1].startVert + 1,
+                           impl.halfedge_[i + 2].startVert + 1);
+  sort(triangles.begin(), triangles.end());
+  for (const auto& tri : triangles)
+    stream << "f " << tri.x << " " << tri.y << " " << tri.z << std::endl;
+  stream << "# ======== end mesh =======" << std::endl;
+  return stream;
+}
+#endif
+
+#ifdef MANIFOLD_EXPORT
+Manifold Manifold::ImportMeshGL64(std::istream& stream) {
+  MeshGL64 mesh;
+  std::optional<double> epsilon;
+  stream.precision(17);
+  while (true) {
+    char c = stream.get();
+    if (stream.eof()) break;
+    switch (c) {
+      case '#': {
+        char c = stream.get();
+        if (c == ' ') {
+          constexpr int SIZE = 10;
+          std::array<char, SIZE> tmp;
+          stream.get(tmp.data(), SIZE, '\n');
+          if (strncmp(tmp.data(), "tolerance", SIZE) == 0) {
+            // skip 3 letters
+            for (int i : {0, 1, 2}) stream.get();
+            stream >> mesh.tolerance;
+          } else if (strncmp(tmp.data(), "epsilon =", SIZE) == 0) {
+            double tmp;
+            stream >> tmp;
+            epsilon = {tmp};
+          } else {
+            // add it back because it is not what we want
+            int end = 0;
+            while (tmp[end] != 0 && end < SIZE) end++;
+            while (--end > -1) stream.putback(tmp[end]);
+          }
+          c = stream.get();
+        }
+        // just skip the remaining comment
+        while (c != '\n' && !stream.eof()) {
+          c = stream.get();
+        }
+        break;
+      }
+      case 'v':
+        for (int i : {0, 1, 2}) {
+          double x;
+          stream >> x;
+          mesh.vertProperties.push_back(x);
+        }
+        break;
+      case 'f':
+        for (int i : {0, 1, 2}) {
+          uint64_t x;
+          stream >> x;
+          mesh.triVerts.push_back(x - 1);
+        }
+        break;
+      case '\n':
+        break;
+      default:
+        DEBUG_ASSERT(false, userErr, "unexpected character in MeshGL64 import");
+    }
+  }
+  auto m = std::make_shared<Manifold::Impl>(mesh);
+  if (epsilon) m->SetEpsilon(*epsilon);
+  return Manifold(m);
+}
+#endif
+
 }  // namespace manifold
diff --git a/thirdparty/manifold/src/impl.h b/thirdparty/manifold/src/impl.h
index 8f7d700098f1..d37a9a7d27e9 100644
--- a/thirdparty/manifold/src/impl.h
+++ b/thirdparty/manifold/src/impl.h
@@ -198,10 +198,6 @@ struct Manifold::Impl {
     }
 
     CalculateBBox();
-    if (!IsFinite()) {
-      MarkFailure(Error::NonFiniteVertex);
-      return;
-    }
     SetEpsilon(-1, std::is_same<Precision, float>::value);
 
     SplitPinchedVerts();
@@ -215,8 +211,14 @@ struct Manifold::Impl {
     CreateFaces();
 
     SimplifyTopology();
+    RemoveUnreferencedVerts();
     Finish();
 
+    if (!IsFinite()) {
+      MarkFailure(Error::NonFiniteVertex);
+      return;
+    }
+
     // A Manifold created from an input mesh is never an original - the input is
     // the original.
     meshRelation_.originalID = -1;
@@ -271,7 +273,7 @@ struct Manifold::Impl {
                           : meshRelation_.properties.size() / NumProp();
   }
 
-  // properties.cu
+  // properties.cpp
   enum class Property { Volume, SurfaceArea };
   double GetProperty(Property prop) const;
   void CalculateCurvature(int gaussianIdx, int meanIdx);
@@ -281,11 +283,12 @@ struct Manifold::Impl {
   void SetEpsilon(double minEpsilon = -1, bool useSingle = false);
   bool IsManifold() const;
   bool Is2Manifold() const;
+  bool IsSelfIntersecting() const;
   bool MatchesTriNormals() const;
   int NumDegenerateTris() const;
   double MinGap(const Impl& other, double searchLength) const;
 
-  // sort.cu
+  // sort.cpp
   void Finish();
   void SortVerts();
   void ReindexVerts(const Vec<int>& vertNew2Old, size_t numOldVert);
@@ -295,7 +298,7 @@ struct Manifold::Impl {
   void GatherFaces(const Vec<int>& faceNew2Old);
   void GatherFaces(const Impl& old, const Vec<int>& faceNew2Old);
 
-  // face_op.cu
+  // face_op.cpp
   void Face2Tri(const Vec<int>& faceEdge, const Vec<TriRef>& halfedgeRef);
   PolygonsIdx Face2Polygons(VecView<Halfedge>::IterC start,
                             VecView<Halfedge>::IterC end,
@@ -303,7 +306,7 @@ struct Manifold::Impl {
   Polygons Slice(double height) const;
   Polygons Project() const;
 
-  // edge_op.cu
+  // edge_op.cpp
   void CleanupTopology();
   void SimplifyTopology();
   void DedupeEdge(int edge);
@@ -349,4 +352,9 @@ struct Manifold::Impl {
   // quickhull.cpp
   void Hull(VecView<vec3> vertPos);
 };
+
+#ifdef MANIFOLD_DEBUG
+extern std::mutex dump_lock;
+std::ostream& operator<<(std::ostream& stream, const Manifold::Impl& impl);
+#endif
 }  // namespace manifold
diff --git a/thirdparty/manifold/src/properties.cpp b/thirdparty/manifold/src/properties.cpp
index 911133776ab3..5b6403a6c15d 100644
--- a/thirdparty/manifold/src/properties.cpp
+++ b/thirdparty/manifold/src/properties.cpp
@@ -202,6 +202,66 @@ bool Manifold::Impl::Is2Manifold() const {
       });
 }
 
+#ifdef MANIFOLD_DEBUG
+std::mutex dump_lock;
+#endif
+
+/**
+ * Returns true if this manifold is self-intersecting.
+ * Note that this is not checking for epsilon-validity.
+ */
+bool Manifold::Impl::IsSelfIntersecting() const {
+  const double epsilonSq = epsilon_ * epsilon_;
+  Vec<Box> faceBox;
+  Vec<uint32_t> faceMorton;
+  GetFaceBoxMorton(faceBox, faceMorton);
+  SparseIndices collisions = collider_.Collisions<true>(faceBox.cview());
+
+  const bool verbose = ManifoldParams().verbose;
+  return !all_of(countAt(0), countAt(collisions.size()), [&](size_t i) {
+    size_t x = collisions.Get(i, false);
+    size_t y = collisions.Get(i, true);
+    std::array<vec3, 3> tri_x, tri_y;
+    for (int i : {0, 1, 2}) {
+      tri_x[i] = vertPos_[halfedge_[3 * x + i].startVert];
+      tri_y[i] = vertPos_[halfedge_[3 * y + i].startVert];
+    }
+    // if triangles x and y share a vertex, return true to skip the
+    // check. we relax the sharing criteria a bit to allow for at most
+    // distance epsilon squared
+    for (int i : {0, 1, 2})
+      for (int j : {0, 1, 2})
+        if (distance2(tri_x[i], tri_y[j]) <= epsilonSq) return true;
+
+    if (DistanceTriangleTriangleSquared(tri_x, tri_y) == 0.0) {
+      // try to move the triangles around the normal of the other face
+      std::array<vec3, 3> tmp_x, tmp_y;
+      for (int i : {0, 1, 2}) tmp_x[i] = tri_x[i] + epsilon_ * faceNormal_[y];
+      if (DistanceTriangleTriangleSquared(tmp_x, tri_y) > 0.0) return true;
+      for (int i : {0, 1, 2}) tmp_x[i] = tri_x[i] - epsilon_ * faceNormal_[y];
+      if (DistanceTriangleTriangleSquared(tmp_x, tri_y) > 0.0) return true;
+      for (int i : {0, 1, 2}) tmp_y[i] = tri_y[i] + epsilon_ * faceNormal_[x];
+      if (DistanceTriangleTriangleSquared(tri_x, tmp_y) > 0.0) return true;
+      for (int i : {0, 1, 2}) tmp_y[i] = tri_y[i] - epsilon_ * faceNormal_[x];
+      if (DistanceTriangleTriangleSquared(tri_x, tmp_y) > 0.0) return true;
+
+#ifdef MANIFOLD_DEBUG
+      if (verbose) {
+        dump_lock.lock();
+        std::cout << "intersecting:" << std::endl;
+        for (int i : {0, 1, 2}) std::cout << tri_x[i] << " ";
+        std::cout << std::endl;
+        for (int i : {0, 1, 2}) std::cout << tri_y[i] << " ";
+        std::cout << std::endl;
+        dump_lock.unlock();
+      }
+#endif
+      return false;
+    }
+    return true;
+  });
+}
+
 /**
  * Returns true if all triangles are CCW relative to their triNormals_.
  */
diff --git a/thirdparty/manifold/src/sdf.cpp b/thirdparty/manifold/src/sdf.cpp
index 64dd7094d7f4..ad95b9b6fe40 100644
--- a/thirdparty/manifold/src/sdf.cpp
+++ b/thirdparty/manifold/src/sdf.cpp
@@ -526,6 +526,7 @@ Manifold Manifold::LevelSet(std::function<double(vec3)> sdf, Box bounds,
 
   pImpl_->CreateHalfedges(triVerts);
   pImpl_->CleanupTopology();
+  pImpl_->RemoveUnreferencedVerts();
   pImpl_->Finish();
   pImpl_->InitializeOriginal();
   return Manifold(pImpl_);