Experiment w/ fast-union heuristic to break local overlapping barriers

benchmark results.txt

# Context: https://github.com/openscad/openscad/pull/3636

Notice the N=3 vs. N=5 below. All test on a MacBook Pro w/ 2.6GHz Core i7.

# Baseline
time openscad spheres_single_for_loop.scad  -o out.stl -DN=3 -Doverlap=true # 52s
time openscad spheres_nested_for_loops.scad -o out.stl -DN=3 -Doverlap=true # 45s
time openscad spheres_single_for_loop.scad  -o out.stl -DN=5 -Doverlap=true # 3min32
time openscad spheres_nested_for_loops.scad -o out.stl -DN=5 -Doverlap=true # 2min41

# fast-union before this commit
time openscad spheres_single_for_loop.scad  -o out.stl --enable=fast-union -DN=3 -Doverlap=true # 48s
time openscad spheres_nested_for_loops.scad -o out.stl --enable=fast-union -DN=3 -Doverlap=true # 47s
time openscad spheres_single_for_loop.scad  -o out.stl --enable=fast-union -DN=5 -Doverlap=true # 2m38
time openscad spheres_nested_for_loops.scad -o out.stl --enable=fast-union -DN=5 -Doverlap=true # 2m13

# fast-union after this commit
time openscad spheres_single_for_loop.scad  -o out.stl --enable=fast-union -DN=3 -Doverlap=true # 35s
time openscad spheres_nested_for_loops.scad -o out.stl --enable=fast-union -DN=3 -Doverlap=true # 36s
time openscad spheres_single_for_loop.scad  -o out.stl --enable=fast-union -DN=5 -Doverlap=true # 1m42s
time openscad spheres_nested_for_loops.scad -o out.stl --enable=fast-union -DN=5 -Doverlap=true # 1m46

spheres_single_for_loop.scad:
  N = 3;
  overlap = false;

  union() {
    for (i=[0:N-1], j=[0:N-1])
      translate([i,j,0])
        sphere(d=(overlap ? 1.1 : 0.9), $fn=50);
  }

spheres_nested_for_loops.scad:
  N = 3;
  overlap = false;

  union() {
    for (i=[0:N-1]) translate([i,0,0])
    for (j=[0:N-1]) translate([0,j,0])
      sphere(d=(overlap ? 1.1 : 0.9), $fn=50);
  }

cgalutils-applyops.cc

// Copyright 2021 Google LLC.
// SPDX-License-Identifier: Apache-2.0

// Edits relative to https://github.com/openscad/openscad/pull/3636 on 20210130:

	shared_ptr<const LazyGeometry> applyUnion3DFast(
		Geometry::Geometries::iterator chbegin, Geometry::Geometries::iterator chend,
		const Tree* tree)
	{
		try {
			typedef shared_ptr<const LazyGeometry> QueueConstItem;
			struct QueueItemGreater {
				// stable sort for priority_queue by facets, then progress mark
				bool operator()(const QueueConstItem &lhs, const QueueConstItem& rhs) const
				{
					size_t l = lhs->getNumberOfFacets();
					size_t r = rhs->getNumberOfFacets();
					return (l > r) || (l == r && lhs->progress_mark > rhs->progress_mark);
				}
			};

			std::vector<shared_ptr<const LazyGeometry>> geometries;

			for (auto it = chbegin; it != chend; ++it) {
				auto chgeom = it->second;
				auto &chnode = it->first;
				if (!dynamic_pointer_cast<const PolySet>(chgeom) &&
					!dynamic_pointer_cast<const CGAL_Nef_polyhedron>(chgeom)) {
					continue;
				}
				if (chgeom && !chgeom->isEmpty()) {
					auto progress_mark = chnode ? chnode->progress_mark : -1;
					auto location = chnode && chnode->modinst ? chnode->modinst->location() : Location::NONE;
					auto cacheKey = chnode && tree ? tree->getIdString(*chnode) : "";
					geometries.push_back(make_shared<LazyGeometry>(chgeom, location, cacheKey, progress_mark));
				}
			}

			// This does wonders on chainmail as it allows fast-union of sequences
			// of locally-overlapping geometries with lots of global disjoin groups
			// of geometries.
			LazyGeometry::reduceFastUnions(geometries);

			// This build the queue in linear time.
			std::priority_queue<QueueConstItem, std::vector<QueueConstItem>, QueueItemGreater> q(geometries.begin(), geometries.end());

			progress_tick();
			while (q.size() > 1) {
				auto p1 = q.top();
				q.pop();
				assert(!q.empty());
				auto p2 = q.top();
				q.pop();
				q.emplace(make_shared<LazyGeometry>(*p1 + *p2));
				progress_tick();
			}
			assert(q.size() <= 1);
			return q.empty() ? nullptr : q.top();
		}
		catch (const CGAL::Failure_exception &e) {
			LOG(message_group::Error, Location::NONE, "", "CGAL error in CGALUtils::applyUnion3DFast: %1$s", e.what());
		}
		return nullptr;
	}

lazy_geometry.cc

// Copyright 2021 Google LLC.
// SPDX-License-Identifier: Apache-2.0
//
// This goes in lazy_geometry.cc
//

// This bugdet controls the overlap detection logic.
// We do up to this amound of unsuccessful tests to find "easy" unions for each
// geometry being unioned. Since we don't end up testing all the other
// geometries, it's not a quadratic search (even if the constant is big). But
// that's not to say it's linear either.
#define FAST_JOIN_EXPLORATION_BUDGET 10

/*! A cheap heuristic to find "easy" fast union (of disjoint geometries).
 * This allows fast-union of large sub-parts of patterns where neighbouring
 * tiles overlap, like chaimail or linked chains.
 *
 * TODO(ochafik): Update logic to fast-union eligible solids in the
 * smaller-edge-first order as for general unions.
 */
void LazyGeometry::reduceFastUnions(std::vector<shared_ptr<const LazyGeometry>> &geometries)
{
	std::list<shared_ptr<const LazyGeometry>> list;
	for (auto geom : geometries) list.push_back(geom);

	for (auto it1 = list.begin(); it1 != list.end(); it1++) {
		auto it2 = it1;
		it2++;

		auto remaining_budget = FAST_JOIN_EXPLORATION_BUDGET;
		while (it2 != list.end() && remaining_budget > 0) {
			if (!(*it1)->getBoundingBoxes()->intersects(*(*it2)->getBoundingBoxes())) {
				// Found an easy union!
				*it1 = make_shared<const LazyGeometry>(**it1 + **it2);
				it2 = list.erase(it2);
				remaining_budget = FAST_JOIN_EXPLORATION_BUDGET;
			}
			else {
				remaining_budget--;
				it2++;
			}
		}

		if (remaining_budget == 0) {
			LOG(message_group::Echo, Location::NONE, "",
					"Exhausted search budget for fast-union pairs finding");
		}
	}

	geometries.clear();
	geometries.insert(geometries.end(), list.begin(), list.end());
}