IfcOpenShell/src/ifcgeom/iterator.cpp

#include "iterator.h"

/**
* Initialize iterator's list of tasks.
*
* Will automatically process first element, if 'defer-processing-first-element' is not set to `true`.
*
* @return Returns true if the iterator is initialized with any elements, false otherwise.
*
* @note
* - A true return value does not guarantee successful initialization of all elements.
*   Some elements may have failed to initialize. Check had_error_processing_elements()
*   to see whether there were errors during the initialization.
*
* - For non-concurrent iterators, a false return may occur if initialization of the first
*   element fails, even if subsequent elements could be initialized successfully.
*/
bool ifcopenshell::geom::iterator::initialize() {
	using std::chrono::high_resolution_clock;

	if (initialization_outcome_) {
		return *initialization_outcome_;
	}

	time_points[0] = high_resolution_clock::now();
	std::vector<ifcopenshell::geom::geometry_conversion_task> reps;
	if (num_threads_ != 1) {
		// @todo this shouldn't be necessary with properly immutable taxonomy items
		converter_->mapping()->use_caching() = false;
	}
	try {
		converter_->mapping()->get_representations(reps, filters_);
	} catch (const std::exception& e) {
		logger_.error("GEO", 50, e);
	}
	time_points[1] = high_resolution_clock::now();

	for (auto& task : reps) {
		geometry_conversion_result res;
		res.index = task.index;
		if (!settings_.get<ifcopenshell::geom::settings::NoParallelMapping>().get()) {
			res.representation = task.representation;
			res.products_2 = task.products;
		} else {
			res.item = converter_->mapping()->map(task.representation);
			if (!res.item) {
				continue;
			}
			std::transform(task.products.begin(), task.products.end(), std::back_inserter(res.products), [this, &res](const express::base& prod) {
				auto prod_item = converter_->mapping()->map(prod);
				return std::make_pair(prod, ifcopenshell::geom::taxonomy::cast<ifcopenshell::geom::taxonomy::geom_item>(prod_item)->matrix);
			});
		}
		tasks_.push_back(std::move(res));
	}

	if (settings_.get<ifcopenshell::geom::settings::NoParallelMapping>().get() && settings_.get<ifcopenshell::geom::settings::PermissiveShapeReuse>().get()) {
		std::unordered_map<
			ifcopenshell::geom::taxonomy::item::ptr,
			std::vector<std::pair<express::base, ifcopenshell::geom::taxonomy::matrix4::ptr>>> folded;

		for (auto& r : tasks_) {
			auto i = r.item;

			Eigen::Matrix4d m4 = Eigen::Matrix4d::Identity();

			while (auto col = std::dynamic_pointer_cast<ifcopenshell::geom::taxonomy::collection>(i)) {
				if (col->children.size() == 1) {
					if (col->matrix) {
						m4 *= col->matrix->ccomponents();
					}
					i = col->children[0];
				} else {
					break;
				}
			}

			for (auto& p : r.products) {
				auto pl = ifcopenshell::geom::taxonomy::matrix4::ptr(p.second->clone_());
				pl->components() *= m4;
				folded[i].push_back(
					{ p.first, pl }
				);
			}
		}

		if (folded.size() < tasks_.size()) {
			auto old_size = tasks_.size();
			tasks_.clear();
			int i = 0;
			for (auto& p : folded) {
				tasks_.emplace_back();
				tasks_.back().index = i++;
				tasks_.back().item = p.first;
				tasks_.back().products = p.second;
			}
			logger_.notice("SYS", 26, "Merged " + std::to_string(old_size) + " tasks into " + std::to_string(tasks_.size()) + " tasks due to permissive shape reuse");
		}
	}

	if (settings_.get<ifcopenshell::geom::settings::NoParallelMapping>().get()) {
		remove_offset_();
	}

	size_t num_products = 0;
	for (auto& r : tasks_) {
		num_products += !settings_.get<ifcopenshell::geom::settings::NoParallelMapping>().get() ? r.products_2.size() : r.products.size();
	}

	time_points[2] = high_resolution_clock::now();

	/*
	// What to do, map representation and product individually?
	// There needs to be two options, mapped item respecting (does that still work?), and optimized based on topology sorting.
	// Or is the sorting not necessary if we just cache?

	std::vector<taxonomy::ptr> items;
	std::map<taxonomy::ptr, taxonomy::matrix4> placements;
	std::transform(products.begin(), products.end(), std::back_inserter(items), [this, &placements](express::base p) {
	auto item = converter_->mapping()->map(p);
	// Product placements do not affect item reuse and should temporarily be swapped to identity
	if (item) {
	std::swap(placements[item], ((taxonomy::geom_ptr)item)->matrix);
	}
	return item;
	});
	items.erase(std::remove(items.begin(), items.end(), nullptr), items.end());
	std::sort(items.begin(), items.end(), taxonomy::less);
	auto it = items.begin();
	while (it < items.end()) {
	auto jt = std::upper_bound(it, items.end(), *it, taxonomy::less);
	geometry_conversion_result r;
	r.item = *it;
	std::transform(it, jt, std::back_inserter(r.products), [&r, &placements](taxonomy::ptr product_node) {
	return std::make_pair((express::base) product_node->instance, placements[product_node]);
	});
	tasks_.push_back(r);
	it = jt;
	}
	*/

	logger_.notice("SYS", 27, "Created " + boost::lexical_cast<std::string>(tasks_.size()) + " tasks for " + boost::lexical_cast<std::string>(num_products) + " products");

	if (tasks_.size() == 0) {
		logger_.warning("GEO", 51, "No representations encountered, aborting");
		initialization_outcome_ = false;
	} else if (!settings_.get<ifcopenshell::geom::settings::DeferProcessingFirstElement>().get()) {

		task_iterator_ = tasks_.begin();

		done = 0;
		total = (int)tasks_.size();

		if (num_threads_ != 1) {
			init_future_ = std::async(std::launch::async, [this]() { process_concurrently(); });

			// wait for the first element, because after init(), get() can be called.
			// so the element conversion must succeed
			initialization_outcome_ = wait_for_element();
		} else {
			initialization_outcome_ = create();
		}
	} else {
		initialization_outcome_.emplace(true);
	}

	return *initialization_outcome_;
}

void ifcopenshell::geom::iterator::flush_worker_log(ifcopenshell::geom::converter* kernel) {
	if (kernel && &kernel->logger() != &logger_) {
		logger_.append(kernel->logger());
	}
}

void ifcopenshell::geom::iterator::process_finished_rep(geometry_conversion_result* rep, ifcopenshell::geom::converter* kernel) {
	flush_worker_log(kernel);

	if (rep->elements.empty()) {
		return;
	}

	std::lock_guard<std::mutex> lk(element_ready_mutex_);

	all_processed_elements_.insert(
		all_processed_elements_.end(),
		std::make_move_iterator(rep->elements.begin()),
		std::make_move_iterator(rep->elements.end()));
	all_processed_native_elements_.insert(
		all_processed_native_elements_.end(),
		std::make_move_iterator(rep->native_elements.begin()),
		std::make_move_iterator(rep->native_elements.end()));
	rep->elements.clear();
	rep->native_elements.clear();

	if (!task_result_ptr_initialized) {
		task_result_iterator_ = all_processed_elements_.begin();
		native_task_result_iterator_ = all_processed_native_elements_.begin();
		task_result_ptr_initialized = true;
	}

	progress_ = (int)(++processed_ * 100 / tasks_.size());
}

void ifcopenshell::geom::iterator::process_concurrently() {
	size_t conc_threads = num_threads_;
	if (conc_threads > tasks_.size()) {
		conc_threads = tasks_.size();
	}

	kernel_pool.reserve(conc_threads);
	worker_loggers_.reserve(conc_threads);
	for (unsigned i = 0; i < conc_threads; ++i) {
		worker_loggers_.emplace_back(std::make_unique<logger>());
		ifcopenshell::logger& worker_logger = *worker_loggers_.back();
		worker_logger.verbosity(logger_.verbosity());
		worker_logger.output_format(logger_.output_format());
		worker_logger.print_performance_stats_on_element(logger_.print_performance_stats_on_element());
		if (worker_logger.output_format() != ifcopenshell::logger::FMT_INMEMORY) {
			worker_logger.set_output(static_cast<std::ostream*>(nullptr), static_cast<std::ostream*>(nullptr));
		}
		kernel_pool.push_back(new ifcopenshell::geom::converter(std::unique_ptr<ifcopenshell::geom::kernels::abstract_kernel>(converter_->kernel()->clone(worker_logger)), ifc_file, settings_, worker_logger));
	}

	std::vector<std::future<geometry_conversion_result*>> threadpool;

	for (auto& rep : tasks_) {
		ifcopenshell::geom::converter* K = nullptr;
		if (threadpool.size() < kernel_pool.size()) {
			K = kernel_pool[threadpool.size()];
		}

		while (threadpool.size() == conc_threads) {
			for (int i = 0; i < (int)threadpool.size(); i++) {
				auto& fu = threadpool[i];
				std::future_status status;
				status = fu.wait_for(std::chrono::seconds(0));
				if (status == std::future_status::ready) {
					process_finished_rep(fu.get(), kernel_pool[i]);

					std::swap(threadpool[i], threadpool.back());
					threadpool.pop_back();
					std::swap(kernel_pool[i], kernel_pool.back());
					std::swap(worker_loggers_[i], worker_loggers_.back());
					K = kernel_pool.back();
					break;
				} // if
			}   // for
		}	 // while

		std::future<geometry_conversion_result*> fu = std::async(
			std::launch::async, [this](
				ifcopenshell::geom::converter* kernel,
				ifcopenshell::geom::settings settings,
				geometry_conversion_result* rep) {
			// Catch exceptions to be safe from freezing the iterator.
			try {
				this->create_element_(kernel, settings, rep);
			} catch (const std::exception& e) {
				kernel->logger().error("GEO", 52,
					std::string("Exception '") + e.what() +
					std::string("' occurred while iterator was creating a shape: "),
					rep->item->instance
				);
				had_error_processing_elements_ = true;
			} catch (...) {
				kernel->logger().error("GEO", 53,
					"Unknown exception occurred while iteartor was creating a shape: ",
					rep->item->instance
				);
				had_error_processing_elements_ = true;
			}
			return rep;
		},
			K,
			std::ref(settings_),
			&rep);

		if (terminating_) {
			break;
		}

		threadpool.emplace_back(std::move(fu));
	}

	for (size_t i = 0; i < threadpool.size(); ++i) {
		process_finished_rep(threadpool[i].get(), kernel_pool[i]);
	}

	finished_ = true;

	logger_.set_product(std::optional<express::base>{});

	if (!terminating_) {
		logger_.status("\rDone creating geometry (" + boost::lexical_cast<std::string>(all_processed_elements_.size()) +
			" objects)								");
	}
}

/// Computes model's bounding box (bounds_min and bounds_max).
/// @note Can take several minutes for large files.
void ifcopenshell::geom::iterator::compute_bounds(bool with_geometry)
{
	for (int i = 0; i < 3; ++i) {
		bounds_min_.components()(i) = std::numeric_limits<double>::infinity();
		bounds_max_.components()(i) = -std::numeric_limits<double>::infinity();
	}

	if (with_geometry) {
		size_t num_created = 0;
		do {
			auto geom_object = get();
			const ifcopenshell::geom::triangulation_element* o = static_cast<const ifcopenshell::geom::triangulation_element*>(geom_object.get());
			const ifcopenshell::geom::triangulation& mesh = o->geometry();
			auto mat = o->transformation().data()->ccomponents();
			Eigen::Vector4d vec, transformed;

			for (typename std::vector<double>::const_iterator it = mesh.verts().begin(); it != mesh.verts().end();) {
				const double& x = *(it++);
				const double& y = *(it++);
				const double& z = *(it++);
				vec << x, y, z, 1.;
				transformed = mat * vec;

				for (int i = 0; i < 3; ++i) {
					bounds_min_.components()(i) = std::min(bounds_min_.components()(i), transformed(i));
					bounds_max_.components()(i) = std::max(bounds_max_.components()(i), transformed(i));
				}
			}
		} while (++num_created, next());
	} else {
		std::vector<ifcopenshell::geom::geometry_conversion_task> reps;
		converter_->mapping()->get_representations(reps, filters_);

		std::vector<express::base> products;
		for (auto& r : reps) {
			std::copy(r.products.begin(), r.products.end(), std::back_inserter(products));
		}

		for (auto& product : products) {
			auto prod_item = converter_->mapping()->map(product);
			auto vec = ifcopenshell::geom::taxonomy::cast<ifcopenshell::geom::taxonomy::geom_item>(prod_item)->matrix->translation_part();

			for (int i = 0; i < 3; ++i) {
				bounds_min_.components()(i) = std::min(bounds_min_.components()(i), vec(i));
				bounds_max_.components()(i) = std::max(bounds_max_.components()(i), vec(i));
			}
		}
	}
}

express::base ifcopenshell::geom::iterator::create_shape_model_for_next_entity() {
	geometry_conversion_result* task = nullptr;
	for (; task_iterator_ < tasks_.end();) {
		task = &*task_iterator_++;
		create_element_(converter_, settings_, task);
		if (task->elements.empty()) {
			task = nullptr;
		} else {
			break;
		}
	}
	if (task) {
		process_finished_rep(task);
        return task->item->instance;
	} else {
        return express::base{};
	}
}

void ifcopenshell::geom::iterator::create_element_(ifcopenshell::geom::converter* kernel, ifcopenshell::geom::settings settings, geometry_conversion_result* rep)
{
	ifcopenshell::logger& kernel_logger = kernel->logger();

	if (!settings_.get<ifcopenshell::geom::settings::NoParallelMapping>().get()) {
		rep->item = kernel->mapping()->map(rep->representation);
		if (!rep->item) {
			return;
		}
		std::transform(rep->products_2.begin(), rep->products_2.end(), std::back_inserter(rep->products), [this, &rep, kernel](const express::base& prod) {
			auto prod_item = kernel->mapping()->map(prod);
			return std::make_pair(prod, ifcopenshell::geom::taxonomy::cast<ifcopenshell::geom::taxonomy::geom_item>(prod_item)->matrix);
		});
	} else {
	}

	auto product_node = rep->products.front();
	const express::base product = product_node.first;
	const auto& place = product_node.second;

	kernel_logger.set_product(product);

	std::unique_ptr<ifcopenshell::geom::native_element> brep(static_cast<ifcopenshell::geom::native_element*>(create_processed_element_([kernel, settings, product, place, rep]() {
		return kernel->create_brep_for_representation_and_product(rep->item, product, place);
	})));

	if (!brep) {
        kernel_logger.set_product(std::optional<express::base>{});
		return;
	}

	auto* brep_for_reuse = brep.get();
	std::unique_ptr<ifcopenshell::geom::element> elem;
	if (settings.get<ifcopenshell::geom::settings::IteratorOutput>().get() == ifcopenshell::geom::settings::NATIVE) {
		elem = std::move(brep);
	} else {
		elem = process_based_on_settings(settings, brep.get(), kernel_logger);
	}
	if (!elem) {
        kernel_logger.set_product(std::optional<express::base>{});
		return;
	}

	rep->native_elements.push_back(std::move(brep));
	rep->elements.push_back(std::move(elem));

	for (auto it = rep->products.begin() + 1; it != rep->products.end(); ++it) {
		const auto& p = *it;
		const express::base product2 = p.first;
		const auto& place2 = p.second;

		kernel_logger.set_product(product2);

		std::unique_ptr<ifcopenshell::geom::native_element> brep2(static_cast<ifcopenshell::geom::native_element*>(create_processed_element_([kernel, settings, product2, place2, brep_for_reuse]() {
			return kernel->create_brep_for_processed_representation(product2, place2, brep_for_reuse);
		})));
		if (brep2) {
			std::unique_ptr<ifcopenshell::geom::element> elem2;
			if (settings.get<ifcopenshell::geom::settings::IteratorOutput>().get() == ifcopenshell::geom::settings::NATIVE) {
				elem2 = std::move(brep2);
			} else {
				elem2 = process_based_on_settings(
					settings,
					brep2.get(),
					kernel_logger,
					dynamic_cast<ifcopenshell::geom::triangulation_element*>(rep->elements.front().get()));
			}
			if (elem2) {
				rep->native_elements.push_back(std::move(brep2));
				rep->elements.push_back(std::move(elem2));
			}
		}
	}

	kernel_logger.set_product(std::optional<express::base>{});
}

std::unique_ptr<ifcopenshell::geom::element> ifcopenshell::geom::iterator::process_based_on_settings(ifcopenshell::geom::settings settings, ifcopenshell::geom::native_element* elem, ifcopenshell::logger& logger, ifcopenshell::geom::triangulation_element* previous)
{
	if (settings.get<ifcopenshell::geom::settings::IteratorOutput>().get() == ifcopenshell::geom::settings::SERIALIZED) {
		try {
			return std::make_unique<ifcopenshell::geom::serialized_element>(*elem);
		} catch (...) {
			logger.message(ifcopenshell::logger::LOG_ERROR, "GEO", 54, "Getting a serialized element from model failed.");
			return nullptr;
		}
	} else if (settings.get<ifcopenshell::geom::settings::IteratorOutput>().get() == ifcopenshell::geom::settings::TRIANGULATED) {
		try {
			if (!previous) {
				return std::make_unique<triangulation_element>(*elem);
			} else {
				return std::make_unique<triangulation_element>(*elem, previous->geometry_pointer());
			}
		} catch (...) {
			logger.message(ifcopenshell::logger::LOG_ERROR, "GEO", 55, "Getting a triangulation element from model failed.");
			return nullptr;
		}
	} else {
		throw std::runtime_error("native iterator output must be moved directly");
	}
}

bool ifcopenshell::geom::iterator::wait_for_element() {
	while (true) {
		size_t s;
		{
			std::lock_guard<std::mutex> lk(element_ready_mutex_);
			s = all_processed_elements_.size();
		}
		if (s > async_elements_returned_) {
			++async_elements_returned_;
			return true;
		} else if (finished_) {
			return false;
		} else {
			std::this_thread::sleep_for(std::chrono::milliseconds(10));
		}
	}
}

void ifcopenshell::geom::iterator::log_timepoints() const {
	using std::chrono::high_resolution_clock;
	using std::chrono::duration;
	using namespace std::string_literals;

	std::array<std::string, 3> labels = {
		"Initializing mapping"s,
		"Performing mapping"s,
		"Geometry interpretation"s
	};

	for (auto it = time_points.begin() + 1; it != time_points.end(); ++it) {
		auto jt = it - 1;
		duration<double, std::milli> ms_double = (*it) - (*jt);
		logger_.notice("SYS", 28, labels[std::distance(time_points.begin(), jt)] + " took " + std::to_string(ms_double.count()) + "ms");
	}
}

void ifcopenshell::geom::iterator::validate_iterator_state() const {
	if (!initialization_outcome_) {
		throw std::runtime_error("iterator not initialized");
	}

	// Causes:
	// - iterator was initialized but there were no elements to process
	// - iterator was initialized but 'defer-processing-first-element' setting is enabled
	// and some element should be processed manually first
	if (!task_result_ptr_initialized) {
		throw std::runtime_error("No elements processed");
	}

	if (task_result_ptr_exhausted) {
		throw std::runtime_error("iterator is exhausted");
	}
}

/// Moves to the next shape representation, create its geometry, and returns the associated product.
/// Use get() to retrieve the created geometry.
express::base ifcopenshell::geom::iterator::next() {
	using std::chrono::high_resolution_clock;
	validate_iterator_state();

	{
		std::lock_guard<std::mutex> lock(element_ready_mutex_);
		task_result_iterator_->reset();
		native_task_result_iterator_->reset();
	}

	if (num_threads_ != 1) {
		if (!wait_for_element()) {
			logger_.set_product(std::optional<express::base>{});
			time_points[3] = high_resolution_clock::now();
			log_timepoints();
			task_result_ptr_exhausted = true;
            return express::base{};
		}

		std::lock_guard<std::mutex> lock(element_ready_mutex_);
		task_result_iterator_++;
		native_task_result_iterator_++;

		return task_result_iterator_->get()->product();
	} else {
		// Increment the iterator over the list of products using the current
		// shape representation
		if (task_result_iterator_ == --all_processed_elements_.end()) {
			if (!create()) {
				logger_.set_product(std::optional<express::base>{});
				time_points[3] = high_resolution_clock::now();
				log_timepoints();
				task_result_ptr_exhausted = true;
                return express::base{};
			}
		}

		std::lock_guard<std::mutex> lock(element_ready_mutex_);
		task_result_iterator_++;
		native_task_result_iterator_++;

		return task_result_iterator_->get()->product();
	}
}

/// Gets the representation of the current geometrical entity.
std::unique_ptr<ifcopenshell::geom::element> ifcopenshell::geom::iterator::get()
{
	validate_iterator_state();

	std::unique_ptr<ifcopenshell::geom::element> ret;
	{
		std::lock_guard<std::mutex> lock(element_ready_mutex_);
		ret = std::move(*task_result_iterator_);
		if (!ret) {
			throw std::runtime_error("current element has already been retrieved");
		}
	}

	// If we want to organize the element considering their hierarchy
	if (settings_.get<ifcopenshell::geom::settings::UseElementHierarchy>().get()) {
		// We are going to build a vector with the element parents.
		// First, create the parent vector
		std::vector<std::unique_ptr<ifcopenshell::geom::element>> parents;

		// if the element has a parent
		if (ret->parent_id() != -1) {
			ifcopenshell::geom::element* parent_object = NULL;
			bool hasParent = true;

			// get the parent
			try {
				auto parent = get_object(ret->parent_id());
				parent_object = parent.get();
				parents.insert(parents.begin(), std::move(parent));
			} catch (const std::exception& e) {
				logger_.error("GEO", 56, e);
				hasParent = false;
			}

			// Add the previously found parent to the vector
			// We need to find all the parents
			while (parent_object != NULL && hasParent && parent_object->parent_id() != -1) {
				// Find the next parent
                auto pid = parent_object->parent_id();
                auto ifc_product = ifc_file->instance_by_id(pid);
				if (ifc_product.declaration().name() == "IfcProject") {
                    hasParent = false;
                } else {
					try {
						auto parent = get_object(pid);
						parent_object = parent.get();
						parents.insert(parents.begin(), std::move(parent));
					} catch (const std::exception& e) {
						logger_.error("GEO", 57, e);
						hasParent = false;
					}
                }
				
				// Add the previously found parent to the vector
				hasParent = hasParent && parent_object->parent_id() != -1;
			}

			// when done push the parent list in the element object
			ret->set_parents(std::move(parents));
		}
	}

	return ret;
}

std::unique_ptr<ifcopenshell::geom::element> ifcopenshell::geom::iterator::get_object(int id) {
	ifcopenshell::geom::taxonomy::matrix4::ptr m4;
	int parent_id = -1;
	std::string instance_type, product_name, product_guid;
    express::base ifc_product;

	try {
        ifc_product = ifc_file->instance_by_id(id);
		instance_type = ifc_product.declaration().name();

		if (ifc_product.declaration().is("IfcRoot")) {
			product_guid = ifc_product.as<express::entity>().get_value<std::string>("GlobalId");
            product_name = ifc_product.as<express::entity>().get_value<std::string>("Name", "");
		}

		auto parent_object = converter_->mapping()->get_decomposing_entity(ifc_product);
		if (parent_object) {
			parent_id = parent_object.id();
		}

		// fails in case of IfcProject
		auto mapped = converter_->mapping()->map(ifc_product);
		auto casted = mapped ? ifcopenshell::geom::taxonomy::dcast<ifcopenshell::geom::taxonomy::geom_item>(mapped) : nullptr;

		if (casted) {
			m4 = casted->matrix;
		}
	} catch (const std::exception& e) {
		logger_.error("GEO", 58, e);
	} catch (...) {
		logger_.error("GEO", 59, "Unknown error returning product");
	}

	return std::make_unique<element>(settings_, id, parent_id, product_name, instance_type, product_guid, "", m4, ifc_product.as<express::entity>());
}

express::base ifcopenshell::geom::iterator::create() {
    express::base product;
	try {
		product = create_shape_model_for_next_entity();
	} catch (const std::exception& e) {
		logger_.error("GEO", 60, e);
		had_error_processing_elements_ = true;
	} catch (...) {
		logger_.error("GEO", 61, "Unknown error creating geometry");
		had_error_processing_elements_ = true;
	}
	return product;
}

ifcopenshell::geom::taxonomy::direction3::ptr ifcopenshell::geom::iterator::remove_offset_() {
	
	using namespace ifcopenshell::geom::taxonomy;
	
	if (!settings_.get<ifcopenshell::geom::settings::MaxOffset>().has()) {
		return nullptr;
	}
	
	if (!settings_.get<ifcopenshell::geom::settings::NoParallelMapping>().get()) {
		throw std::runtime_error("remove_offset() can only be called with defer-processing-first-element and no-parallel-mapping settings");
	}

	auto collect_offset = [&](const item::ptr& itm, const std::vector<std::pair<express::base, matrix4::ptr>>& pr) -> std::pair<double, Eigen::Vector3d> {
		std::function<std::pair<double, Eigen::Vector3d>(const item::ptr&, Eigen::Matrix4d)> traverse;
		traverse = [&](const item::ptr& node, Eigen::Matrix4d m4) -> std::pair<double, Eigen::Vector3d> {
			if (auto shl = std::dynamic_pointer_cast<shell>(node)) {
				auto p = shl->centroid();
				Eigen::Vector4d v;
				v << p->components()(0), p->components()(1), p->components()(2), 1.0;
				Eigen::Vector3d translation_part = (m4 * v).head<3>();
				double translation_amnt = translation_part.norm();
				if (translation_amnt > settings_.get<ifcopenshell::geom::settings::MaxOffset>().get()) {
					return { translation_amnt, translation_part };
				} else {
					return { 0.0, Eigen::Vector3d::Zero() };
				}
			} else {
				if (auto gi = std::dynamic_pointer_cast<geom_item>(node)) {
					if (gi->matrix) {
						m4 = m4 * gi->matrix->ccomponents();
					}
				}
				Eigen::Vector3d translation_part = m4.block<3, 1>(0, 3);
				double translation_amnt = translation_part.norm();
				if (translation_amnt > settings_.get<ifcopenshell::geom::settings::MaxOffset>().get()) {
					return { translation_amnt, translation_part };
				} else if (auto col = std::dynamic_pointer_cast<collection>(node)) {
					std::vector<std::pair<double, Eigen::Vector3d>> child_transforms;
					for (const auto& child : col->children) {
						child_transforms.push_back(traverse(child, m4));
					}
					if (!child_transforms.empty()) {
						return *std::max_element(child_transforms.begin(), child_transforms.end(),
							[](const auto& a, const auto& b) { return a.first < b.first; });
					}
				}
				return { 0.0, Eigen::Vector3d::Zero() };
			}
		};

		Eigen::Matrix4d m4 = Eigen::Matrix4d::Identity();
		if (pr.size() == 1 && pr[0].second) {
			m4 = pr[0].second->ccomponents();
		}
		return traverse(itm, m4);
	};

	Eigen::Vector3d vec;

	if (settings_.get<ifcopenshell::geom::settings::ApplyOffset>().has()) {
		auto vs = settings_.get<ifcopenshell::geom::settings::ApplyOffset>().get();
		if (vs.size() != 3) {
			throw std::runtime_error("ApplyOffset setting must be a vector of size 3");
		}
		vec = Eigen::Vector3d(vs[0], vs[1], vs[2]);
	} else {
		// Collect all norms and vectors
		std::vector<double> norms;
		std::vector<Eigen::Vector3d> vectors;
		for (const auto& task : tasks_) {
			auto result = collect_offset(task.item, task.products);
			norms.push_back(result.first);
			vectors.push_back(result.second);
		}

		// Find the median norm index
		std::vector<double> sorted_norms = norms;
		std::nth_element(sorted_norms.begin(), sorted_norms.begin() + sorted_norms.size() / 2, sorted_norms.end());
		double median = sorted_norms[sorted_norms.size() / 2];
		auto median_it = std::find(norms.begin(), norms.end(), median);
		size_t median_index = std::distance(norms.begin(), median_it);

		if (median_index >= vectors.size()) {
			return nullptr;
		}

		vec = -vectors[median_index];
	}

	Eigen::Matrix4d translation_matrix = Eigen::Matrix4d::Identity();
	translation_matrix.block<3, 1>(0, 3) = vec;

	auto remove_offset = [&](const item::ptr& itm, const std::vector<std::pair<express::base, matrix4::ptr>>& pr) -> bool {
		std::function<bool(const item::ptr&, Eigen::Matrix4d)> traverse;
		traverse = [&](const item::ptr& node, Eigen::Matrix4d m4) -> bool {
			if (auto shl = std::dynamic_pointer_cast<shell>(node)) {
				auto p = shl->centroid();
				Eigen::Vector4d v;
				v << p->components()(0), p->components()(1), p->components()(2), 1.0;
				Eigen::Vector3d translation_part = (m4 * v).head<3>();
				double translation_amnt = translation_part.norm();
				if (translation_amnt > settings_.get<ifcopenshell::geom::settings::MaxOffset>().get()) {
					shl->matrix = make<matrix4>(translation_matrix);
				}
				return true;
			} else {
				auto m4b = m4;
				if (auto gi = std::dynamic_pointer_cast<geom_item>(node)) {
					if (gi->matrix) {
						m4b = m4 * gi->matrix->ccomponents();
					}
					Eigen::Vector3d translation_part = m4b.block<3, 1>(0, 3);
					double translation_amnt = translation_part.norm();
					if (translation_amnt > settings_.get<ifcopenshell::geom::settings::MaxOffset>().get()) {
						auto inverted_rot_scale3 = m4.block<3, 3>(0, 0).inverse();
						Eigen::Matrix4d inverted_rot_scale = Eigen::Matrix4d::Identity();
						inverted_rot_scale.block<3, 3>(0, 0) = inverted_rot_scale3;
						if (!gi->matrix) {
							gi->matrix = make<matrix4>();
						}
						gi->matrix->components() = (inverted_rot_scale * translation_matrix) * gi->matrix->ccomponents();
						return true;
					}
				}
				bool b = true;
				if (auto col = std::dynamic_pointer_cast<collection>(node)) {
					for (const auto& child : col->children) {
						if (!traverse(child, m4b)) {
							b = false;
						}
					}
				}
				return b;
			}
		};

		Eigen::Matrix4d m4 = Eigen::Matrix4d::Identity();
		if (pr.size() == 1 && pr[0].second) {
			m4 = pr[0].second->ccomponents();
		}
		return traverse(itm, m4);
	};

	size_t num_offset_applied = 0;
	for (auto& task : tasks_) {
		bool all_applied = true;
		for (auto& p : task.products) {
			auto bb = p.second->components().block<3, 1>(0, 3);
			double translation_amnt = bb.norm();
			if (translation_amnt > settings_.get<ifcopenshell::geom::settings::MaxOffset>().get()) {
				// block has an underlying mutable ref to the matrix
				bb += vec;
			} else {
				all_applied = false;
			}
		}
		if (all_applied) {
			num_offset_applied += 1;
			continue;
		}
		if (remove_offset(task.item, task.products)) {
			num_offset_applied += 1;
		}
	}

	logger_.notice("SYS", 29, "Removed large offsets within " + std::to_string(num_offset_applied) + " products");
	logger_.notice("SYS", 30, "Offset applied (" + std::to_string(vec(0)) + "," + std::to_string(vec(1)) + "," + std::to_string(vec(2)) + ")");

	return make<direction3>(vec);
}

ifcopenshell::geom::iterator::~iterator() {
	if (num_threads_ != 1) {
		terminating_ = true;

		if (init_future_.valid()) {
			init_future_.wait();
		}
	}

	for (auto& k : kernel_pool) {
		flush_worker_log(k);
		delete k;
	}

	delete converter_;
}