Simplify the CsgOpNodes as we build them, rather than in GetChildren (elalish#368)

* Flatten the CsgOpNodes as we build them, rather than in GetChildren

* Don't inline reused nodes in CsgNode::Boolean

* Reinstate difference -> union rewrite

* No need for opportunistic flattening anymore

* Fix typo

* Remove CsgOpNode::Impl::simplified_

* Don't overskip referenced nodes in simplification

* Avoid needless copy of shared pointers

* Disambiguate CsgOpNode::Impl::flatten_ -> forcedToLeafNodes_ and align GetChildren arg naming

* Simplified CsgOpNode::GetChildren

* Simpler (a-(b+c)) -> (a-b-c) transform in CsgOpNode::Boolean to work better w/ CsgOpNode::ToLeafNode

* Simpler build hints for large_scene_test.cpp

* [flatten] Refine rules about tree flattening + added test case

* [manifold] Propagate transforms when flattening op trees

* [flatten] move tests from manifold_test to boolean_test

Author: ochafik

extras/CMakeLists.txt

@@ -16,18 +16,20 @@ project(extras)
 
 add_executable(perfTest perf_test.cpp)
 target_link_libraries(perfTest manifold)
-
 target_compile_options(perfTest PRIVATE ${MANIFOLD_FLAGS})
 target_compile_features(perfTest PUBLIC cxx_std_17)
 
+add_executable(largeSceneTest large_scene_test.cpp)
+target_link_libraries(largeSceneTest manifold)
+target_compile_options(largeSceneTest PRIVATE ${MANIFOLD_FLAGS})
+target_compile_features(largeSceneTest PUBLIC cxx_std_17)
+
 if(BUILD_TEST_CGAL)
     add_executable(perfTestCGAL perf_test_cgal.cpp)
     find_package(CGAL REQUIRED COMPONENTS Core)
     target_compile_definitions(perfTestCGAL PRIVATE CGAL_USE_GMPXX)
-
     # target_compile_definitions(perfTestCGAL PRIVATE CGAL_DEBUG)
     target_link_libraries(perfTestCGAL manifold CGAL::CGAL CGAL::CGAL_Core)
-
     target_compile_options(perfTestCGAL PRIVATE ${MANIFOLD_FLAGS})
     target_compile_features(perfTestCGAL PUBLIC cxx_std_17)
 endif()

extras/large_scene_test.cpp

@@ -0,0 +1,54 @@
+// Copyright 2020 The Manifold Authors.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//      http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include <chrono>
+#include <iostream>
+
+#include "manifold.h"
+
+using namespace manifold;
+
+/*
+  Build & execute with the following command:
+
+  ( mkdir -p build && cd build && \
+    cmake -DCMAKE_BUILD_TYPE=Release -DMANIFOLD_PAR=TBB .. && \
+    make -j && \
+    time ./extras/largeSceneTest 50 )
+*/
+int main(int argc, char **argv) {
+  int n = 20;
+  if (argc == 2) n = atoi(argv[1]);
+
+  std::cout << "n = " << n << std::endl;
+
+  auto start = std::chrono::high_resolution_clock::now();
+  Manifold scene;
+
+  for (int i = 0; i < n; ++i) {
+    for (int j = 0; j < n; ++j) {
+      for (int k = 0; k < n; ++k) {
+        if (i == 0 && j == 0 && k == 0) continue;
+
+        Manifold sphere = Manifold::Sphere(1).Translate(glm::vec3(i, j, k));
+        scene = scene.Boolean(sphere, OpType::Add);
+      }
+    }
+  }
+  scene.NumTri();
+  auto end = std::chrono::high_resolution_clock::now();
+  std::chrono::duration<double> elapsed = end - start;
+  std::cout << "nTri = " << scene.NumTri() << ", time = " << elapsed.count()
+            << " sec" << std::endl;
+}

src/manifold/src/csg_tree.cpp

@@ -73,6 +73,19 @@ struct CheckOverlap {
 }  // namespace
 namespace manifold {
 
+std::shared_ptr<CsgNode> CsgNode::Boolean(
+    const std::shared_ptr<CsgNode> &second, OpType op) {
+  if (auto opNode = std::dynamic_pointer_cast<CsgOpNode>(second)) {
+    // "this" is not a CsgOpNode (which overrides Boolean), but if "second" is
+    // and the operation is commutative, we let it built the tree.
+    if ((op == OpType::Add || op == OpType::Intersect)) {
+      return opNode->Boolean(shared_from_this(), op);
+    }
+  }
+  std::vector<std::shared_ptr<CsgNode>> children({shared_from_this(), second});
+  return std::make_shared<CsgOpNode>(children, op);
+}
+
 std::shared_ptr<CsgNode> CsgNode::Translate(const glm::vec3 &t) const {
   glm::mat4x3 transform(1.0f);
   transform[3] += t;
@@ -227,8 +240,6 @@ CsgOpNode::CsgOpNode(const std::vector<std::shared_ptr<CsgNode>> &children,
   auto impl = impl_.GetGuard();
   impl->children_ = children;
   SetOp(op);
-  // opportunistically flatten the tree without costly evaluation
-  GetChildren(false);
 }
 
 CsgOpNode::CsgOpNode(std::vector<std::shared_ptr<CsgNode>> &&children,
@@ -237,8 +248,50 @@ CsgOpNode::CsgOpNode(std::vector<std::shared_ptr<CsgNode>> &&children,
   auto impl = impl_.GetGuard();
   impl->children_ = children;
   SetOp(op);
-  // opportunistically flatten the tree without costly evaluation
-  GetChildren(false);
+}
+
+std::shared_ptr<CsgNode> CsgOpNode::Boolean(
+    const std::shared_ptr<CsgNode> &second, OpType op) {
+  std::vector<std::shared_ptr<CsgNode>> children;
+
+  auto isReused = [](const auto &node) { return node->impl_.UseCount() > 1; };
+
+  auto copyChildren = [&](const auto &list, const glm::mat4x3 &transform) {
+    for (const auto &child : list) {
+      children.push_back(child->Transform(transform));
+    }
+  };
+
+  auto self = std::dynamic_pointer_cast<CsgOpNode>(shared_from_this());
+  assert(self);
+  if (IsOp(op) && !isReused(self)) {
+    auto impl = impl_.GetGuard();
+    copyChildren(impl->children_, transform_);
+  } else {
+    children.push_back(self);
+  }
+
+  auto secondOp = std::dynamic_pointer_cast<CsgOpNode>(second);
+  auto canInlineSecondOp = [&]() {
+    switch (op) {
+      case OpType::Add:
+      case OpType::Intersect:
+        return secondOp->IsOp(op);
+      case OpType::Subtract:
+        return secondOp->IsOp(OpType::Add);
+      default:
+        return false;
+    }
+  };
+
+  if (secondOp && canInlineSecondOp() && !isReused(secondOp)) {
+    auto secondImpl = secondOp->impl_.GetGuard();
+    copyChildren(secondImpl->children_, secondOp->transform_);
+  } else {
+    children.push_back(second);
+  }
+
+  return std::make_shared<CsgOpNode>(children, op);
 }
 
 std::shared_ptr<CsgNode> CsgOpNode::Transform(const glm::mat4x3 &m) const {
@@ -410,44 +463,24 @@ void CsgOpNode::BatchUnion() const {
 
 /**
  * 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.
+ * If forceToLeafNodes 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 {
+    bool forceToLeafNodes) const {
   auto impl = impl_.GetGuard();
-  auto &children_ = impl->children_;
-  if (children_.empty() || (impl->simplified_ && !finalize) || impl->flattened_)
-    return children_;
-  impl->simplified_ = true;
-  impl->flattened_ = finalize;
-  std::vector<std::shared_ptr<CsgNode>> newChildren;
-
-  CsgNodeType op = op_;
-  for (auto &child : children_) {
-    if (child->GetNodeType() == op && child.use_count() == 1 &&
-        std::dynamic_pointer_cast<CsgOpNode>(child)->impl_.UseCount() == 1) {
-      auto grandchildren =
-          std::dynamic_pointer_cast<CsgOpNode>(child)->GetChildren(finalize);
-      int start = children_.size();
-      for (auto &grandchild : grandchildren) {
-        newChildren.push_back(grandchild->Transform(child->GetTransform()));
-      }
-    } else {
-      if (!finalize || child->GetNodeType() == CsgNodeType::Leaf) {
-        newChildren.push_back(child);
-      } else {
-        newChildren.push_back(child->ToLeafNode());
+
+  if (forceToLeafNodes && !impl->forcedToLeafNodes_) {
+    impl->forcedToLeafNodes_ = true;
+    for (auto &child : impl->children_) {
+      if (child->GetNodeType() != CsgNodeType::Leaf) {
+        child = 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_;
+  return impl->children_;
 }
 
 void CsgOpNode::SetOp(OpType op) {
@@ -464,6 +497,19 @@ void CsgOpNode::SetOp(OpType op) {
   }
 }
 
+bool CsgOpNode::IsOp(OpType op) {
+  switch (op) {
+    case OpType::Add:
+      return op_ == CsgNodeType::Union;
+    case OpType::Subtract:
+      return op_ == CsgNodeType::Difference;
+    case OpType::Intersect:
+      return op_ == CsgNodeType::Intersection;
+    default:
+      return false;
+  }
+}
+
 glm::mat4x3 CsgOpNode::GetTransform() const { return transform_; }
 
 }  // namespace manifold

src/manifold/src/csg_tree.h

@@ -22,13 +22,16 @@ enum class CsgNodeType { Union, Intersection, Difference, Leaf };
 
 class CsgLeafNode;
 
-class CsgNode {
+class CsgNode : public std::enable_shared_from_this<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;
 
+  virtual std::shared_ptr<CsgNode> Boolean(
+      const std::shared_ptr<CsgNode> &second, OpType op);
+
   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,
@@ -68,6 +71,9 @@ class CsgOpNode final : public CsgNode {
 
   CsgOpNode(std::vector<std::shared_ptr<CsgNode>> &&children, OpType op);
 
+  std::shared_ptr<CsgNode> Boolean(const std::shared_ptr<CsgNode> &second,
+                                   OpType op) override;
+
   std::shared_ptr<CsgNode> Transform(const glm::mat4x3 &m) const override;
 
   std::shared_ptr<CsgLeafNode> ToLeafNode() const override;
@@ -79,8 +85,7 @@ class CsgOpNode final : public CsgNode {
  private:
   struct Impl {
     std::vector<std::shared_ptr<CsgNode>> children_;
-    bool simplified_ = false;
-    bool flattened_ = false;
+    bool forcedToLeafNodes_ = false;
   };
   mutable ConcurrentSharedPtr<Impl> impl_ = ConcurrentSharedPtr<Impl>(Impl{});
   CsgNodeType op_;
@@ -89,6 +94,7 @@ class CsgOpNode final : public CsgNode {
   mutable std::shared_ptr<CsgLeafNode> cache_ = nullptr;
 
   void SetOp(OpType);
+  bool IsOp(OpType op);
 
   static std::shared_ptr<Manifold::Impl> BatchBoolean(
       OpType operation,
@@ -97,7 +103,7 @@ class CsgOpNode final : public CsgNode {
   void BatchUnion() const;
 
   std::vector<std::shared_ptr<CsgNode>> &GetChildren(
-      bool finalize = true) const;
+      bool forceToLeafNodes = true) const;
 };
 
 }  // namespace manifold

src/manifold/src/manifold.cpp

@@ -585,8 +585,7 @@ Manifold Manifold::Refine(int n) const {
  * @param op The type of operation to perform.
  */
 Manifold Manifold::Boolean(const Manifold& second, OpType op) const {
-  std::vector<std::shared_ptr<CsgNode>> children({pNode_, second.pNode_});
-  return Manifold(std::make_shared<CsgOpNode>(children, op));
+  return Manifold(pNode_->Boolean(second.pNode_, op));
 }
 
 Manifold Manifold::BatchBoolean(const std::vector<Manifold>& manifolds,

test/boolean_test.cpp

@@ -589,3 +589,39 @@ TEST(Boolean, UnionDifference) {
   float blocksize = block.GetProperties().volume;
   EXPECT_NEAR(resultsize, blocksize * 2, 0.0001);
 }
+
+TEST(Boolean, BooleanVolumes) {
+  glm::mat4 m = glm::translate(glm::mat4(1.0f), glm::vec3(1.0f));
+
+  // Define solids which volumes are easy to compute w/ bit arithmetics:
+  // m1, m2, m4 are unique, non intersecting "bits" (of volume 1, 2, 4)
+  // m3 = m1 + m2
+  // m7 = m1 + m2 + m3
+  auto m1 = Manifold::Cube({1, 1, 1});
+  auto m2 = Manifold::Cube({2, 1, 1}).Transform(
+      glm::translate(glm::mat4(1.0f), glm::vec3(1.0f, 0, 0)));
+  auto m4 = Manifold::Cube({4, 1, 1}).Transform(
+      glm::translate(glm::mat4(1.0f), glm::vec3(3.0f, 0, 0)));
+  auto m3 = Manifold::Cube({3, 1, 1});
+  auto m7 = Manifold::Cube({7, 1, 1});
+
+  EXPECT_FLOAT_EQ((m1 ^ m2).GetProperties().volume, 0);
+  EXPECT_FLOAT_EQ((m1 + m2 + m4).GetProperties().volume, 7);
+  EXPECT_FLOAT_EQ((m1 + m2 - m4).GetProperties().volume, 3);
+  EXPECT_FLOAT_EQ((m1 + (m2 ^ m4)).GetProperties().volume, 1);
+  EXPECT_FLOAT_EQ((m7 ^ m4).GetProperties().volume, 4);
+  EXPECT_FLOAT_EQ((m7 ^ m3 ^ m1).GetProperties().volume, 1);
+  EXPECT_FLOAT_EQ((m7 ^ (m1 + m2)).GetProperties().volume, 3);
+  EXPECT_FLOAT_EQ((m7 - m4).GetProperties().volume, 3);
+  EXPECT_FLOAT_EQ((m7 - m4 - m2).GetProperties().volume, 1);
+  EXPECT_FLOAT_EQ((m7 - (m7 - m1)).GetProperties().volume, 1);
+  EXPECT_FLOAT_EQ((m7 - (m1 + m2)).GetProperties().volume, 4);
+}
+
+TEST(Boolean, TreeTransforms) {
+  auto a = (Manifold::Cube({1, 1, 1}) + Manifold::Cube({1, 1, 1}))
+               .Translate({1, 0, 0});
+  auto b = (Manifold::Cube({1, 1, 1}) + Manifold::Cube({1, 1, 1}));
+
+  EXPECT_FLOAT_EQ((a + b).GetProperties().volume, 2);
+}