diff --git a/.github/workflows/manifold.yml b/.github/workflows/manifold.yml
index 90fd1bf0a..358a1c93e 100644
--- a/.github/workflows/manifold.yml
+++ b/.github/workflows/manifold.yml
@@ -17,6 +17,7 @@ jobs:
         cuda_support: [ON, OFF]
         parallel_backend: [NONE, OMP, TBB]
     runs-on: ubuntu-20.04
+    if: github.event.pull_request.draft == false
     container:
       image: docker://nvidia/cuda:11.6.0-devel-ubuntu20.04
     steps:
@@ -47,6 +48,7 @@ jobs:
 
   build_wasm:
     runs-on: ubuntu-20.04
+    if: github.event.pull_request.draft == false
     steps:
     - name: Install dependencies
       run: |
@@ -83,6 +85,7 @@ jobs:
         cuda_support: [OFF]
       max-parallel: 1
     runs-on: windows-2019
+    if: github.event.pull_request.draft == false
     steps:
     - uses: actions/checkout@v3
       with:
@@ -126,6 +129,7 @@ jobs:
       matrix:
         variant: [none, omp, tbb, none-cuda, omp-cuda, tbb-cuda, js]
     runs-on: ubuntu-latest
+    if: github.event.pull_request.draft == false
     steps:
     - uses: actions/checkout@v2.4.0
     - uses: cachix/install-nix-action@v15
diff --git a/clang-format.sh b/clang-format.sh
index 9b512a6a4..b2650d256 100755
--- a/clang-format.sh
+++ b/clang-format.sh
@@ -1,4 +1,9 @@
 #!/usr/bin/env -S bash -e
 for f in {collider,manifold,meshIO,polygon,samples,utilities}/**/*.{h,cpp} {test,tools}/*.cpp; do
-  clang-format --dry-run --Werror $f
+  clang-format --dry-run --Werror $f &
+  pids[${i}]=$!
+done
+
+for pid in ${pids[*]}; do
+    wait $pid
 done
diff --git a/manifold/CMakeLists.txt b/manifold/CMakeLists.txt
index 482166f1e..e41527499 100644
--- a/manifold/CMakeLists.txt
+++ b/manifold/CMakeLists.txt
@@ -2,7 +2,7 @@ project (manifold)
 
 add_library(${PROJECT_NAME} src/manifold.cpp src/constructors.cpp src/impl.cpp
     src/properties.cpp src/sort.cpp src/edge_op.cpp src/face_op.cpp
-    src/smoothing.cpp src/boolean3.cpp src/boolean_result.cpp
+    src/smoothing.cpp src/boolean3.cpp src/boolean_result.cpp src/csg_tree.cpp
 )
 
 if(MANIFOLD_USE_CUDA)
@@ -10,7 +10,7 @@ if(MANIFOLD_USE_CUDA)
         src/manifold.cpp src/constructors.cpp src/impl.cpp
         src/properties.cpp src/sort.cpp src/edge_op.cpp src/face_op.cpp
         src/smoothing.cpp src/boolean3.cpp src/boolean_result.cpp
-        src/manifold_set.cpp
+        src/manifold_set.cpp src/csg_tree.cpp
 
         PROPERTIES LANGUAGE CUDA
     )
diff --git a/manifold/include/manifold.h b/manifold/include/manifold.h
index b101c34a6..d73d19503 100644
--- a/manifold/include/manifold.h
+++ b/manifold/include/manifold.h
@@ -20,6 +20,9 @@
 
 namespace manifold {
 
+class CsgNode;
+class CsgLeafNode;
+
 /** @defgroup Core
  *  @brief The central classes of the library
  *  @{
@@ -104,20 +107,19 @@ class Manifold {
   ///@{
   MeshRelation GetMeshRelation() const;
   std::vector<int> GetMeshIDs() const;
-  int SetAsOriginal();
+  Manifold AsOriginal() const;
   ///@}
 
   /** @name Modification
-   *  Change this manifold in-place.
    */
   ///@{
-  Manifold& Translate(glm::vec3);
-  Manifold& Scale(glm::vec3);
-  Manifold& Rotate(float xDegrees, float yDegrees = 0.0f,
-                   float zDegrees = 0.0f);
-  Manifold& Transform(const glm::mat4x3&);
-  Manifold& Warp(std::function<void(glm::vec3&)>);
-  Manifold& Refine(int);
+  Manifold Translate(glm::vec3) const;
+  Manifold Scale(glm::vec3) const;
+  Manifold Rotate(float xDegrees, float yDegrees = 0.0f,
+                  float zDegrees = 0.0f) const;
+  Manifold Transform(const glm::mat4x3&) const;
+  Manifold Warp(std::function<void(glm::vec3&)>) const;
+  Manifold Refine(int) const;
   // Manifold RefineToLength(float);
   // Manifold RefineToPrecision(float);
   ///@}
@@ -157,7 +159,13 @@ class Manifold {
   struct Impl;
 
  private:
-  std::unique_ptr<Impl> pImpl_;
+  Manifold(std::shared_ptr<CsgNode> pNode_);
+  Manifold(std::shared_ptr<Impl> pImpl_);
+
+  mutable std::shared_ptr<CsgNode> pNode_;
+
+  CsgLeafNode& GetCsgLeafNode() const;
+
   static int circularSegments_;
   static float circularAngle_;
   static float circularEdgeLength_;
diff --git a/manifold/src/constructors.cpp b/manifold/src/constructors.cpp
index 8fdef1099..7f4b07353 100644
--- a/manifold/src/constructors.cpp
+++ b/manifold/src/constructors.cpp
@@ -14,6 +14,7 @@
 
 #include <thrust/sequence.h>
 
+#include "csg_tree.h"
 #include "graph.h"
 #include "impl.h"
 #include "par.h"
@@ -32,30 +33,6 @@ struct ToSphere {
   }
 };
 
-struct UpdateTriBary {
-  const int nextBary;
-
-  __host__ __device__ BaryRef operator()(BaryRef ref) {
-    for (int i : {0, 1, 2})
-      if (ref.vertBary[i] >= 0) ref.vertBary[i] += nextBary;
-    return ref;
-  }
-};
-
-struct UpdateHalfedge {
-  const int nextVert;
-  const int nextEdge;
-  const int nextFace;
-
-  __host__ __device__ Halfedge operator()(Halfedge edge) {
-    edge.startVert += nextVert;
-    edge.endVert += nextVert;
-    edge.pairedHalfedge += nextEdge;
-    edge.face += nextFace;
-    return edge;
-  }
-};
-
 struct Equals {
   int val;
   __host__ __device__ bool operator()(int x) { return x == val; }
@@ -108,9 +85,9 @@ Manifold Manifold::Smooth(const Mesh& mesh,
       "when supplying tangents, the normal constructor should be used "
       "rather than Smooth().");
 
-  Manifold manifold(mesh);
-  manifold.pImpl_->CreateTangents(sharpenedEdges);
-  return manifold;
+  std::shared_ptr<Impl> impl = std::make_shared<Impl>(mesh);
+  impl->CreateTangents(sharpenedEdges);
+  return Manifold(impl);
 }
 
 /**
@@ -118,9 +95,7 @@ Manifold Manifold::Smooth(const Mesh& mesh,
  * and the rest at similarly symmetric points.
  */
 Manifold Manifold::Tetrahedron() {
-  Manifold tetrahedron;
-  tetrahedron.pImpl_ = std::make_unique<Impl>(Impl::Shape::TETRAHEDRON);
-  return tetrahedron;
+  return Manifold(std::make_shared<Impl>(Impl::Shape::TETRAHEDRON));
 }
 
 /**
@@ -131,10 +106,9 @@ Manifold Manifold::Tetrahedron() {
  * @param center Set to true to shift the center to the origin.
  */
 Manifold Manifold::Cube(glm::vec3 size, bool center) {
-  Manifold cube;
-  cube.pImpl_ = std::make_unique<Impl>(Impl::Shape::CUBE);
-  cube.Scale(size);
-  if (center) cube.Translate(-size / 2.0f);
+  auto cube = Manifold(std::make_shared<Impl>(Impl::Shape::CUBE));
+  cube = cube.Scale(size);
+  if (center) cube = cube.Translate(-size / 2.0f);
   return cube;
 }
 
@@ -164,7 +138,8 @@ Manifold Manifold::Cylinder(float height, float radiusLow, float radiusHigh,
   }
   Manifold cylinder =
       Manifold::Extrude(circle, height, 0, 0.0f, glm::vec2(scale));
-  if (center) cylinder.Translate(glm::vec3(0.0f, 0.0f, -height / 2.0f));
+  if (center)
+    cylinder = cylinder.Translate(glm::vec3(0.0f, 0.0f, -height / 2.0f));
   return cylinder;
 }
 
@@ -181,15 +156,14 @@ Manifold Manifold::Cylinder(float height, float radiusLow, float radiusHigh,
 Manifold Manifold::Sphere(float radius, int circularSegments) {
   int n = circularSegments > 0 ? (circularSegments + 3) / 4
                                : GetCircularSegments(radius) / 4;
-  Manifold sphere;
-  sphere.pImpl_ = std::make_unique<Impl>(Impl::Shape::OCTAHEDRON);
-  sphere.pImpl_->Subdivide(n);
-  for_each_n(autoPolicy(sphere.NumVert()), sphere.pImpl_->vertPos_.begin(),
-             sphere.NumVert(), ToSphere({radius}));
-  sphere.pImpl_->Finish();
+  auto pImpl_ = std::make_shared<Impl>(Impl::Shape::OCTAHEDRON);
+  pImpl_->Subdivide(n);
+  for_each_n(autoPolicy(pImpl_->NumVert()), pImpl_->vertPos_.begin(),
+             pImpl_->NumVert(), ToSphere({radius}));
+  pImpl_->Finish();
   // Ignore preceding octahedron.
-  sphere.pImpl_->ReinitializeReference(Impl::meshIDCounter_.fetch_add(1));
-  return sphere;
+  pImpl_->ReinitializeReference(Impl::meshIDCounter_.fetch_add(1));
+  return Manifold(pImpl_);
 }
 
 /**
@@ -211,9 +185,9 @@ Manifold Manifold::Extrude(Polygons crossSection, float height, int nDivisions,
                            float twistDegrees, glm::vec2 scaleTop) {
   ALWAYS_ASSERT(scaleTop.x >= 0 && scaleTop.y >= 0, userErr,
                 "scale values cannot be negative");
-  Manifold extrusion;
+  auto pImpl_ = std::make_shared<Impl>();
   ++nDivisions;
-  auto& vertPos = extrusion.pImpl_->vertPos_;
+  auto& vertPos = pImpl_->vertPos_;
   VecDH<glm::ivec3> triVertsDH;
   auto& triVerts = triVertsDH;
   int nCrossSection = 0;
@@ -263,10 +237,10 @@ Manifold Manifold::Extrude(Polygons crossSection, float height, int nDivisions,
     if (!isCone) triVerts.push_back(tri + nCrossSection * nDivisions);
   }
 
-  extrusion.pImpl_->CreateHalfedges(triVertsDH);
-  extrusion.pImpl_->Finish();
-  extrusion.pImpl_->InitializeNewReference();
-  return extrusion;
+  pImpl_->CreateHalfedges(triVertsDH);
+  pImpl_->Finish();
+  pImpl_->InitializeNewReference();
+  return Manifold(pImpl_);
 }
 
 /**
@@ -289,8 +263,8 @@ Manifold Manifold::Revolve(const Polygons& crossSection, int circularSegments) {
   }
   int nDivisions =
       circularSegments > 2 ? circularSegments : GetCircularSegments(radius);
-  Manifold revoloid;
-  auto& vertPos = revoloid.pImpl_->vertPos_;
+  auto pImpl_ = std::make_shared<Impl>();
+  auto& vertPos = pImpl_->vertPos_;
   VecDH<glm::ivec3> triVertsDH;
   auto& triVerts = triVertsDH;
   float dPhi = 360.0f / nDivisions;
@@ -362,10 +336,10 @@ Manifold Manifold::Revolve(const Polygons& crossSection, int circularSegments) {
     }
   }
 
-  revoloid.pImpl_->CreateHalfedges(triVertsDH);
-  revoloid.pImpl_->Finish();
-  revoloid.pImpl_->InitializeNewReference();
-  return revoloid;
+  pImpl_->CreateHalfedges(triVertsDH);
+  pImpl_->Finish();
+  pImpl_->InitializeNewReference();
+  return Manifold(pImpl_);
 }
 
 /**
@@ -376,73 +350,11 @@ Manifold Manifold::Revolve(const Polygons& crossSection, int circularSegments) {
  * @param manifolds A vector of Manifolds to lazy-union together.
  */
 Manifold Manifold::Compose(const std::vector<Manifold>& manifolds) {
-  int numVert = 0;
-  int numEdge = 0;
-  int numTri = 0;
-  int numBary = 0;
-  for (const Manifold& manifold : manifolds) {
-    numVert += manifold.NumVert();
-    numEdge += manifold.NumEdge();
-    numTri += manifold.NumTri();
-    numBary += manifold.pImpl_->meshRelation_.barycentric.size();
+  std::vector<std::shared_ptr<CsgLeafNode>> children;
+  for (const auto& manifold : manifolds) {
+    children.push_back(manifold.pNode_->ToLeafNode());
   }
-
-  Manifold out;
-  Impl& combined = *(out.pImpl_);
-  combined.vertPos_.resize(numVert);
-  combined.halfedge_.resize(2 * numEdge);
-  combined.faceNormal_.resize(numTri);
-  combined.halfedgeTangent_.resize(2 * numEdge);
-  combined.meshRelation_.barycentric.resize(numBary);
-  combined.meshRelation_.triBary.resize(numTri);
-  auto policy = autoPolicy(numTri);
-
-  int nextVert = 0;
-  int nextEdge = 0;
-  int nextTri = 0;
-  int nextBary = 0;
-  for (const Manifold& manifold : manifolds) {
-    const Impl& impl = *(manifold.pImpl_);
-    impl.ApplyTransform();
-
-    copy(policy, impl.vertPos_.begin(), impl.vertPos_.end(),
-         combined.vertPos_.begin() + nextVert);
-    copy(policy, impl.faceNormal_.begin(), impl.faceNormal_.end(),
-         combined.faceNormal_.begin() + nextTri);
-    copy(policy, impl.halfedgeTangent_.begin(), impl.halfedgeTangent_.end(),
-         combined.halfedgeTangent_.begin() + nextEdge);
-    copy(policy, impl.meshRelation_.barycentric.begin(),
-         impl.meshRelation_.barycentric.end(),
-         combined.meshRelation_.barycentric.begin() + nextBary);
-    transform(policy, impl.meshRelation_.triBary.begin(),
-              impl.meshRelation_.triBary.end(),
-              combined.meshRelation_.triBary.begin() + nextTri,
-              UpdateTriBary({nextBary}));
-    transform(policy, impl.halfedge_.begin(), impl.halfedge_.end(),
-              combined.halfedge_.begin() + nextEdge,
-              UpdateHalfedge({nextVert, nextEdge, nextTri}));
-
-    // Assign new IDs to triangles added in this iteration, to differentiate
-    // triangles coming from different manifolds.
-    // See the end of `boolean_result.cpp` for details.
-    VecDH<int> meshIDs;
-    VecDH<int> original;
-    for (auto& entry : impl.meshRelation_.originalID) {
-      meshIDs.push_back(entry.first);
-      original.push_back(entry.second);
-    }
-    int meshIDStart = combined.meshRelation_.originalID.size();
-    combined.UpdateMeshIDs(meshIDs, original, nextTri,
-                           impl.meshRelation_.triBary.size(), meshIDStart);
-
-    nextVert += manifold.NumVert();
-    nextEdge += 2 * manifold.NumEdge();
-    nextTri += manifold.NumTri();
-    nextBary += impl.meshRelation_.barycentric.size();
-  }
-
-  combined.Finish();
-  return out;
+  return Manifold(std::make_shared<Impl>(CsgLeafNode::Compose(children)));
 }
 
 /**
@@ -452,6 +364,7 @@ Manifold Manifold::Compose(const std::vector<Manifold>& manifolds) {
  */
 std::vector<Manifold> Manifold::Decompose() const {
   Graph graph;
+  auto pImpl_ = GetCsgLeafNode().GetImpl();
   for (int i = 0; i < NumVert(); ++i) {
     graph.add_nodes(i);
   }
@@ -476,20 +389,22 @@ std::vector<Manifold> Manifold::Decompose() const {
     original.push_back(entry.second);
   }
 
-  std::vector<Manifold> meshes(numLabel);
+  std::vector<Manifold> meshes;
   for (int i = 0; i < numLabel; ++i) {
-    meshes[i].pImpl_->vertPos_.resize(NumVert());
+    auto impl = std::make_shared<Impl>();
+    // inherit original object's precision
+    impl->precision_ = pImpl_->precision_;
+    impl->vertPos_.resize(NumVert());
     VecDH<int> vertNew2Old(NumVert());
     auto policy = autoPolicy(NumVert());
-    auto start = zip(meshes[i].pImpl_->vertPos_.begin(), vertNew2Old.begin());
+    auto start = zip(impl->vertPos_.begin(), vertNew2Old.begin());
     int nVert =
         copy_if<decltype(start)>(
             policy, zip(pImpl_->vertPos_.begin(), countAt(0)),
             zip(pImpl_->vertPos_.end(), countAt(NumVert())), vertLabel.begin(),
-            zip(meshes[i].pImpl_->vertPos_.begin(), vertNew2Old.begin()),
-            Equals({i})) -
+            zip(impl->vertPos_.begin(), vertNew2Old.begin()), Equals({i})) -
         start;
-    meshes[i].pImpl_->vertPos_.resize(nVert);
+    impl->vertPos_.resize(nVert);
 
     VecDH<int> faceNew2Old(NumTri());
     sequence(policy, faceNew2Old.begin(), faceNew2Old.end());
@@ -501,16 +416,15 @@ std::vector<Manifold> Manifold::Decompose() const {
         faceNew2Old.begin();
     faceNew2Old.resize(nFace);
 
-    meshes[i].pImpl_->GatherFaces(*pImpl_, faceNew2Old);
-    meshes[i].pImpl_->ReindexVerts(vertNew2Old, pImpl_->NumVert());
+    impl->GatherFaces(*pImpl_, faceNew2Old);
+    impl->ReindexVerts(vertNew2Old, pImpl_->NumVert());
 
-    meshes[i].pImpl_->Finish();
+    impl->Finish();
 
     // meshIDs and original will only be sorted after successful updates, so we
     // can keep using the old one.
-    meshes[i].pImpl_->UpdateMeshIDs(meshIDs, original);
-
-    meshes[i].pImpl_->transform_ = pImpl_->transform_;
+    impl->UpdateMeshIDs(meshIDs, original);
+    meshes.push_back(Manifold(impl));
   }
   return meshes;
 }
diff --git a/manifold/src/csg_tree.cpp b/manifold/src/csg_tree.cpp
new file mode 100644
index 000000000..74a5a4adf
--- /dev/null
+++ b/manifold/src/csg_tree.cpp
@@ -0,0 +1,447 @@
+// Copyright 2022 Emmett Lalish, Chun Kit LAM
+//
+// 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 "csg_tree.h"
+
+#include <algorithm>
+
+#include "boolean3.h"
+#include "impl.h"
+#include "par.h"
+
+namespace {
+using namespace manifold;
+struct Transform4x3 {
+  const glm::mat4x3 transform;
+
+  __host__ __device__ glm::vec3 operator()(glm::vec3 position) {
+    return transform * glm::vec4(position, 1.0f);
+  }
+};
+
+struct TransformNormals {
+  const glm::mat3 transform;
+
+  __host__ __device__ glm::vec3 operator()(glm::vec3 normal) {
+    normal = glm::normalize(transform * normal);
+    if (isnan(normal.x)) normal = glm::vec3(0.0f);
+    return normal;
+  }
+};
+
+struct UpdateTriBary {
+  const int nextBary;
+
+  __host__ __device__ BaryRef operator()(BaryRef ref) {
+    for (int i : {0, 1, 2})
+      if (ref.vertBary[i] >= 0) ref.vertBary[i] += nextBary;
+    return ref;
+  }
+};
+
+struct UpdateHalfedge {
+  const int nextVert;
+  const int nextEdge;
+  const int nextFace;
+
+  __host__ __device__ Halfedge operator()(Halfedge edge) {
+    edge.startVert += nextVert;
+    edge.endVert += nextVert;
+    edge.pairedHalfedge += nextEdge;
+    edge.face += nextFace;
+    return edge;
+  }
+};
+
+struct CheckOverlap {
+  const Box *boxes;
+  const size_t i;
+  __host__ __device__ bool operator()(int j) {
+    return boxes[i].DoesOverlap(boxes[j]);
+  }
+};
+}  // namespace
+namespace manifold {
+
+std::shared_ptr<CsgNode> CsgNode::Translate(const glm::vec3 &t) const {
+  glm::mat4x3 transform(1.0f);
+  transform[3] += t;
+  return Transform(transform);
+}
+
+std::shared_ptr<CsgNode> CsgNode::Scale(const glm::vec3 &v) const {
+  glm::mat4x3 transform(1.0f);
+  for (int i : {0, 1, 2}) transform[i] *= v;
+  return Transform(transform);
+}
+
+std::shared_ptr<CsgNode> CsgNode::Rotate(float xDegrees, float yDegrees,
+                                         float zDegrees) const {
+  glm::mat3 rX(1.0f, 0.0f, 0.0f,                      //
+               0.0f, cosd(xDegrees), sind(xDegrees),  //
+               0.0f, -sind(xDegrees), cosd(xDegrees));
+  glm::mat3 rY(cosd(yDegrees), 0.0f, -sind(yDegrees),  //
+               0.0f, 1.0f, 0.0f,                       //
+               sind(yDegrees), 0.0f, cosd(yDegrees));
+  glm::mat3 rZ(cosd(zDegrees), sind(zDegrees), 0.0f,   //
+               -sind(zDegrees), cosd(zDegrees), 0.0f,  //
+               0.0f, 0.0f, 1.0f);
+  glm::mat4x3 transform(rZ * rY * rX);
+  return Transform(transform);
+}
+
+CsgLeafNode::CsgLeafNode() : pImpl_(std::make_shared<Manifold::Impl>()) {}
+
+CsgLeafNode::CsgLeafNode(std::shared_ptr<const Manifold::Impl> pImpl_)
+    : pImpl_(pImpl_) {}
+
+CsgLeafNode::CsgLeafNode(std::shared_ptr<const Manifold::Impl> pImpl_,
+                         glm::mat4x3 transform_)
+    : pImpl_(pImpl_), transform_(transform_) {}
+
+std::shared_ptr<const Manifold::Impl> CsgLeafNode::GetImpl() const {
+  if (transform_ == glm::mat4x3(1.0f)) return pImpl_;
+  pImpl_ =
+      std::make_shared<const Manifold::Impl>(pImpl_->Transform(transform_));
+  transform_ = glm::mat4x3(1.0f);
+  return pImpl_;
+}
+
+glm::mat4x3 CsgLeafNode::GetTransform() const { return transform_; }
+
+Box CsgLeafNode::GetBoundingBox() const {
+  return pImpl_->bBox_.Transform(transform_);
+}
+
+std::shared_ptr<CsgLeafNode> CsgLeafNode::ToLeafNode() const {
+  return std::make_shared<CsgLeafNode>(*this);
+}
+
+std::shared_ptr<CsgNode> CsgLeafNode::Transform(const glm::mat4x3 &m) const {
+  return std::make_shared<CsgLeafNode>(pImpl_, m * glm::mat4(transform_));
+}
+
+CsgNodeType CsgLeafNode::GetNodeType() const { return CsgNodeType::LEAF; }
+
+/**
+ * Efficient union of a set of pairwise disjoint meshes.
+ */
+Manifold::Impl CsgLeafNode::Compose(
+    const std::vector<std::shared_ptr<CsgLeafNode>> &nodes) {
+  float precision = -1;
+  int numVert = 0;
+  int numEdge = 0;
+  int numTri = 0;
+  int numBary = 0;
+  for (auto &node : nodes) {
+    float nodeOldScale = node->pImpl_->bBox_.Scale();
+    float nodeNewScale = node->GetBoundingBox().Scale();
+    float nodePrecision = node->pImpl_->precision_;
+    nodePrecision *= glm::max(1.0f, nodeNewScale / nodeOldScale);
+    nodePrecision = glm::max(nodePrecision, kTolerance * nodeNewScale);
+    if (!glm::isfinite(nodePrecision)) nodePrecision = -1;
+    precision = glm::max(precision, nodePrecision);
+
+    numVert += node->pImpl_->NumVert();
+    numEdge += node->pImpl_->NumEdge();
+    numTri += node->pImpl_->NumTri();
+    numBary += node->pImpl_->meshRelation_.barycentric.size();
+  }
+
+  Manifold::Impl combined;
+  combined.precision_ = precision;
+  combined.vertPos_.resize(numVert);
+  combined.halfedge_.resize(2 * numEdge);
+  combined.faceNormal_.resize(numTri);
+  combined.halfedgeTangent_.resize(2 * numEdge);
+  combined.meshRelation_.barycentric.resize(numBary);
+  combined.meshRelation_.triBary.resize(numTri);
+  auto policy = autoPolicy(numTri);
+
+  int nextVert = 0;
+  int nextEdge = 0;
+  int nextTri = 0;
+  int nextBary = 0;
+  for (auto &node : nodes) {
+    if (node->transform_ == glm::mat4x3(1.0f)) {
+      copy(policy, node->pImpl_->vertPos_.begin(), node->pImpl_->vertPos_.end(),
+           combined.vertPos_.begin() + nextVert);
+      copy(policy, node->pImpl_->faceNormal_.begin(),
+           node->pImpl_->faceNormal_.end(),
+           combined.faceNormal_.begin() + nextTri);
+    } else {
+      // no need to apply the transform to the node, just copy the vertices and
+      // face normals and apply transform on the fly
+      auto vertPosBegin = thrust::make_transform_iterator(
+          node->pImpl_->vertPos_.begin(), Transform4x3({node->transform_}));
+      glm::mat3 normalTransform =
+          glm::inverse(glm::transpose(glm::mat3(node->transform_)));
+      auto faceNormalBegin =
+          thrust::make_transform_iterator(node->pImpl_->faceNormal_.begin(),
+                                          TransformNormals({normalTransform}));
+      copy_n(policy, vertPosBegin, node->pImpl_->vertPos_.size(),
+             combined.vertPos_.begin() + nextVert);
+      copy_n(policy, faceNormalBegin, node->pImpl_->faceNormal_.size(),
+             combined.faceNormal_.begin() + nextTri);
+    }
+    copy(policy, node->pImpl_->halfedgeTangent_.begin(),
+         node->pImpl_->halfedgeTangent_.end(),
+         combined.halfedgeTangent_.begin() + nextEdge);
+    copy(policy, node->pImpl_->meshRelation_.barycentric.begin(),
+         node->pImpl_->meshRelation_.barycentric.end(),
+         combined.meshRelation_.barycentric.begin() + nextBary);
+    transform(policy, node->pImpl_->meshRelation_.triBary.begin(),
+              node->pImpl_->meshRelation_.triBary.end(),
+              combined.meshRelation_.triBary.begin() + nextTri,
+              UpdateTriBary({nextBary}));
+    transform(policy, node->pImpl_->halfedge_.begin(),
+              node->pImpl_->halfedge_.end(),
+              combined.halfedge_.begin() + nextEdge,
+              UpdateHalfedge({nextVert, nextEdge, nextTri}));
+
+    // Assign new IDs to triangles added in this iteration, to differentiate
+    // triangles coming from different manifolds.
+    VecDH<int> meshIDs;
+    VecDH<int> original;
+    for (auto &entry : node->pImpl_->meshRelation_.originalID) {
+      meshIDs.push_back(entry.first);
+      original.push_back(entry.second);
+    }
+    int meshIDStart = combined.meshRelation_.originalID.size();
+    combined.UpdateMeshIDs(meshIDs, original, nextTri,
+                           node->pImpl_->meshRelation_.triBary.size(),
+                           meshIDStart);
+
+    nextVert += node->pImpl_->NumVert();
+    nextEdge += 2 * node->pImpl_->NumEdge();
+    nextTri += node->pImpl_->NumTri();
+    nextBary += node->pImpl_->meshRelation_.barycentric.size();
+  }
+  // required to remove parts that are smaller than the precision
+  combined.SimplifyTopology();
+  combined.Finish();
+  return combined;
+}
+
+CsgOpNode::CsgOpNode() {}
+
+CsgOpNode::CsgOpNode(const std::vector<std::shared_ptr<CsgNode>> &children,
+                     Manifold::OpType op)
+    : children_(children) {
+  SetOp(op);
+  // opportunisticly flatten the tree without costly evaluation
+  GetChildren(false);
+}
+
+CsgOpNode::CsgOpNode(std::vector<std::shared_ptr<CsgNode>> &&children,
+                     Manifold::OpType op)
+    : children_(children) {
+  SetOp(op);
+  // opportunisticly flatten the tree without costly evaluation
+  GetChildren(false);
+}
+
+std::shared_ptr<CsgNode> CsgOpNode::Transform(const glm::mat4x3 &m) const {
+  auto node = std::make_shared<CsgOpNode>();
+  node->children_ = children_;
+  node->op_ = op_;
+  node->transform_ = m * glm::mat4(transform_);
+  node->simplified_ = simplified_;
+  node->flattened_ = flattened_;
+  return node;
+}
+
+std::shared_ptr<CsgLeafNode> CsgOpNode::ToLeafNode() const {
+  if (cache_ != nullptr) return cache_;
+  if (children_.empty()) return nullptr;
+  // turn the children into leaf nodes
+  GetChildren();
+  Manifold::OpType op;
+  switch (op_) {
+    case CsgNodeType::UNION:
+      op = Manifold::OpType::ADD;
+      break;
+    case CsgNodeType::DIFFERENCE:
+      op = Manifold::OpType::SUBTRACT;
+      break;
+    case CsgNodeType::INTERSECTION:
+      op = Manifold::OpType::INTERSECT;
+      break;
+    default:
+      throw std::runtime_error("unreachable CSG operation");
+      break;
+  }
+  auto a = std::static_pointer_cast<CsgLeafNode>(children_[0])->GetImpl();
+  for (int i = 1; i < children_.size(); i++) {
+    auto b = std::static_pointer_cast<CsgLeafNode>(children_[i])->GetImpl();
+    Boolean3 boolean(*a, *b, op);
+    a = std::make_shared<Manifold::Impl>(boolean.Result(op));
+  }
+  children_.clear();
+  children_.push_back(std::make_shared<CsgLeafNode>(a));
+  // children_ must contain only one CsgLeafNode now, and its Transform will
+  // give CsgLeafNode as well
+  cache_ = std::dynamic_pointer_cast<CsgLeafNode>(
+      children_.front()->Transform(transform_));
+  return cache_;
+}
+
+/**
+ * Efficient boolean operation on a set of nodes utilizing commutativity of the
+ * operation. Only supports union and intersection.
+ */
+void CsgOpNode::BatchBoolean(
+    Manifold::OpType operation,
+    std::vector<std::shared_ptr<const Manifold::Impl>> &results) {
+  assert(operation != Manifold::OpType::SUBTRACT);
+  auto cmpFn = [](std::shared_ptr<const Manifold::Impl> a,
+                  std::shared_ptr<const Manifold::Impl> b) {
+    // invert the order because we want a min heap
+    return a->NumVert() > b->NumVert();
+  };
+
+  // apply boolean operations starting from smaller meshes
+  // the assumption is that boolean operations on smaller meshes is faster,
+  // due to less data being copied and processed
+  std::make_heap(results.begin(), results.end(), cmpFn);
+  while (results.size() > 1) {
+    std::pop_heap(results.begin(), results.end(), cmpFn);
+    auto a = std::move(results.back());
+    results.pop_back();
+    std::pop_heap(results.begin(), results.end(), cmpFn);
+    auto b = std::move(results.back());
+    results.pop_back();
+    // boolean operation
+    Boolean3 boolean(*a, *b, operation);
+    results.push_back(
+        std::make_shared<const Manifold::Impl>(boolean.Result(operation)));
+    std::push_heap(results.begin(), results.end(), cmpFn);
+  }
+}
+
+/**
+ * Efficient union operation on a set of nodes by doing Compose as much as
+ * possible.
+ */
+void CsgOpNode::BatchUnion() const {
+  // INVARIANT: children_ is a vector of leaf nodes
+  // this kMaxUnionSize is a heuristic to avoid the pairwise disjoint check
+  // with O(n^2) complexity to take too long.
+  // If the number of children exceeded this limit, we will operate on chunks
+  // with size kMaxUnionSize.
+  constexpr int kMaxUnionSize = 1000;
+  while (children_.size() > 1) {
+    int start;
+    if (children_.size() > kMaxUnionSize) {
+      start = children_.size() - kMaxUnionSize;
+    } else {
+      start = 0;
+    }
+    VecDH<Box> boxes;
+    boxes.reserve(children_.size() - start);
+    for (int i = start; i < children_.size(); i++) {
+      boxes.push_back(std::dynamic_pointer_cast<CsgLeafNode>(children_[i])
+                          ->GetBoundingBox());
+    }
+    const Box *boxesD = boxes.cptrD();
+    // partition the children into a set of disjoint sets
+    // each set contains a set of children that are pairwise disjoint
+    std::vector<VecDH<size_t>> disjointSets;
+    for (size_t i = 0; i < boxes.size(); i++) {
+      auto lambda = [boxesD, i](const VecDH<size_t> &set) {
+        return find_if<decltype(set.end())>(
+                   autoPolicy(set.size()), set.begin(), set.end(),
+                   CheckOverlap({boxesD, i})) == set.end();
+      };
+      decltype(disjointSets.end()) it =
+          std::find_if(disjointSets.begin(), disjointSets.end(), lambda);
+      if (it == disjointSets.end()) {
+        disjointSets.push_back(std::vector<size_t>{i});
+      } else {
+        it->push_back(i);
+      }
+    }
+    // compose each set of disjoint children
+    std::vector<std::shared_ptr<const Manifold::Impl>> impls;
+    for (const auto &set : disjointSets) {
+      if (set.size() == 1) {
+        impls.push_back(
+            std::dynamic_pointer_cast<CsgLeafNode>(children_[start + set[0]])
+                ->GetImpl());
+      } else {
+        std::vector<std::shared_ptr<CsgLeafNode>> tmp;
+        for (size_t j : set) {
+          tmp.push_back(
+              std::dynamic_pointer_cast<CsgLeafNode>(children_[start + j]));
+        }
+        impls.push_back(
+            std::make_shared<const Manifold::Impl>(CsgLeafNode::Compose(tmp)));
+      }
+    }
+    BatchBoolean(Manifold::OpType::ADD, impls);
+    children_.erase(children_.begin() + start, children_.end());
+    children_.push_back(std::make_shared<CsgLeafNode>(impls.back()));
+    // move it to the front as we process from the back, and the newly added
+    // child should be quite complicated
+    std::swap(children_.front(), children_.back());
+  }
+}
+
+/**
+ * Flatten the children to a list of leaf nodes and return them.
+ * If finalize is true, the list will be guaranteed to be a list of leaf nodes
+ * (i.e. no ops). Otherwise, the list may contain ops.
+ * Note that this function will not apply the transform to children, as they may
+ * be shared with other nodes.
+ */
+std::vector<std::shared_ptr<CsgNode>> &CsgOpNode::GetChildren(
+    bool finalize) const {
+  if (children_.empty() || (simplified_ && !finalize) || flattened_)
+    return children_;
+  simplified_ = true;
+  flattened_ = finalize;
+  std::vector<std::shared_ptr<CsgNode>> newChildren;
+
+  CsgNodeType op = op_;
+  for (auto &child : children_) {
+    if (!finalize || child->GetNodeType() == CsgNodeType::LEAF) {
+      newChildren.push_back(child);
+    } else {
+      newChildren.push_back(child->ToLeafNode());
+    }
+    // special handling for difference: we treat it as first - (second + third +
+    // ...) so op = UNION after the first node
+    if (op == CsgNodeType::DIFFERENCE) op = CsgNodeType::UNION;
+  }
+  children_ = newChildren;
+  return children_;
+}
+
+void CsgOpNode::SetOp(Manifold::OpType op) {
+  switch (op) {
+    case Manifold::OpType::ADD:
+      op_ = CsgNodeType::UNION;
+      break;
+    case Manifold::OpType::SUBTRACT:
+      op_ = CsgNodeType::DIFFERENCE;
+      break;
+    case Manifold::OpType::INTERSECT:
+      op_ = CsgNodeType::INTERSECTION;
+      break;
+  }
+}
+
+glm::mat4x3 CsgOpNode::GetTransform() const { return transform_; }
+
+}  // namespace manifold
diff --git a/manifold/src/csg_tree.h b/manifold/src/csg_tree.h
new file mode 100644
index 000000000..04dadee6c
--- /dev/null
+++ b/manifold/src/csg_tree.h
@@ -0,0 +1,89 @@
+#pragma once
+#include "manifold.h"
+
+namespace manifold {
+
+enum class CsgNodeType { UNION, INTERSECTION, DIFFERENCE, LEAF };
+
+class CsgLeafNode;
+
+class CsgNode {
+ public:
+  virtual std::shared_ptr<CsgLeafNode> ToLeafNode() const = 0;
+  virtual std::shared_ptr<CsgNode> Transform(const glm::mat4x3 &m) const = 0;
+  virtual CsgNodeType GetNodeType() const = 0;
+  virtual glm::mat4x3 GetTransform() const = 0;
+
+  std::shared_ptr<CsgNode> Translate(const glm::vec3 &t) const;
+  std::shared_ptr<CsgNode> Scale(const glm::vec3 &s) const;
+  std::shared_ptr<CsgNode> Rotate(float xDegrees = 0, float yDegrees = 0,
+                                  float zDegrees = 0) const;
+};
+
+class CsgLeafNode final : public CsgNode {
+ public:
+  CsgLeafNode();
+  CsgLeafNode(std::shared_ptr<const Manifold::Impl> pImpl_);
+  CsgLeafNode(std::shared_ptr<const Manifold::Impl> pImpl_,
+              glm::mat4x3 transform_);
+
+  std::shared_ptr<const Manifold::Impl> GetImpl() const;
+
+  Box GetBoundingBox() const;
+
+  std::shared_ptr<CsgLeafNode> ToLeafNode() const override;
+
+  std::shared_ptr<CsgNode> Transform(const glm::mat4x3 &m) const override;
+
+  CsgNodeType GetNodeType() const override;
+
+  glm::mat4x3 GetTransform() const override;
+
+  static Manifold::Impl Compose(
+      const std::vector<std::shared_ptr<CsgLeafNode>> &nodes);
+
+ private:
+  mutable std::shared_ptr<const Manifold::Impl> pImpl_;
+  mutable glm::mat4x3 transform_ = glm::mat4x3(1.0f);
+};
+
+class CsgOpNode final : public CsgNode {
+ public:
+  CsgOpNode();
+
+  CsgOpNode(const std::vector<std::shared_ptr<CsgNode>> &children,
+            Manifold::OpType op);
+
+  CsgOpNode(std::vector<std::shared_ptr<CsgNode>> &&children,
+            Manifold::OpType op);
+
+  std::shared_ptr<CsgNode> Transform(const glm::mat4x3 &m) const override;
+
+  std::shared_ptr<CsgLeafNode> ToLeafNode() const override;
+
+  CsgNodeType GetNodeType() const override { return op_; }
+
+  glm::mat4x3 GetTransform() const override;
+
+ private:
+  CsgNodeType op_;
+  glm::mat4x3 transform_ = glm::mat4x3(1.0f);
+  // the following fields are for lazy evaluation, so they are mutable
+  mutable std::vector<std::shared_ptr<CsgNode>> children_;
+  mutable std::shared_ptr<CsgLeafNode> cache_ = nullptr;
+  mutable bool simplified_ = false;
+  mutable bool flattened_ = false;
+
+  void SetOp(Manifold::OpType);
+
+  static void BatchBoolean(
+      Manifold::OpType operation,
+      std::vector<std::shared_ptr<const Manifold::Impl>> &results);
+
+  void BatchUnion() const;
+
+  std::vector<std::shared_ptr<CsgNode>> &GetChildren(
+      bool finalize = true) const;
+};
+
+}  // namespace manifold
diff --git a/manifold/src/impl.cpp b/manifold/src/impl.cpp
index a522799b0..8fde13d1d 100644
--- a/manifold/src/impl.cpp
+++ b/manifold/src/impl.cpp
@@ -47,17 +47,18 @@ struct Normalize {
 struct Transform4x3 {
   const glm::mat4x3 transform;
 
-  __host__ __device__ void operator()(glm::vec3& position) {
-    position = transform * glm::vec4(position, 1.0f);
+  __host__ __device__ glm::vec3 operator()(glm::vec3 position) {
+    return transform * glm::vec4(position, 1.0f);
   }
 };
 
 struct TransformNormals {
   const glm::mat3 transform;
 
-  __host__ __device__ void operator()(glm::vec3& normal) {
+  __host__ __device__ glm::vec3 operator()(glm::vec3 normal) {
     normal = glm::normalize(transform * normal);
     if (isnan(normal.x)) normal = glm::vec3(0.0f);
+    return normal;
   }
 };
 
@@ -244,6 +245,7 @@ struct EdgeBox {
     thrust::get<0>(inout) = Box(vertPos[edge.first], vertPos[edge.second]);
   }
 };
+
 }  // namespace
 
 namespace manifold {
@@ -473,45 +475,42 @@ void Manifold::Impl::Update() {
   collider_.UpdateBoxes(faceBox);
 }
 
-void Manifold::Impl::ApplyTransform() const {
-  // This const_cast is here because these operations cancel out, leaving the
-  // state conceptually unchanged. This enables lazy transformation evaluation.
-  const_cast<Impl*>(this)->ApplyTransform();
-}
-
-/**
- * Bake the manifold's transform into its vertices. This function allows lazy
- * evaluation, which is important because often several transforms are applied
- * between operations.
- */
-void Manifold::Impl::ApplyTransform() {
-  if (transform_ == glm::mat4x3(1.0f)) return;
-  auto policy = autoPolicy(vertPos_.size());
-  for_each(policy, vertPos_.begin(), vertPos_.end(),
-           Transform4x3({transform_}));
+Manifold::Impl Manifold::Impl::Transform(const glm::mat4x3& transform_) const {
+  if (transform_ == glm::mat4x3(1.0f)) return *this;
+  auto policy = autoPolicy(NumVert());
+  Impl result;
+  result.collider_ = collider_;
+  result.meshRelation_ = meshRelation_;
+  result.precision_ = precision_;
+  result.bBox_ = bBox_;
+  result.halfedge_ = halfedge_;
+  result.halfedgeTangent_ = halfedgeTangent_;
+
+  result.vertPos_.resize(NumVert());
+  result.faceNormal_.resize(faceNormal_.size());
+  result.vertNormal_.resize(vertNormal_.size());
+  transform(policy, vertPos_.begin(), vertPos_.end(), result.vertPos_.begin(),
+            Transform4x3({transform_}));
 
   glm::mat3 normalTransform =
       glm::inverse(glm::transpose(glm::mat3(transform_)));
-  for_each(policy, faceNormal_.begin(), faceNormal_.end(),
-           TransformNormals({normalTransform}));
-  for_each(policy, vertNormal_.begin(), vertNormal_.end(),
-           TransformNormals({normalTransform}));
+  transform(policy, faceNormal_.begin(), faceNormal_.end(),
+            result.faceNormal_.begin(), TransformNormals({normalTransform}));
+  transform(policy, vertNormal_.begin(), vertNormal_.end(),
+            result.vertNormal_.begin(), TransformNormals({normalTransform}));
   // This optimization does a cheap collider update if the transform is
   // axis-aligned.
-  if (!collider_.Transform(transform_)) Update();
+  if (!result.collider_.Transform(transform_)) result.Update();
 
-  const float oldScale = bBox_.Scale();
-  transform_ = glm::mat4x3(1.0f);
-  CalculateBBox();
+  const float oldScale = result.bBox_.Scale();
+  result.CalculateBBox();
 
-  const float newScale = bBox_.Scale();
-  precision_ *= glm::max(1.0f, newScale / oldScale) *
-                glm::max(glm::length(transform_[0]),
-                         glm::max(glm::length(transform_[1]),
-                                  glm::length(transform_[2])));
+  const float newScale = result.bBox_.Scale();
+  result.precision_ *= glm::max(1.0f, newScale / oldScale);
 
   // Maximum of inherited precision loss and translational precision loss.
-  SetPrecision(precision_);
+  result.SetPrecision(result.precision_);
+  return result;
 }
 
 /**
diff --git a/manifold/src/impl.h b/manifold/src/impl.h
index a2e9ad395..5f21ebc98 100644
--- a/manifold/src/impl.h
+++ b/manifold/src/impl.h
@@ -48,7 +48,6 @@ struct Manifold::Impl {
   VecDH<glm::vec4> halfedgeTangent_;
   MeshRelationD meshRelation_;
   Collider collider_;
-  glm::mat4x3 transform_ = glm::mat4x3(1.0f);
 
   static std::atomic<int> meshIDCounter_;
 
@@ -73,8 +72,7 @@ struct Manifold::Impl {
                      int startTri = 0, int n = -1, int startID = 0);
 
   void Update();
-  void ApplyTransform() const;
-  void ApplyTransform();
+  Impl Transform(const glm::mat4x3& transform) const;
   SparseIndices EdgeCollisions(const Impl& B) const;
   SparseIndices VertexCollisionsZ(const VecDH<glm::vec3>& vertsIn) const;
 
diff --git a/manifold/src/manifold.cpp b/manifold/src/manifold.cpp
index d8678b425..99ee1b3cf 100644
--- a/manifold/src/manifold.cpp
+++ b/manifold/src/manifold.cpp
@@ -13,6 +13,7 @@
 // limitations under the License.
 
 #include "boolean3.h"
+#include "csg_tree.h"
 #include "impl.h"
 #include "par.h"
 
@@ -45,7 +46,7 @@ Manifold Halfspace(Box bBox, glm::vec3 normal, float originOffset) {
       Manifold::Cube(glm::vec3(2.0f), true).Translate({1.0f, 0.0f, 0.0f});
   float size = glm::length(bBox.Center() - normal * originOffset) +
                0.5f * glm::length(bBox.Size());
-  cutter.Scale(glm::vec3(size)).Translate({originOffset, 0.0f, 0.0f});
+  cutter = cutter.Scale(glm::vec3(size)).Translate({originOffset, 0.0f, 0.0f});
   float yDeg = glm::degrees(-glm::asin(normal.z));
   float zDeg = glm::degrees(glm::atan(normal.y, normal.x));
   return cutter.Rotate(0.0f, yDeg, zDeg);
@@ -58,20 +59,32 @@ namespace manifold {
  * Construct an empty Manifold.
  *
  */
-Manifold::Manifold() : pImpl_{std::make_unique<Impl>()} {}
+Manifold::Manifold() : pNode_{std::make_shared<CsgLeafNode>()} {}
 Manifold::~Manifold() = default;
 Manifold::Manifold(Manifold&&) noexcept = default;
 Manifold& Manifold::operator=(Manifold&&) noexcept = default;
 
-Manifold::Manifold(const Manifold& other) : pImpl_(new Impl(*other.pImpl_)) {}
+Manifold::Manifold(const Manifold& other) : pNode_(other.pNode_) {}
+
+Manifold::Manifold(std::shared_ptr<CsgNode> pNode) : pNode_(pNode) {}
+
+Manifold::Manifold(std::shared_ptr<Impl> pImpl_)
+    : pNode_(std::make_shared<CsgLeafNode>(pImpl_)) {}
 
 Manifold& Manifold::operator=(const Manifold& other) {
   if (this != &other) {
-    pImpl_.reset(new Impl(*other.pImpl_));
+    pNode_ = other.pNode_;
   }
   return *this;
 }
 
+CsgLeafNode& Manifold::GetCsgLeafNode() const {
+  if (pNode_->GetNodeType() != CsgNodeType::LEAF) {
+    pNode_ = pNode_->ToLeafNode();
+  }
+  return *std::static_pointer_cast<CsgLeafNode>(pNode_);
+}
+
 /**
  * Convert a Mesh into a Manifold. Will throw a topologyErr exception if the
  * input is not an oriented 2-manifold. Will collapse degenerate triangles and
@@ -102,29 +115,29 @@ Manifold::Manifold(const Mesh& mesh,
                    const std::vector<glm::ivec3>& triProperties,
                    const std::vector<float>& properties,
                    const std::vector<float>& propertyTolerance)
-    : pImpl_{std::make_unique<Impl>(mesh, triProperties, properties,
-                                    propertyTolerance)} {}
+    : pNode_(std::make_shared<CsgLeafNode>(std::make_shared<Impl>(
+          mesh, triProperties, properties, propertyTolerance))) {}
 
 /**
  * This returns a Mesh of simple vectors of vertices and triangles suitable for
  * saving or other operations outside of the context of this library.
  */
 Mesh Manifold::GetMesh() const {
-  pImpl_->ApplyTransform();
+  const Impl& impl = *GetCsgLeafNode().GetImpl();
 
   Mesh result;
-  result.vertPos.insert(result.vertPos.end(), pImpl_->vertPos_.begin(),
-                        pImpl_->vertPos_.end());
-  result.vertNormal.insert(result.vertNormal.end(), pImpl_->vertNormal_.begin(),
-                           pImpl_->vertNormal_.end());
+  result.vertPos.insert(result.vertPos.end(), impl.vertPos_.begin(),
+                        impl.vertPos_.end());
+  result.vertNormal.insert(result.vertNormal.end(), impl.vertNormal_.begin(),
+                           impl.vertNormal_.end());
   result.halfedgeTangent.insert(result.halfedgeTangent.end(),
-                                pImpl_->halfedgeTangent_.begin(),
-                                pImpl_->halfedgeTangent_.end());
+                                impl.halfedgeTangent_.begin(),
+                                impl.halfedgeTangent_.end());
 
   result.triVerts.resize(NumTri());
   // note that `triVerts` is `std::vector`, so we cannot use thrust::device
   thrust::for_each_n(thrust::host, zip(result.triVerts.begin(), countAt(0)),
-                     NumTri(), MakeTri({pImpl_->halfedge_.cptrH()}));
+                     NumTri(), MakeTri({impl.halfedge_.cptrH()}));
 
   return result;
 }
@@ -193,26 +206,24 @@ int Manifold::GetCircularSegments(float radius) {
 /**
  * Does the Manifold have any triangles?
  */
-bool Manifold::IsEmpty() const { return pImpl_->IsEmpty(); }
+bool Manifold::IsEmpty() const { return GetCsgLeafNode().GetImpl()->IsEmpty(); }
 /**
  * The number of vertices in the Manifold.
  */
-int Manifold::NumVert() const { return pImpl_->NumVert(); }
+int Manifold::NumVert() const { return GetCsgLeafNode().GetImpl()->NumVert(); }
 /**
  * The number of edges in the Manifold.
  */
-int Manifold::NumEdge() const { return pImpl_->NumEdge(); }
+int Manifold::NumEdge() const { return GetCsgLeafNode().GetImpl()->NumEdge(); }
 /**
  * The number of triangles in the Manifold.
  */
-int Manifold::NumTri() const { return pImpl_->NumTri(); }
+int Manifold::NumTri() const { return GetCsgLeafNode().GetImpl()->NumTri(); }
 
 /**
  * Returns the axis-aligned bounding box of all the Manifold's vertices.
  */
-Box Manifold::BoundingBox() const {
-  return pImpl_->bBox_.Transform(pImpl_->transform_);
-}
+Box Manifold::BoundingBox() const { return GetCsgLeafNode().GetBoundingBox(); }
 
 /**
  * Returns the precision of this Manifold's vertices, which tracks the
@@ -222,8 +233,7 @@ Box Manifold::BoundingBox() const {
  * [&epsilon;-valid](https://github.com/elalish/manifold/wiki/Manifold-Library#definition-of-%CE%B5-valid).
  */
 float Manifold::Precision() const {
-  pImpl_->ApplyTransform();
-  return pImpl_->precision_;
+  return GetCsgLeafNode().GetImpl()->precision_;
 }
 
 /**
@@ -242,7 +252,9 @@ int Manifold::Genus() const {
  * means degenerate manifolds can by identified by testing these properties as
  * == 0.
  */
-Properties Manifold::GetProperties() const { return pImpl_->GetProperties(); }
+Properties Manifold::GetProperties() const {
+  return GetCsgLeafNode().GetImpl()->GetProperties();
+}
 
 /**
  * Curvature is the inverse of the radius of curvature, and signed such that
@@ -252,7 +264,9 @@ Properties Manifold::GetProperties() const { return pImpl_->GetProperties(); }
  * curvature is their sum. This approximates them for every vertex (returned as
  * vectors in the structure) and also returns their minimum and maximum values.
  */
-Curvature Manifold::GetCurvature() const { return pImpl_->GetCurvature(); }
+Curvature Manifold::GetCurvature() const {
+  return GetCsgLeafNode().GetImpl()->GetCurvature();
+}
 
 /**
  * Gets the relationship to the previous mesh, for the purpose of assinging
@@ -265,7 +279,7 @@ Curvature Manifold::GetCurvature() const { return pImpl_->GetCurvature(); }
  */
 MeshRelation Manifold::GetMeshRelation() const {
   MeshRelation out;
-  const auto& relation = pImpl_->meshRelation_;
+  const auto& relation = GetCsgLeafNode().GetImpl()->meshRelation_;
   out.triBary.insert(out.triBary.end(), relation.triBary.begin(),
                      relation.triBary.end());
   for (auto& bary : out.triBary) {
@@ -284,8 +298,8 @@ MeshRelation Manifold::GetMeshRelation() const {
  */
 std::vector<int> Manifold::GetMeshIDs() const {
   std::vector<int> out;
-  out.reserve(pImpl_->meshRelation_.originalID.size());
-  for (auto& entry : pImpl_->meshRelation_.originalID) {
+  out.reserve(GetCsgLeafNode().GetImpl()->meshRelation_.originalID.size());
+  for (auto& entry : GetCsgLeafNode().GetImpl()->meshRelation_.originalID) {
     out.push_back(entry.second);
   }
   return out;
@@ -304,31 +318,38 @@ std::vector<int> Manifold::GetMeshIDs() const {
  * should instead call GetMesh(), calculate your properties and use these to
  * construct a new manifold.
  *
- * @returns New MeshID
+ * @returns New Mesh
  */
-int Manifold::SetAsOriginal() {
-  int meshID = pImpl_->InitializeNewReference();
-  return meshID;
+Manifold Manifold::AsOriginal() const {
+  auto newImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+  newImpl->InitializeNewReference();
+  return Manifold(std::make_shared<CsgLeafNode>(newImpl));
 }
 
 /**
  * Should always be true. Also checks saneness of the internal data structures.
  */
-bool Manifold::IsManifold() const { return pImpl_->IsManifold(); }
+bool Manifold::IsManifold() const {
+  return GetCsgLeafNode().GetImpl()->IsManifold();
+}
 
 /**
  * The triangle normal vectors are saved over the course of operations rather
  * than recalculated to avoid rounding error. This checks that triangles still
  * match their normal vectors within Precision().
  */
-bool Manifold::MatchesTriNormals() const { return pImpl_->MatchesTriNormals(); }
+bool Manifold::MatchesTriNormals() const {
+  return GetCsgLeafNode().GetImpl()->MatchesTriNormals();
+}
 
 /**
  * The number of triangles that are colinear within Precision(). This library
  * attempts to remove all of these, but it cannot always remove all of them
  * without changing the mesh by too much.
  */
-int Manifold::NumDegenerateTris() const { return pImpl_->NumDegenerateTris(); }
+int Manifold::NumDegenerateTris() const {
+  return GetCsgLeafNode().GetImpl()->NumDegenerateTris();
+}
 
 /**
  * This is a checksum-style verification of the collider, simply returning the
@@ -337,13 +358,12 @@ int Manifold::NumDegenerateTris() const { return pImpl_->NumDegenerateTris(); }
  * @param other A Manifold to overlap with.
  */
 int Manifold::NumOverlaps(const Manifold& other) const {
-  pImpl_->ApplyTransform();
-  other.pImpl_->ApplyTransform();
-
-  SparseIndices overlaps = pImpl_->EdgeCollisions(*other.pImpl_);
+  SparseIndices overlaps = GetCsgLeafNode().GetImpl()->EdgeCollisions(
+      *other.GetCsgLeafNode().GetImpl());
   int num_overlaps = overlaps.size();
 
-  overlaps = other.pImpl_->EdgeCollisions(*pImpl_);
+  overlaps = other.GetCsgLeafNode().GetImpl()->EdgeCollisions(
+      *GetCsgLeafNode().GetImpl());
   return num_overlaps += overlaps.size();
 }
 
@@ -353,9 +373,8 @@ int Manifold::NumOverlaps(const Manifold& other) const {
  *
  * @param v The vector to add to every vertex.
  */
-Manifold& Manifold::Translate(glm::vec3 v) {
-  pImpl_->transform_[3] += v;
-  return *this;
+Manifold Manifold::Translate(glm::vec3 v) const {
+  return Manifold(pNode_->Translate(v));
 }
 
 /**
@@ -364,11 +383,8 @@ Manifold& Manifold::Translate(glm::vec3 v) {
  *
  * @param v The vector to multiply every vertex by per component.
  */
-Manifold& Manifold::Scale(glm::vec3 v) {
-  glm::mat3 s(1.0f);
-  for (int i : {0, 1, 2}) s[i] *= v;
-  pImpl_->transform_ = s * pImpl_->transform_;
-  return *this;
+Manifold Manifold::Scale(glm::vec3 v) const {
+  return Manifold(pNode_->Scale(v));
 }
 
 /**
@@ -383,18 +399,9 @@ Manifold& Manifold::Scale(glm::vec3 v) {
  * @param yDegrees Second rotation, degrees about the Y-axis.
  * @param zDegrees Third rotation, degrees about the Z-axis.
  */
-Manifold& Manifold::Rotate(float xDegrees, float yDegrees, float zDegrees) {
-  glm::mat3 rX(1.0f, 0.0f, 0.0f,                      //
-               0.0f, cosd(xDegrees), sind(xDegrees),  //
-               0.0f, -sind(xDegrees), cosd(xDegrees));
-  glm::mat3 rY(cosd(yDegrees), 0.0f, -sind(yDegrees),  //
-               0.0f, 1.0f, 0.0f,                       //
-               sind(yDegrees), 0.0f, cosd(yDegrees));
-  glm::mat3 rZ(cosd(zDegrees), sind(zDegrees), 0.0f,   //
-               -sind(zDegrees), cosd(zDegrees), 0.0f,  //
-               0.0f, 0.0f, 1.0f);
-  pImpl_->transform_ = rZ * rY * rX * pImpl_->transform_;
-  return *this;
+Manifold Manifold::Rotate(float xDegrees, float yDegrees,
+                          float zDegrees) const {
+  return Manifold(pNode_->Rotate(xDegrees, yDegrees, zDegrees));
 }
 
 /**
@@ -404,10 +411,8 @@ Manifold& Manifold::Rotate(float xDegrees, float yDegrees, float zDegrees) {
  *
  * @param m The affine transform matrix to apply to all the vertices.
  */
-Manifold& Manifold::Transform(const glm::mat4x3& m) {
-  glm::mat4 old(pImpl_->transform_);
-  pImpl_->transform_ = m * old;
-  return *this;
+Manifold Manifold::Transform(const glm::mat4x3& m) const {
+  return Manifold(pNode_->Transform(m));
 }
 
 /**
@@ -419,15 +424,15 @@ Manifold& Manifold::Transform(const glm::mat4x3& m) {
  *
  * @param warpFunc A function that modifies a given vertex position.
  */
-Manifold& Manifold::Warp(std::function<void(glm::vec3&)> warpFunc) {
-  pImpl_->ApplyTransform();
-  thrust::for_each_n(thrust::host, pImpl_->vertPos_.begin(), NumVert(),
+Manifold Manifold::Warp(std::function<void(glm::vec3&)> warpFunc) const {
+  auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+  thrust::for_each_n(thrust::host, pImpl->vertPos_.begin(), NumVert(),
                      warpFunc);
-  pImpl_->Update();
-  pImpl_->faceNormal_.resize(0);  // force recalculation of triNormal
-  pImpl_->CalculateNormals();
-  pImpl_->SetPrecision();
-  return *this;
+  pImpl->Update();
+  pImpl->faceNormal_.resize(0);  // force recalculation of triNormal
+  pImpl->CalculateNormals();
+  pImpl->SetPrecision();
+  return Manifold(std::make_shared<CsgLeafNode>(pImpl));
 }
 
 /**
@@ -440,9 +445,10 @@ Manifold& Manifold::Warp(std::function<void(glm::vec3&)> warpFunc) {
  *
  * @param n The number of pieces to split every edge into. Must be > 1.
  */
-Manifold& Manifold::Refine(int n) {
-  pImpl_->Refine(n);
-  return *this;
+Manifold Manifold::Refine(int n) const {
+  auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+  pImpl->Refine(n);
+  return Manifold(std::make_shared<CsgLeafNode>(pImpl));
 }
 
 /**
@@ -460,12 +466,8 @@ Manifold& Manifold::Refine(int n) {
  * @param op The type of operation to perform.
  */
 Manifold Manifold::Boolean(const Manifold& second, OpType op) const {
-  pImpl_->ApplyTransform();
-  second.pImpl_->ApplyTransform();
-  Boolean3 boolean(*pImpl_, *second.pImpl_, op);
-  Manifold result;
-  result.pImpl_ = std::make_unique<Impl>(boolean.Result(op));
-  return result;
+  std::vector<std::shared_ptr<CsgNode>> children({pNode_, second.pNode_});
+  return Manifold(std::make_shared<CsgOpNode>(children, op));
 }
 
 /**
@@ -521,15 +523,15 @@ Manifold& Manifold::operator^=(const Manifold& Q) {
  * @param cutter
  */
 std::pair<Manifold, Manifold> Manifold::Split(const Manifold& cutter) const {
-  pImpl_->ApplyTransform();
-  cutter.pImpl_->ApplyTransform();
-  Boolean3 boolean(*pImpl_, *cutter.pImpl_, OpType::SUBTRACT);
-  std::pair<Manifold, Manifold> result;
-  result.first.pImpl_ =
-      std::make_unique<Impl>(boolean.Result(OpType::INTERSECT));
-  result.second.pImpl_ =
-      std::make_unique<Impl>(boolean.Result(OpType::SUBTRACT));
-  return result;
+  auto impl1 = GetCsgLeafNode().GetImpl();
+  auto impl2 = cutter.GetCsgLeafNode().GetImpl();
+
+  Boolean3 boolean(*impl1, *impl2, OpType::SUBTRACT);
+  auto result1 = std::make_shared<CsgLeafNode>(
+      std::make_unique<Impl>(boolean.Result(OpType::INTERSECT)));
+  auto result2 = std::make_shared<CsgLeafNode>(
+      std::make_unique<Impl>(boolean.Result(OpType::SUBTRACT)));
+  return std::make_pair(Manifold(result1), Manifold(result2));
 }
 
 /**
@@ -556,7 +558,7 @@ std::pair<Manifold, Manifold> Manifold::SplitByPlane(glm::vec3 normal,
  * direction of the normal vector.
  */
 Manifold Manifold::TrimByPlane(glm::vec3 normal, float originOffset) const {
-  pImpl_->ApplyTransform();
   return *this ^ Halfspace(BoundingBox(), normal, originOffset);
 }
+
 }  // namespace manifold
diff --git a/manifold/src/properties.cpp b/manifold/src/properties.cpp
index 9377ae07f..931b85997 100644
--- a/manifold/src/properties.cpp
+++ b/manifold/src/properties.cpp
@@ -251,7 +251,6 @@ int Manifold::Impl::NumDegenerateTris() const {
 
 Properties Manifold::Impl::GetProperties() const {
   if (IsEmpty()) return {0, 0};
-  ApplyTransform();
   auto areaVolume = transform_reduce<thrust::pair<float, float>>(
       autoPolicy(NumTri()), countAt(0), countAt(NumTri()),
       FaceAreaVolume({halfedge_.cptrD(), vertPos_.cptrD(), precision_}),
@@ -262,7 +261,6 @@ Properties Manifold::Impl::GetProperties() const {
 Curvature Manifold::Impl::GetCurvature() const {
   Curvature result;
   if (IsEmpty()) return result;
-  ApplyTransform();
   VecDH<float> vertMeanCurvature(NumVert(), 0);
   VecDH<float> vertGaussianCurvature(NumVert(), glm::two_pi<float>());
   VecDH<float> vertArea(NumVert(), 0);
diff --git a/samples/src/bracelet.cpp b/samples/src/bracelet.cpp
index c8b048550..2ad28a15f 100644
--- a/samples/src/bracelet.cpp
+++ b/samples/src/bracelet.cpp
@@ -31,12 +31,12 @@ Manifold Base(float width, float radius, float decorRadius, float twistRadius,
                    decorRadius * sind(dPhiDeg * i)),
          0});
   }
-  Manifold decor = std::move(Manifold::Extrude(circle, width, nDivision, 180)
-                                 .Scale({1.0f, 0.5f, 1.0f})
-                                 .Translate({0.0f, radius, 0.0f}));
+  Manifold decor = Manifold::Extrude(circle, width, nDivision, 180)
+                       .Scale({1.0f, 0.5f, 1.0f})
+                       .Translate({0.0f, radius, 0.0f});
 
   for (int i = 0; i < nDecor; ++i) {
-    base += decor.Rotate(0, 0, 360.0f / nDecor);
+    base += decor.Rotate(0, 0, (360.0f / nDecor) * i);
   }
 
   Polygons stretch(1);
@@ -52,7 +52,7 @@ Manifold Base(float width, float radius, float decorRadius, float twistRadius,
   }
 
   base = Manifold::Extrude(stretch, width) ^ base;
-  base.SetAsOriginal();
+  base = base.AsOriginal();
 
   return base;
 }
diff --git a/samples/src/knot.cpp b/samples/src/knot.cpp
index dcc8bc569..a0a009348 100644
--- a/samples/src/knot.cpp
+++ b/samples/src/knot.cpp
@@ -58,31 +58,32 @@ Manifold TorusKnot(int p, int q, float majorRadius, float minorRadius,
 
   Manifold knot = Manifold::Revolve(circle, m);
 
-  knot.Warp([p, q, majorRadius, minorRadius, threadRadius](glm::vec3& v) {
-    float psi = q * atan2(v.x, v.y);
-    float theta = psi * p / q;
-    glm::vec2 xy = glm::vec2(v);
-    float x1 = sqrt(glm::dot(xy, xy));
-    float phi = atan2(x1 - 2, v.z);
-    v = glm::vec3(cos(phi), 0.0f, sin(phi));
-    v *= threadRadius;
-    float r = majorRadius + minorRadius * cos(theta);
-    v = glm::rotateX(v, -float(atan2(p * minorRadius, q * r)));
-    v.x += minorRadius;
-    v = glm::rotateY(v, theta);
-    v.x += majorRadius;
-    v = glm::rotateZ(v, psi);
-  });
+  knot =
+      knot.Warp([p, q, majorRadius, minorRadius, threadRadius](glm::vec3& v) {
+        float psi = q * atan2(v.x, v.y);
+        float theta = psi * p / q;
+        glm::vec2 xy = glm::vec2(v);
+        float x1 = sqrt(glm::dot(xy, xy));
+        float phi = atan2(x1 - 2, v.z);
+        v = glm::vec3(cos(phi), 0.0f, sin(phi));
+        v *= threadRadius;
+        float r = majorRadius + minorRadius * cos(theta);
+        v = glm::rotateX(v, -float(atan2(p * minorRadius, q * r)));
+        v.x += minorRadius;
+        v = glm::rotateY(v, theta);
+        v.x += majorRadius;
+        v = glm::rotateZ(v, psi);
+      });
 
   if (kLoops > 1) {
     std::vector<Manifold> knots;
     for (float k = 0; k < kLoops; ++k) {
-      knots.emplace_back(knot);
-      knots.back().Rotate(0, 0, 360.0f * (k / kLoops) * (q / float(p)));
+      knots.push_back(
+          knot.Rotate(0, 0, 360.0f * (k / kLoops) * (q / float(p))));
     }
     knot = Manifold::Compose(knots);
   }
 
   return knot;
 }
-}  // namespace manifold
\ No newline at end of file
+}  // namespace manifold
diff --git a/samples/src/menger_sponge.cpp b/samples/src/menger_sponge.cpp
index 2cb4798ed..a9168cf27 100644
--- a/samples/src/menger_sponge.cpp
+++ b/samples/src/menger_sponge.cpp
@@ -21,8 +21,8 @@ using namespace manifold;
 void Fractal(std::vector<Manifold>& holes, Manifold& hole, float w,
              glm::vec2 position, int depth, int maxDepth) {
   w /= 3;
-  holes.emplace_back(hole);
-  holes.back().Scale({w, w, 1.0f}).Translate(glm::vec3(position, 0.0f));
+  holes.push_back(
+      hole.Scale({w, w, 1.0f}).Translate(glm::vec3(position, 0.0f)));
   if (depth == maxDepth) return;
 
   glm::vec2 offsets[8] = {{-w, -w}, {-w, 0.0f}, {-w, w}, {0.0f, w},
@@ -50,9 +50,11 @@ Manifold MengerSponge(int n) {
   Manifold hole = Manifold::Compose(holes);
 
   result -= hole;
-  result -= hole.Rotate(90);
-  result -= hole.Rotate(0, 0, 90);
-  result.SetAsOriginal();
+  hole = hole.Rotate(90);
+  result -= hole;
+  hole = hole.Rotate(0, 0, 90);
+  result -= hole;
+  result = result.AsOriginal();
   return result;
 }
 }  // namespace manifold
diff --git a/samples/src/rounded_frame.cpp b/samples/src/rounded_frame.cpp
index 3e8ee1cb5..229c3f0ce 100644
--- a/samples/src/rounded_frame.cpp
+++ b/samples/src/rounded_frame.cpp
@@ -31,20 +31,18 @@ Manifold RoundedFrame(float edgeLength, float radius, int circularSegments) {
   Manifold corner = Manifold::Sphere(radius, circularSegments);
 
   Manifold edge1 = corner + edge;
-  edge1.Rotate(-90).Translate({-edgeLength / 2, -edgeLength / 2, 0});
+  edge1 = edge1.Rotate(-90).Translate({-edgeLength / 2, -edgeLength / 2, 0});
 
-  Manifold edge2 = edge1;
-  edge2.Rotate(0, 0, 180);
+  Manifold edge2 = edge1.Rotate(0, 0, 180);
   edge2 += edge1;
   edge2 += edge.Translate({-edgeLength / 2, -edgeLength / 2, 0});
 
-  Manifold edge4 = edge2;
-  edge4.Rotate(0, 0, 90);
+  Manifold edge4 = edge2.Rotate(0, 0, 90);
   edge4 += edge2;
 
   Manifold frame = edge4.Translate({0, 0, -edgeLength / 2});
-  frame += edge4.Rotate(180);
+  frame += frame.Rotate(180);
 
   return frame;
 }
-}  // namespace manifold
\ No newline at end of file
+}  // namespace manifold
diff --git a/test/mesh_test.cpp b/test/mesh_test.cpp
index 218257ee4..a6097e682 100644
--- a/test/mesh_test.cpp
+++ b/test/mesh_test.cpp
@@ -173,8 +173,7 @@ TEST(Manifold, DISABLED_Determinism) {
   Manifold manifold(ImportMesh("data/gyroidpuzzle.ply"));
   EXPECT_TRUE(manifold.IsManifold());
 
-  Manifold manifold1 = manifold;
-  manifold1.Translate(glm::vec3(5.0f));
+  Manifold manifold1 = manifold.Translate(glm::vec3(5.0f));
   int num_overlaps = manifold.NumOverlaps(manifold1);
   ASSERT_EQ(num_overlaps, 229611);
 
@@ -273,7 +272,7 @@ TEST(Manifold, Revolve2) {
 TEST(Manifold, Smooth) {
   Manifold tet = Manifold::Tetrahedron();
   Manifold smooth = Manifold::Smooth(tet.GetMesh());
-  smooth.Refine(100);
+  smooth = smooth.Refine(100);
   ExpectMeshes(smooth, {{20002, 40000}});
   auto prop = smooth.GetProperties();
   EXPECT_NEAR(prop.volume, 17.38, 0.1);
@@ -312,7 +311,7 @@ TEST(Manifold, ManualSmooth) {
   smooth.halfedgeTangent[16] = {0, 0, 0, 1};
   smooth.halfedgeTangent[18] = {0, 0, 0, 1};
   Manifold interp(smooth);
-  interp.Refine(100);
+  interp = interp.Refine(100);
 
   ExpectMeshes(interp, {{40002, 80000}});
   auto prop = interp.GetProperties();
@@ -343,7 +342,7 @@ TEST(Manifold, ManualSmooth) {
 
 TEST(Manifold, Csaszar) {
   Manifold csaszar = Manifold::Smooth(Csaszar());
-  csaszar.Refine(100);
+  csaszar = csaszar.Refine(100);
   ExpectMeshes(csaszar, {{70000, 140000}});
   auto prop = csaszar.GetProperties();
   EXPECT_NEAR(prop.volume, 84699, 10);
@@ -381,7 +380,7 @@ TEST(Manifold, GetProperties) {
   EXPECT_FLOAT_EQ(prop.volume, 1.0f);
   EXPECT_FLOAT_EQ(prop.surfaceArea, 6.0f);
 
-  cube.Scale(glm::vec3(-1.0f));
+  cube = cube.Scale(glm::vec3(-1.0f));
   prop = cube.GetProperties();
   EXPECT_FLOAT_EQ(prop.volume, -1.0f);
   EXPECT_FLOAT_EQ(prop.surfaceArea, 6.0f);
@@ -390,9 +389,9 @@ TEST(Manifold, GetProperties) {
 TEST(Manifold, Precision) {
   Manifold cube = Manifold::Cube();
   EXPECT_FLOAT_EQ(cube.Precision(), kTolerance);
-  cube.Scale({0.1, 1, 10});
+  cube = cube.Scale({0.1, 1, 10});
   EXPECT_FLOAT_EQ(cube.Precision(), 10 * kTolerance);
-  cube.Translate({-100, -10, -1});
+  cube = cube.Translate({-100, -10, -1});
   EXPECT_FLOAT_EQ(cube.Precision(), 100 * kTolerance);
 }
 
@@ -415,7 +414,7 @@ TEST(Manifold, GetCurvature) {
     EXPECT_NEAR(curvature.minGaussianCurvature, 1, precision);
     EXPECT_NEAR(curvature.maxGaussianCurvature, 1, precision);
 
-    sphere.Scale(glm::vec3(2.0f));
+    sphere = sphere.Scale(glm::vec3(2.0f));
     curvature = sphere.GetCurvature();
     EXPECT_NEAR(curvature.minMeanCurvature, 1, precision);
     EXPECT_NEAR(curvature.maxMeanCurvature, 1, precision);
@@ -431,7 +430,7 @@ TEST(Manifold, GetCurvature) {
 TEST(Manifold, Transform) {
   Manifold cube = Manifold::Cube({1, 2, 3});
   Manifold cube2 = cube;
-  cube.Rotate(30, 40, 50).Scale({6, 5, 4}).Translate({1, 2, 3});
+  cube = cube.Rotate(30, 40, 50).Scale({6, 5, 4}).Translate({1, 2, 3});
 
   glm::mat3 rX(1.0f, 0.0f, 0.0f,          //
                0.0f, cosd(30), sind(30),  //
@@ -448,7 +447,7 @@ TEST(Manifold, Transform) {
   s[2][2] = 4;
   glm::mat4x3 transform = glm::mat4x3(s * rZ * rY * rX);
   transform[3] = glm::vec3(1, 2, 3);
-  cube2.Transform(transform);
+  cube2 = cube2.Transform(transform);
 
   Identical(cube.GetMesh(), cube2.GetMesh());
 }
@@ -491,8 +490,7 @@ TEST(Boolean, Tetra) {
   EXPECT_TRUE(tetra.IsManifold());
 
   Manifold tetra2 = tetra;
-  tetra2.Translate(glm::vec3(0.5f));
-  tetra2.SetAsOriginal();
+  tetra2 = tetra2.Translate(glm::vec3(0.5f)).AsOriginal();
   Manifold result = tetra2 - tetra;
 
   ExpectMeshes(result, {{8, 12}});
@@ -533,11 +531,11 @@ TEST(Boolean, Perturb) {
 
 TEST(Boolean, Coplanar) {
   Manifold cylinder = Manifold::Cylinder(1.0f, 1.0f);
-  Manifold cylinder2 = cylinder;
-  cylinder2.SetAsOriginal();
-  Manifold out = cylinder - cylinder2.Scale({0.5f, 0.5f, 1.0f})
-                                .Rotate(0, 0, 15)
-                                .Translate({0.25f, 0.25f, 0.0f});
+  Manifold cylinder2 = cylinder.AsOriginal();
+  cylinder2 = cylinder2.Scale({0.5f, 0.5f, 1.0f})
+                  .Rotate(0, 0, 15)
+                  .Translate({0.25f, 0.25f, 0.0f});
+  Manifold out = cylinder - cylinder2;
   ExpectMeshes(out, {{32, 64}});
   EXPECT_EQ(out.NumDegenerateTris(), 0);
   EXPECT_EQ(out.Genus(), 1);
@@ -549,9 +547,8 @@ TEST(Boolean, Coplanar) {
 
 TEST(Boolean, MultiCoplanar) {
   Manifold cube = Manifold::Cube();
-  Manifold cube2 = cube;
-  Manifold first = cube - cube2.Translate({0.3f, 0.3f, 0.0f});
-  cube.Translate({-0.3f, -0.3f, 0.0f});
+  Manifold first = cube - cube.Translate({0.3f, 0.3f, 0.0f});
+  cube = cube.Translate({-0.3f, -0.3f, 0.0f});
   Manifold out = first - cube;
   CheckStrictly(out);
   EXPECT_EQ(out.Genus(), -1);
@@ -562,8 +559,7 @@ TEST(Boolean, MultiCoplanar) {
 
 TEST(Boolean, FaceUnion) {
   Manifold cubes = Manifold::Cube();
-  Manifold cube2 = cubes;
-  cubes += cube2.Translate({1, 0, 0});
+  cubes += cubes.Translate({1, 0, 0});
   EXPECT_EQ(cubes.Genus(), 0);
   ExpectMeshes(cubes, {{12, 20}});
   auto prop = cubes.GetProperties();
@@ -575,23 +571,21 @@ TEST(Boolean, FaceUnion) {
 
 TEST(Boolean, EdgeUnion) {
   Manifold cubes = Manifold::Cube();
-  Manifold cube2 = cubes;
-  cubes += cube2.Translate({1, 1, 0});
+  cubes += cubes.Translate({1, 1, 0});
   ExpectMeshes(cubes, {{8, 12}, {8, 12}});
 }
 
 TEST(Boolean, EdgeUnion2) {
   Manifold tets = Manifold::Tetrahedron();
   Manifold cube2 = tets;
-  tets.Translate({0, 0, -1});
+  tets = tets.Translate({0, 0, -1});
   tets += cube2.Translate({0, 0, 1}).Rotate(0, 0, 90);
   ExpectMeshes(tets, {{4, 4}, {4, 4}});
 }
 
 TEST(Boolean, CornerUnion) {
   Manifold cubes = Manifold::Cube();
-  Manifold cube2 = cubes;
-  cubes += cube2.Translate({1, 1, 1});
+  cubes += cubes.Translate({1, 1, 1});
   ExpectMeshes(cubes, {{8, 12}, {8, 12}});
 }
 
@@ -601,8 +595,7 @@ TEST(Boolean, CornerUnion) {
  */
 TEST(Boolean, Split) {
   Manifold cube = Manifold::Cube(glm::vec3(2.0f), true);
-  Manifold oct = Manifold::Sphere(1, 4);
-  oct.Translate(glm::vec3(0.0f, 0.0f, 1.0f));
+  Manifold oct = Manifold::Sphere(1, 4).Translate(glm::vec3(0.0f, 0.0f, 1.0f));
   std::pair<Manifold, Manifold> splits = cube.Split(oct);
   CheckStrictly(splits.first);
   CheckStrictly(splits.second);
@@ -613,8 +606,8 @@ TEST(Boolean, Split) {
 
 TEST(Boolean, SplitByPlane) {
   Manifold cube = Manifold::Cube(glm::vec3(2.0f), true);
-  cube.Translate({0.0f, 1.0f, 0.0f});
-  cube.Rotate(90.0f, 0.0f, 0.0f);
+  cube = cube.Translate({0.0f, 1.0f, 0.0f});
+  cube = cube.Rotate(90.0f, 0.0f, 0.0f);
   std::pair<Manifold, Manifold> splits =
       cube.SplitByPlane({0.0f, 0.0f, 1.0f}, 1.0f);
   CheckStrictly(splits.first);
@@ -631,9 +624,9 @@ TEST(Boolean, SplitByPlane) {
 
 TEST(Boolean, SplitByPlane60) {
   Manifold cube = Manifold::Cube(glm::vec3(2.0f), true);
-  cube.Translate({0.0f, 1.0f, 0.0f});
-  cube.Rotate(0.0f, 0.0f, -60.0f);
-  cube.Translate({2.0f, 0.0f, 0.0f});
+  cube = cube.Translate({0.0f, 1.0f, 0.0f});
+  cube = cube.Rotate(0.0f, 0.0f, -60.0f);
+  cube = cube.Translate({2.0f, 0.0f, 0.0f});
   float phi = 30.0f;
   std::pair<Manifold, Manifold> splits =
       cube.SplitByPlane({sind(phi), -cosd(phi), 0.0f}, 1.0f);
@@ -685,8 +678,7 @@ TEST(Boolean, Winding) {
 TEST(Boolean, NonIntersecting) {
   Manifold cube1 = Manifold::Cube();
   float vol1 = cube1.GetProperties().volume;
-  Manifold cube2 = cube1;
-  cube2.Scale(glm::vec3(2)).Translate({3, 0, 0});
+  Manifold cube2 = cube1.Scale(glm::vec3(2)).Translate({3, 0, 0});
   float vol2 = cube2.GetProperties().volume;
 
   EXPECT_EQ((cube1 + cube2).GetProperties().volume, vol1 + vol2);
@@ -700,12 +692,12 @@ TEST(Boolean, Precision) {
   Manifold cube3 = cube;
   float distance = 100;
   float scale = distance * kTolerance;
-  cube2.Scale(glm::vec3(scale)).Translate({distance, 0, 0});
+  cube2 = cube2.Scale(glm::vec3(scale)).Translate({distance, 0, 0});
 
   cube += cube2;
   ExpectMeshes(cube, {{8, 12}});
 
-  cube3.Scale(glm::vec3(2 * scale)).Translate({distance, 0, 0});
+  cube3 = cube3.Scale(glm::vec3(2 * scale)).Translate({distance, 0, 0});
   cube += cube3;
   ExpectMeshes(cube, {{8, 12}, {8, 12}});
 }
@@ -716,10 +708,10 @@ TEST(Boolean, Precision2) {
   Manifold cube2 = cube;
   float distance = scale * (1 - kTolerance / 2);
 
-  cube2.Translate(glm::vec3(-distance));
+  cube2 = cube2.Translate(glm::vec3(-distance));
   EXPECT_TRUE((cube ^ cube2).IsEmpty());
 
-  cube2.Translate(glm::vec3(scale * kTolerance));
+  cube2 = cube2.Translate(glm::vec3(scale * kTolerance));
   EXPECT_FALSE((cube ^ cube2).IsEmpty());
 }
 
@@ -730,8 +722,8 @@ TEST(Boolean, Precision2) {
 TEST(Boolean, Sphere) {
   Manifold sphere = Manifold::Sphere(1.0f, 12);
   Manifold sphere2 = sphere;
-  sphere2.Translate(glm::vec3(0.5));
-  sphere2.SetAsOriginal();
+  sphere2 = sphere2.Translate(glm::vec3(0.5));
+  sphere2 = sphere2.AsOriginal();
   Manifold result = sphere - sphere2;
 
   ExpectMeshes(result, {{74, 144}});
@@ -849,7 +841,7 @@ TEST(Boolean, Cylinders) {
         mat[i][j] = array[j * 4 + i];
       }
     }
-    m1 += Manifold(rod).Transform(mat);
+    m1 += rod.Transform(mat);
   }
 
   Manifold m2;
@@ -860,7 +852,7 @@ TEST(Boolean, Cylinders) {
         mat[i][j] = array[j * 4 + i];
       }
     }
-    m2 += Manifold(rod).Transform(mat);
+    m2 += rod.Transform(mat);
   }
   m1 += m2;
 
@@ -882,4 +874,4 @@ TEST(Boolean, Cubes) {
   EXPECT_NEAR(prop.surfaceArea, 9.2, 0.01);
 
   if (options.exportModels) ExportMesh("cubes.glb", result.GetMesh(), {});
-}
\ No newline at end of file
+}
diff --git a/test/samples_test.cpp b/test/samples_test.cpp
index 4004a8225..0a7374444 100644
--- a/test/samples_test.cpp
+++ b/test/samples_test.cpp
@@ -100,7 +100,7 @@ TEST(Samples, Scallop) {
     ExportMesh("scallopFacets.glb", in, options);
   }
 
-  scallop.Refine(50);
+  scallop = scallop.Refine(50);
   CheckManifold(scallop);
   auto prop = scallop.GetProperties();
   EXPECT_NEAR(prop.volume, 41.3, 0.1);
@@ -126,10 +126,9 @@ TEST(Samples, TetPuzzle) {
   Manifold puzzle = TetPuzzle(50, 0.2, 50);
   CheckManifold(puzzle);
   EXPECT_LE(puzzle.NumDegenerateTris(), 2);
-  Manifold puzzle2 = puzzle;
-  puzzle2.Rotate(0, 0, 180);
+  Manifold puzzle2 = puzzle.Rotate(0, 0, 180);
   EXPECT_TRUE((puzzle ^ puzzle2).IsEmpty());
-  puzzle.Transform(RotateUp({1, -1, -1}));
+  puzzle = puzzle.Transform(RotateUp({1, -1, -1}));
   if (options.exportModels) ExportMesh("tetPuzzle.glb", puzzle.GetMesh(), {});
 }
 
diff --git a/tools/perf_test.cpp b/tools/perf_test.cpp
index 2ad385a6b..c2ac36538 100644
--- a/tools/perf_test.cpp
+++ b/tools/perf_test.cpp
@@ -22,10 +22,10 @@ using namespace manifold;
 int main(int argc, char **argv) {
   for (int i = 0; i < 8; ++i) {
     Manifold sphere = Manifold::Sphere(1, (8 << i) * 4);
-    Manifold sphere2 = sphere;
-    sphere2.Translate(glm::vec3(0.5));
+    Manifold sphere2 = sphere.Translate(glm::vec3(0.5));
     auto start = std::chrono::high_resolution_clock::now();
     Manifold diff = sphere - sphere2;
+    diff.NumTri();
     auto end = std::chrono::high_resolution_clock::now();
     std::chrono::duration<double> elapsed = end - start;
     std::cout << "nTri = " << sphere.NumTri() << ", time = " << elapsed.count()
diff --git a/utilities/include/par.h b/utilities/include/par.h
index 0fac6dafc..0a3e55848 100644
--- a/utilities/include/par.h
+++ b/utilities/include/par.h
@@ -132,6 +132,7 @@ THRUST_DYNAMIC_BACKEND_VOID(inclusive_scan)
 THRUST_DYNAMIC_BACKEND_VOID(exclusive_scan)
 THRUST_DYNAMIC_BACKEND_VOID(uninitialized_fill)
 THRUST_DYNAMIC_BACKEND_VOID(uninitialized_copy)
+THRUST_DYNAMIC_BACKEND_VOID(copy_n)
 
 THRUST_DYNAMIC_BACKEND(all_of, bool)
 THRUST_DYNAMIC_BACKEND(is_sorted, bool)
@@ -145,6 +146,7 @@ THRUST_DYNAMIC_BACKEND(copy_if, void)
 THRUST_DYNAMIC_BACKEND(remove_if, void)
 THRUST_DYNAMIC_BACKEND(unique, void)
 THRUST_DYNAMIC_BACKEND(find, void)
+THRUST_DYNAMIC_BACKEND(find_if, void)
 THRUST_DYNAMIC_BACKEND(reduce_by_key, void)
 THRUST_DYNAMIC_BACKEND(transform_reduce, void)
 THRUST_DYNAMIC_BACKEND(lower_bound, void)
diff --git a/utilities/include/sparse.h b/utilities/include/sparse.h
index 2562e561d..d34625e1c 100644
--- a/utilities/include/sparse.h
+++ b/utilities/include/sparse.h
@@ -14,11 +14,6 @@
 
 #pragma once
 #include <math.h>
-#include <thrust/binary_search.h>
-#include <thrust/gather.h>
-#include <thrust/remove.h>
-#include <thrust/sort.h>
-#include <thrust/unique.h>
 
 #include "par.h"
 #include "structs.h"
diff --git a/utilities/include/structs.h b/utilities/include/structs.h
index abafd1e10..d27c44c66 100644
--- a/utilities/include/structs.h
+++ b/utilities/include/structs.h
@@ -415,18 +415,6 @@ struct Box {
   }
 };
 
-/**
- * Print the contents of this vector to standard output.
- */
-template <typename T>
-void Dump(const std::vector<T>& vec) {
-  std::cout << "Vec = " << std::endl;
-  for (int i = 0; i < vec.size(); ++i) {
-    std::cout << i << ", " << vec[i] << ", " << std::endl;
-  }
-  std::cout << std::endl;
-}
-
 inline std::ostream& operator<<(std::ostream& stream, const Box& box) {
   return stream << "min: " << box.min.x << ", " << box.min.y << ", "
                 << box.min.z << ", "
@@ -468,6 +456,18 @@ inline std::ostream& operator<<(std::ostream& stream, const BaryRef& ref) {
   return stream << "meshID: " << ref.meshID << ", tri: " << ref.tri
                 << ", uvw idx: " << ref.vertBary;
 }
+
+/**
+ * Print the contents of this vector to standard output.
+ */
+template <typename T>
+void Dump(const std::vector<T>& vec) {
+  std::cout << "Vec = " << std::endl;
+  for (int i = 0; i < vec.size(); ++i) {
+    std::cout << i << ", " << vec[i] << ", " << std::endl;
+  }
+  std::cout << std::endl;
+}
 /** @} */
 }  // namespace manifold