diff --git a/dUgcServer/CMakeLists.txt b/dUgcServer/CMakeLists.txt
index 0f80c9dce..7bfdf1c93 100644
--- a/dUgcServer/CMakeLists.txt
+++ b/dUgcServer/CMakeLists.txt
@@ -1,7 +1,7 @@
 # The UGC server (docs/UgcServer.md), in folders by role:
 #   Bricks/      the client's brick data: primitives, Materials.xml, LU Toolbox's palette, the CDClient's module data
 #   Model/       player models built from LXFML, car and rocket module combinations, glitter
-#   Render/      the software rasterizer, ambient occlusion, hidden-face removal, icon parameters and pose
+#   Render/      the software rasterizer, ambient occlusion, hidden-face removal, rays, icon parameters and pose
 #   Formats/     the files written for the client and the dashboard (NIF, DDS, PNG, sd0)
 #   Processing/  the job queue and workers, the CPU budget, stored files
 #   UgcServer.cpp  the server: settings, the master link, HTTP routes
@@ -20,13 +20,14 @@ set(DUGC_SOURCES
 	"Render/UgcHsr.cpp"
 	"Render/UgcIconParams.cpp"
 	"Render/UgcIconPose.cpp"
+	"Render/UgcRays.cpp"
 	"Render/UgcRender.cpp"
 )
 
 add_library(dUgc STATIC ${DUGC_SOURCES})
 target_include_directories(dUgc PUBLIC "." "Bricks" "Formats" "Model" "Processing" "Render"
 	"${PROJECT_SOURCE_DIR}/thirdparty/MD5" "${PROJECT_SOURCE_DIR}/thirdparty/nlohmann")
-target_link_libraries(dUgc PUBLIC dCommon glm::glm tinyxml2 MD5)
+target_link_libraries(dUgc PUBLIC dCommon glm::glm tinyxml2 MD5 PRIVATE embree)
 
 add_executable(UgcServer "UgcServer.cpp" "Processing/UgcProcessor.cpp" "Bricks/UgcCdClient.cpp")
 
@@ -51,5 +52,5 @@ target_link_libraries(UgcServer ${COMMON_LIBRARIES} dWeb dServer dUgc)
 
 # Rendering is far too slow unoptimized (minutes a model instead of seconds), so it is optimized in every build type
 if(NOT MSVC)
-	set_source_files_properties("Render/UgcRender.cpp" "Render/UgcHsr.cpp" "Model/UgcModel.cpp" "Bricks/UgcPalette.cpp" "Model/UgcGlitter.cpp" PROPERTIES COMPILE_OPTIONS "-O2")
+	set_source_files_properties("Render/UgcRender.cpp" "Render/UgcHsr.cpp" "Render/UgcRays.cpp" "Model/UgcModel.cpp" "Bricks/UgcPalette.cpp" "Model/UgcGlitter.cpp" PROPERTIES COMPILE_OPTIONS "-O2")
 endif()
diff --git a/dUgcServer/Render/UgcHsr.cpp b/dUgcServer/Render/UgcHsr.cpp
index b11cd87c0..84c3d2620 100644
--- a/dUgcServer/Render/UgcHsr.cpp
+++ b/dUgcServer/Render/UgcHsr.cpp
@@ -4,8 +4,7 @@
 #include <array>
 #include <cmath>
 #include <limits>
-#include <numeric>
-
+#include "UgcRays.h"
 #include "UgcThrottle.h"
 
 namespace {
@@ -298,219 +297,6 @@ namespace {
 		float m_Alpha{};
 	};
 
-	struct Hit {
-		float t{ INF };
-		uint32_t triangle{ UINT32_MAX };
-		float u{}, v{}; // weights of the triangle's second and third vertex
-	};
-
-	/**
-	 * A bounding volume hierarchy (binned surface area heuristic) over a mesh's triangles, for the paths' rays: the
-	 * nearest hit. A ray never hits the triangle it leaves (`skip`), as in Cycles.
-	 */
-	class Bvh {
-	public:
-		explicit Bvh(const UgcModel::Mesh& mesh) {
-			const size_t count = mesh.TriangleCount();
-			std::vector<uint32_t> order(count);
-			std::iota(order.begin(), order.end(), 0u);
-			std::vector<glm::vec3> lo(count), hi(count), centre(count);
-			for (size_t t = 0; t < count; t++) {
-				const auto& a = mesh.positions[mesh.indices[t * 3]];
-				const auto& b = mesh.positions[mesh.indices[t * 3 + 1]];
-				const auto& c = mesh.positions[mesh.indices[t * 3 + 2]];
-				lo[t] = glm::min(a, glm::min(b, c));
-				hi[t] = glm::max(a, glm::max(b, c));
-				centre[t] = (lo[t] + hi[t]) * 0.5f;
-			}
-			if (count > 0) Build(order, lo, hi, centre);
-			m_Triangles.reserve(count);
-			for (const auto t : order) {
-				const auto& a = mesh.positions[mesh.indices[t * 3]];
-				m_Triangles.push_back({ a, mesh.positions[mesh.indices[t * 3 + 1]] - a, mesh.positions[mesh.indices[t * 3 + 2]] - a, t });
-			}
-		}
-
-		// The nearest triangle along the ray (unit direction) before `maxT`
-		Hit Closest(const glm::vec3& origin, const glm::vec3& direction, uint32_t skip, float maxT = INF) const {
-			Hit hit;
-			hit.t = maxT;
-			if (m_Nodes.empty()) return hit;
-			const auto inverse = Inverse(direction);
-			uint32_t stack[128];
-			int top = 0;
-			uint32_t index = 0;
-			while (true) {
-				const auto& node = m_Nodes[index];
-				if (node.count > 0) {
-					for (uint32_t i = node.first; i < node.first + node.count; i++) Intersect(m_Triangles[i], origin, direction, skip, hit);
-				} else {
-					const uint32_t left = node.first, right = node.first + 1;
-					const float tl = Enter(m_Nodes[left], origin, inverse, hit.t);
-					const float tr = Enter(m_Nodes[right], origin, inverse, hit.t);
-					if (tl <= tr) {
-						if (tr != INF && top < 128) stack[top++] = right;
-						if (tl != INF) { index = left; continue; }
-					} else {
-						if (tl != INF && top < 128) stack[top++] = left;
-						index = right;
-						continue;
-					}
-				}
-				// Next from the stack, skipping nodes now farther than the nearest hit
-				bool found = false;
-				while (top > 0) {
-					index = stack[--top];
-					if (Enter(m_Nodes[index], origin, inverse, hit.t) != INF) {
-						found = true;
-						break;
-					}
-				}
-				if (!found) break;
-			}
-			return hit;
-		}
-
-	private:
-		struct Node {
-			glm::vec3 min{ INF };
-			uint32_t first{}; // leaf: first triangle; inner: the left child (the right one follows it)
-			glm::vec3 max{ -INF };
-			uint32_t count{}; // triangles, 0 for inner nodes
-		};
-
-		// Where the ray enters the node's box, INF when it misses it before `maxT`
-		static float Enter(const Node& node, const glm::vec3& origin, const glm::vec3& inverse, float maxT) {
-			const auto t0 = (node.min - origin) * inverse;
-			const auto t1 = (node.max - origin) * inverse;
-			const auto near = glm::min(t0, t1), far = glm::max(t0, t1);
-			const float enter = std::max(std::max(near.x, near.y), std::max(near.z, 0.0f));
-			const float exit = std::min(std::min(far.x, far.y), std::min(far.z, maxT));
-			return enter <= exit ? enter : INF;
-		}
-
-		struct Triangle {
-			glm::vec3 a, e1, e2;
-			uint32_t index;
-		};
-
-		static glm::vec3 Inverse(const glm::vec3& d) {
-			const auto safe = [](float x) { return 1.0f / (std::abs(x) > 1e-20f ? x : std::copysign(1e-20f, x)); };
-			return { safe(d.x), safe(d.y), safe(d.z) };
-		}
-
-		// Möller-Trumbore; keeps the hit when it's nearer than hit.t
-		static bool Intersect(const Triangle& tri, const glm::vec3& origin, const glm::vec3& direction, uint32_t skip, Hit& hit) {
-			if (tri.index == skip) return false;
-			const auto p = glm::cross(direction, tri.e2);
-			const float det = glm::dot(tri.e1, p);
-			if (det == 0.0f) return false;
-			const float inv = 1.0f / det;
-			const auto s = origin - tri.a;
-			const float u = glm::dot(s, p) * inv;
-			if (u < 0.0f || u > 1.0f) return false;
-			const auto q = glm::cross(s, tri.e1);
-			const float v = glm::dot(direction, q) * inv;
-			if (v < 0.0f || u + v > 1.0f) return false;
-			const float t = glm::dot(tri.e2, q) * inv;
-			if (!(t > 0.0f) || t >= hit.t) return false;
-			hit.t = t;
-			hit.triangle = tri.index;
-			hit.u = u;
-			hit.v = v;
-			return true;
-		}
-
-		void Build(std::vector<uint32_t>& order, const std::vector<glm::vec3>& lo, const std::vector<glm::vec3>& hi, const std::vector<glm::vec3>& centre) {
-			constexpr int BINS = 16;
-			constexpr uint32_t LEAF = 4;
-			struct Task { uint32_t node, first, count; };
-			const auto area = [](const glm::vec3& min, const glm::vec3& max) {
-				const auto d = glm::max(max - min, glm::vec3(0.0f));
-				return d.x * d.y + d.y * d.z + d.z * d.x;
-			};
-			m_Nodes.reserve(order.size() * 2 / LEAF + 1);
-			m_Nodes.push_back({});
-			std::vector<Task> tasks{ { 0, 0, static_cast<uint32_t>(order.size()) } };
-			while (!tasks.empty()) {
-				const auto task = tasks.back();
-				tasks.pop_back();
-				Node node;
-				glm::vec3 cmin(INF), cmax(-INF);
-				for (uint32_t i = task.first; i < task.first + task.count; i++) {
-					node.min = glm::min(node.min, lo[order[i]]);
-					node.max = glm::max(node.max, hi[order[i]]);
-					cmin = glm::min(cmin, centre[order[i]]);
-					cmax = glm::max(cmax, centre[order[i]]);
-				}
-				node.first = task.first;
-				node.count = task.count;
-				int bestAxis = -1;
-				int bestSplit = 0;
-				float bestCost = static_cast<float>(task.count) * area(node.min, node.max); // not splitting
-				if (task.count > LEAF) {
-					for (int axis = 0; axis < 3; axis++) {
-						const float extent = cmax[axis] - cmin[axis];
-						if (!(extent > 0.0f)) continue;
-						struct Bin { glm::vec3 min{ INF }, max{ -INF }; uint32_t count{}; };
-						std::array<Bin, BINS> bins{};
-						const float scale = BINS / extent;
-						for (uint32_t i = task.first; i < task.first + task.count; i++) {
-							const auto t = order[i];
-							const int b = std::min(BINS - 1, static_cast<int>((centre[t][axis] - cmin[axis]) * scale));
-							bins[b].min = glm::min(bins[b].min, lo[t]);
-							bins[b].max = glm::max(bins[b].max, hi[t]);
-							bins[b].count++;
-						}
-						std::array<float, BINS - 1> leftCost{};
-						glm::vec3 lmin(INF), lmax(-INF);
-						uint32_t lcount = 0;
-						for (int b = 0; b < BINS - 1; b++) {
-							lmin = glm::min(lmin, bins[b].min);
-							lmax = glm::max(lmax, bins[b].max);
-							lcount += bins[b].count;
-							leftCost[b] = lcount ? lcount * area(lmin, lmax) : 0.0f;
-						}
-						glm::vec3 rmin(INF), rmax(-INF);
-						uint32_t rcount = 0;
-						for (int b = BINS - 1; b > 0; b--) {
-							rmin = glm::min(rmin, bins[b].min);
-							rmax = glm::max(rmax, bins[b].max);
-							rcount += bins[b].count;
-							const float cost = leftCost[b - 1] + (rcount ? rcount * area(rmin, rmax) : 0.0f);
-							if (rcount > 0 && rcount < task.count && cost < bestCost) {
-								bestCost = cost;
-								bestAxis = axis;
-								bestSplit = b;
-							}
-						}
-					}
-				}
-				if (bestAxis < 0) {
-					m_Nodes[task.node] = node;
-					continue;
-				}
-				const float extent = cmax[bestAxis] - cmin[bestAxis];
-				const float scale = BINS / extent;
-				auto* begin = order.data() + task.first;
-				auto* middle = std::partition(begin, begin + task.count, [&](uint32_t t) {
-					return std::min(BINS - 1, static_cast<int>((centre[t][bestAxis] - cmin[bestAxis]) * scale)) < bestSplit;
-				});
-				const auto leftCount = static_cast<uint32_t>(middle - begin);
-				node.first = static_cast<uint32_t>(m_Nodes.size());
-				node.count = 0;
-				m_Nodes[task.node] = node;
-				m_Nodes.push_back({});
-				m_Nodes.push_back({});
-				tasks.push_back({ node.first, task.first, leftCount });
-				tasks.push_back({ node.first + 1, task.first + leftCount, task.count - leftCount });
-			}
-		}
-
-		std::vector<Node> m_Nodes;
-		std::vector<Triangle> m_Triangles;
-	};
-
 	// LU Toolbox's ground plane: a box 1000 x 1000 x 100 whose top is at LDD y 0, black (a path hitting it ends)
 	float GroundHit(const glm::vec3& origin, const glm::vec3& direction, float maxT) {
 		static const glm::vec3 MIN(-500.0f, -100.0f, -500.0f), MAX(500.0f, 0.0f, 500.0f);
@@ -532,7 +318,7 @@ namespace {
 
 	class Tracer {
 	public:
-		Tracer(const UgcModel::Mesh& mesh, const UgcHsr::Options& options) : m_Mesh(mesh), m_Bvh(mesh), m_Options(options) {
+		Tracer(const UgcModel::Mesh& mesh, const UgcHsr::Options& options) : m_Mesh(mesh), m_Rays(UgcRays::Make(options.rays, mesh)), m_Options(options) {
 			m_Smooth = mesh.normals.size() == mesh.positions.size();
 		}
 
@@ -588,9 +374,9 @@ namespace {
 				const auto& direction = sample.direction;
 				throughput *= sample.throughput;
 				minRayPdf = std::min(minRayPdf, sample.pdf);
-				const auto hit = m_Bvh.Closest(p, direction, self);
+				const auto hit = m_Rays->Closest(p, direction, self);
 				if (m_Options.groundPlane && GroundHit(p, direction, hit.t) < hit.t) return false;
-				if (hit.triangle == UINT32_MAX) return true;
+				if (hit.triangle == UgcRays::NONE) return true;
 				// Past the bounce limits the next surface doesn't scatter (Cycles: Max Bounces, Glossy 4)
 				if (bounce + 1 > m_Options.bounces) return false;
 				if (sample.glossy && ++glossy > GLOSSY_BOUNCES) return false;
@@ -614,7 +400,7 @@ namespace {
 
 	private:
 		const UgcModel::Mesh& m_Mesh;
-		Bvh m_Bvh;
+		std::unique_ptr<UgcRays::Scene> m_Rays; // the nearest hit (never the triangle a ray leaves, as in Cycles)
 		const UgcHsr::Options& m_Options;
 		bool m_Smooth{};
 	};
diff --git a/dUgcServer/Render/UgcHsr.h b/dUgcServer/Render/UgcHsr.h
index d7b525012..a29d9e448 100644
--- a/dUgcServer/Render/UgcHsr.h
+++ b/dUgcServer/Render/UgcHsr.h
@@ -7,6 +7,7 @@
 #include <glm/glm.hpp>
 
 #include "UgcModel.h"
+#include "UgcRays.h"
 
 /**
  * Hidden surface removal as LU Toolbox's Remove Hidden Faces decides it, without its texture: paths are traced from
@@ -23,6 +24,7 @@ namespace UgcHsr {
 		float spacing{ 0.1143f };       // hsr_sample_spacing: LDD units between points (a stud is 0.8: 7 points along it)
 		int minPoints{ 28 };            // hsr_min_points: points on a triangle at least (LU Toolbox bakes 28 texels a triangle)
 		uint64_t seed{};                // of the paths' random numbers (the same seed gives the same result)
+		UgcRays::eBackend rays{};       // what traces the paths' rays (ugc_ray_backend)
 	};
 
 	struct Result {
diff --git a/dUgcServer/Render/UgcRays.cpp b/dUgcServer/Render/UgcRays.cpp
new file mode 100644
index 000000000..f254cfce8
--- /dev/null
+++ b/dUgcServer/Render/UgcRays.cpp
@@ -0,0 +1,560 @@
+#include "UgcRays.h"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstring>
+#include <limits>
+#include <numeric>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+#include <embree4/rtcore.h>
+
+namespace {
+	using UgcRays::Hit;
+	using UgcRays::INF;
+
+	/**
+	 * A bounding volume hierarchy (binned surface area heuristic) over a mesh's triangles, for the paths' rays: the
+	 * nearest hit. A ray never hits the triangle it leaves (`skip`), as in Cycles.
+	 */
+	class ClosestBvh {
+	public:
+		explicit ClosestBvh(const UgcModel::Mesh& mesh) {
+			const size_t count = mesh.TriangleCount();
+			std::vector<uint32_t> order(count);
+			std::iota(order.begin(), order.end(), 0u);
+			std::vector<glm::vec3> lo(count), hi(count), centre(count);
+			for (size_t t = 0; t < count; t++) {
+				const auto& a = mesh.positions[mesh.indices[t * 3]];
+				const auto& b = mesh.positions[mesh.indices[t * 3 + 1]];
+				const auto& c = mesh.positions[mesh.indices[t * 3 + 2]];
+				lo[t] = glm::min(a, glm::min(b, c));
+				hi[t] = glm::max(a, glm::max(b, c));
+				centre[t] = (lo[t] + hi[t]) * 0.5f;
+			}
+			if (count > 0) Build(order, lo, hi, centre);
+			m_Triangles.reserve(count);
+			for (const auto t : order) {
+				const auto& a = mesh.positions[mesh.indices[t * 3]];
+				m_Triangles.push_back({ a, mesh.positions[mesh.indices[t * 3 + 1]] - a, mesh.positions[mesh.indices[t * 3 + 2]] - a, t });
+			}
+		}
+
+		// The nearest triangle along the ray (unit direction) before `maxT`
+		Hit Closest(const glm::vec3& origin, const glm::vec3& direction, uint32_t skip, float maxT = INF) const {
+			Hit hit;
+			hit.t = maxT;
+			if (m_Nodes.empty()) return hit;
+			const auto inverse = Inverse(direction);
+			uint32_t stack[128];
+			int top = 0;
+			uint32_t index = 0;
+			while (true) {
+				const auto& node = m_Nodes[index];
+				if (node.count > 0) {
+					for (uint32_t i = node.first; i < node.first + node.count; i++) Intersect(m_Triangles[i], origin, direction, skip, hit);
+				} else {
+					const uint32_t left = node.first, right = node.first + 1;
+					const float tl = Enter(m_Nodes[left], origin, inverse, hit.t);
+					const float tr = Enter(m_Nodes[right], origin, inverse, hit.t);
+					if (tl <= tr) {
+						if (tr != INF && top < 128) stack[top++] = right;
+						if (tl != INF) { index = left; continue; }
+					} else {
+						if (tl != INF && top < 128) stack[top++] = left;
+						index = right;
+						continue;
+					}
+				}
+				// Next from the stack, skipping nodes now farther than the nearest hit
+				bool found = false;
+				while (top > 0) {
+					index = stack[--top];
+					if (Enter(m_Nodes[index], origin, inverse, hit.t) != INF) {
+						found = true;
+						break;
+					}
+				}
+				if (!found) break;
+			}
+			return hit;
+		}
+
+	private:
+		struct Node {
+			glm::vec3 min{ INF };
+			uint32_t first{}; // leaf: first triangle; inner: the left child (the right one follows it)
+			glm::vec3 max{ -INF };
+			uint32_t count{}; // triangles, 0 for inner nodes
+		};
+
+		// Where the ray enters the node's box, INF when it misses it before `maxT`
+		static float Enter(const Node& node, const glm::vec3& origin, const glm::vec3& inverse, float maxT) {
+			const auto t0 = (node.min - origin) * inverse;
+			const auto t1 = (node.max - origin) * inverse;
+			const auto near = glm::min(t0, t1), far = glm::max(t0, t1);
+			const float enter = std::max(std::max(near.x, near.y), std::max(near.z, 0.0f));
+			const float exit = std::min(std::min(far.x, far.y), std::min(far.z, maxT));
+			return enter <= exit ? enter : INF;
+		}
+
+		struct Triangle {
+			glm::vec3 a, e1, e2;
+			uint32_t index;
+		};
+
+		static glm::vec3 Inverse(const glm::vec3& d) {
+			const auto safe = [](float x) { return 1.0f / (std::abs(x) > 1e-20f ? x : std::copysign(1e-20f, x)); };
+			return { safe(d.x), safe(d.y), safe(d.z) };
+		}
+
+		// Möller-Trumbore; keeps the hit when it's nearer than hit.t
+		static bool Intersect(const Triangle& tri, const glm::vec3& origin, const glm::vec3& direction, uint32_t skip, Hit& hit) {
+			if (tri.index == skip) return false;
+			const auto p = glm::cross(direction, tri.e2);
+			const float det = glm::dot(tri.e1, p);
+			if (det == 0.0f) return false;
+			const float inv = 1.0f / det;
+			const auto s = origin - tri.a;
+			const float u = glm::dot(s, p) * inv;
+			if (u < 0.0f || u > 1.0f) return false;
+			const auto q = glm::cross(s, tri.e1);
+			const float v = glm::dot(direction, q) * inv;
+			if (v < 0.0f || u + v > 1.0f) return false;
+			const float t = glm::dot(tri.e2, q) * inv;
+			if (!(t > 0.0f) || t >= hit.t) return false;
+			hit.t = t;
+			hit.triangle = tri.index;
+			hit.u = u;
+			hit.v = v;
+			return true;
+		}
+
+		void Build(std::vector<uint32_t>& order, const std::vector<glm::vec3>& lo, const std::vector<glm::vec3>& hi, const std::vector<glm::vec3>& centre) {
+			constexpr int BINS = 16;
+			constexpr uint32_t LEAF = 4;
+			struct Task { uint32_t node, first, count; };
+			const auto area = [](const glm::vec3& min, const glm::vec3& max) {
+				const auto d = glm::max(max - min, glm::vec3(0.0f));
+				return d.x * d.y + d.y * d.z + d.z * d.x;
+			};
+			m_Nodes.reserve(order.size() * 2 / LEAF + 1);
+			m_Nodes.push_back({});
+			std::vector<Task> tasks{ { 0, 0, static_cast<uint32_t>(order.size()) } };
+			while (!tasks.empty()) {
+				const auto task = tasks.back();
+				tasks.pop_back();
+				Node node;
+				glm::vec3 cmin(INF), cmax(-INF);
+				for (uint32_t i = task.first; i < task.first + task.count; i++) {
+					node.min = glm::min(node.min, lo[order[i]]);
+					node.max = glm::max(node.max, hi[order[i]]);
+					cmin = glm::min(cmin, centre[order[i]]);
+					cmax = glm::max(cmax, centre[order[i]]);
+				}
+				node.first = task.first;
+				node.count = task.count;
+				int bestAxis = -1;
+				int bestSplit = 0;
+				float bestCost = static_cast<float>(task.count) * area(node.min, node.max); // not splitting
+				if (task.count > LEAF) {
+					for (int axis = 0; axis < 3; axis++) {
+						const float extent = cmax[axis] - cmin[axis];
+						if (!(extent > 0.0f)) continue;
+						struct Bin { glm::vec3 min{ INF }, max{ -INF }; uint32_t count{}; };
+						std::array<Bin, BINS> bins{};
+						const float scale = BINS / extent;
+						for (uint32_t i = task.first; i < task.first + task.count; i++) {
+							const auto t = order[i];
+							const int b = std::min(BINS - 1, static_cast<int>((centre[t][axis] - cmin[axis]) * scale));
+							bins[b].min = glm::min(bins[b].min, lo[t]);
+							bins[b].max = glm::max(bins[b].max, hi[t]);
+							bins[b].count++;
+						}
+						std::array<float, BINS - 1> leftCost{};
+						glm::vec3 lmin(INF), lmax(-INF);
+						uint32_t lcount = 0;
+						for (int b = 0; b < BINS - 1; b++) {
+							lmin = glm::min(lmin, bins[b].min);
+							lmax = glm::max(lmax, bins[b].max);
+							lcount += bins[b].count;
+							leftCost[b] = lcount ? lcount * area(lmin, lmax) : 0.0f;
+						}
+						glm::vec3 rmin(INF), rmax(-INF);
+						uint32_t rcount = 0;
+						for (int b = BINS - 1; b > 0; b--) {
+							rmin = glm::min(rmin, bins[b].min);
+							rmax = glm::max(rmax, bins[b].max);
+							rcount += bins[b].count;
+							const float cost = leftCost[b - 1] + (rcount ? rcount * area(rmin, rmax) : 0.0f);
+							if (rcount > 0 && rcount < task.count && cost < bestCost) {
+								bestCost = cost;
+								bestAxis = axis;
+								bestSplit = b;
+							}
+						}
+					}
+				}
+				if (bestAxis < 0) {
+					m_Nodes[task.node] = node;
+					continue;
+				}
+				const float extent = cmax[bestAxis] - cmin[bestAxis];
+				const float scale = BINS / extent;
+				auto* begin = order.data() + task.first;
+				auto* middle = std::partition(begin, begin + task.count, [&](uint32_t t) {
+					return std::min(BINS - 1, static_cast<int>((centre[t][bestAxis] - cmin[bestAxis]) * scale)) < bestSplit;
+				});
+				const auto leftCount = static_cast<uint32_t>(middle - begin);
+				node.first = static_cast<uint32_t>(m_Nodes.size());
+				node.count = 0;
+				m_Nodes[task.node] = node;
+				m_Nodes.push_back({});
+				m_Nodes.push_back({});
+				tasks.push_back({ node.first, task.first, leftCount });
+				tasks.push_back({ node.first + 1, task.first + leftCount, task.count - leftCount });
+			}
+		}
+
+		std::vector<Node> m_Nodes;
+		std::vector<Triangle> m_Triangles;
+	};
+
+	// A bounding volume hierarchy over a mesh's triangles (median splits), for the occlusion rays: any hit. The mesh
+	// is only read while it is built.
+	class AnyBvh {
+	public:
+		explicit AnyBvh(const UgcModel::Mesh& mesh) : m_Mesh(mesh) {
+			const size_t count = mesh.TriangleCount();
+			m_Order.resize(count);
+			std::iota(m_Order.begin(), m_Order.end(), 0u);
+			m_Centers.resize(count);
+			for (size_t t = 0; t < count; t++) m_Centers[t] = (Vertex(t, 0) + Vertex(t, 1) + Vertex(t, 2)) / 3.0f;
+			if (count > 0) Build(0, static_cast<uint32_t>(count));
+			Flatten();
+		}
+
+		// Whether a ray from `origin` along `direction` (unit) hits a triangle further than `minDistance` and nearer
+		// than `maxDistance`
+		bool Hits(const glm::vec3& origin, const glm::vec3& direction, float minDistance, float maxDistance) const {
+			if (m_Nodes.empty()) return false;
+			const glm::vec3 inverse(1.0f / (std::abs(direction.x) > 1e-12f ? direction.x : 1e-12f), 1.0f / (std::abs(direction.y) > 1e-12f ? direction.y : 1e-12f),
+				1.0f / (std::abs(direction.z) > 1e-12f ? direction.z : 1e-12f));
+			uint32_t stack[64];
+			int top = 0;
+			stack[top++] = 0;
+			while (top > 0) {
+				const auto& node = m_Nodes[stack[--top]];
+				if (!BoxHit(node, origin, inverse, maxDistance)) continue;
+				if (node.count > 0) {
+					for (uint32_t i = node.first; i < node.first + node.count; i++) {
+						if (TriangleHit(m_Triangles[i], origin, direction, minDistance, maxDistance)) return true;
+					}
+				} else if (top < 62) {
+					stack[top++] = node.first;
+					stack[top++] = node.first + 1;
+				}
+			}
+			return false;
+		}
+
+	private:
+		struct Node {
+			glm::vec3 min{};
+			glm::vec3 max{};
+			uint32_t first{}; // leaf: first triangle in m_Order; inner: the first of two children
+			uint32_t count{}; // triangles, 0 for inner nodes
+		};
+
+		glm::vec3 Vertex(size_t t, int k) const { return m_Mesh.positions[m_Mesh.indices[t * 3 + k]]; }
+
+		void Build(uint32_t first, uint32_t count) {
+			// Iterative, so deep trees don't use the stack
+			struct Task { uint32_t node, first, count; };
+			m_Nodes.push_back({});
+			std::vector<Task> tasks{ { 0, first, count } };
+			while (!tasks.empty()) {
+				const auto task = tasks.back();
+				tasks.pop_back();
+				Node node;
+				node.min = glm::vec3(INF);
+				node.max = glm::vec3(-INF);
+				glm::vec3 centerMin(INF), centerMax(-INF);
+				for (uint32_t i = task.first; i < task.first + task.count; i++) {
+					for (int k = 0; k < 3; k++) {
+						node.min = glm::min(node.min, Vertex(m_Order[i], k));
+						node.max = glm::max(node.max, Vertex(m_Order[i], k));
+					}
+					centerMin = glm::min(centerMin, m_Centers[m_Order[i]]);
+					centerMax = glm::max(centerMax, m_Centers[m_Order[i]]);
+				}
+				const auto extent = centerMax - centerMin;
+				const int axis = extent.x >= extent.y && extent.x >= extent.z ? 0 : extent.y >= extent.z ? 1 : 2;
+				if (task.count <= 4 || extent[axis] <= 0.0f) {
+					node.first = task.first;
+					node.count = task.count;
+					m_Nodes[task.node] = node;
+					continue;
+				}
+				const uint32_t half = task.count / 2;
+				auto* begin = m_Order.data() + task.first;
+				std::nth_element(begin, begin + half, begin + task.count, [&](uint32_t a, uint32_t b) { return m_Centers[a][axis] < m_Centers[b][axis]; });
+				node.first = static_cast<uint32_t>(m_Nodes.size());
+				node.count = 0;
+				m_Nodes[task.node] = node;
+				m_Nodes.push_back({});
+				m_Nodes.push_back({});
+				tasks.push_back({ node.first, task.first, half });
+				tasks.push_back({ node.first + 1, task.first + half, task.count - half });
+			}
+		}
+
+		static bool BoxHit(const Node& node, const glm::vec3& origin, const glm::vec3& inverse, float maxDistance) {
+			const auto t0 = (node.min - origin) * inverse;
+			const auto t1 = (node.max - origin) * inverse;
+			const auto near = glm::min(t0, t1), far = glm::max(t0, t1);
+			const float enter = std::max(std::max(near.x, near.y), std::max(near.z, 0.0f));
+			const float exit = std::min(std::min(far.x, far.y), std::min(far.z, maxDistance));
+			return enter <= exit;
+		}
+
+		struct Triangle {
+			glm::vec3 a, e1, e2;
+		};
+
+		// The triangles in leaf order, edges worked out once (the rays read them far more often than the tree is built)
+		void Flatten() {
+			m_Triangles.reserve(m_Order.size());
+			for (const auto t : m_Order) {
+				const auto a = Vertex(t, 0);
+				m_Triangles.push_back({ a, Vertex(t, 1) - a, Vertex(t, 2) - a });
+			}
+		}
+
+		static bool TriangleHit(const Triangle& triangle, const glm::vec3& origin, const glm::vec3& direction, float minDistance, float maxDistance) {
+			const auto& a = triangle.a;
+			const auto& e1 = triangle.e1;
+			const auto& e2 = triangle.e2;
+			const auto p = glm::cross(direction, e2);
+			const float det = glm::dot(e1, p);
+			if (std::abs(det) < 1e-12f) return false;
+			const float inv = 1.0f / det;
+			const auto s = origin - a;
+			const float u = glm::dot(s, p) * inv;
+			if (u < 0.0f || u > 1.0f) return false;
+			const auto q = glm::cross(s, e1);
+			const float v = glm::dot(direction, q) * inv;
+			if (v < 0.0f || u + v > 1.0f) return false;
+			const float distance = glm::dot(e2, q) * inv;
+			return distance > minDistance && distance < maxDistance;
+		}
+
+		const UgcModel::Mesh& m_Mesh;
+		std::vector<Triangle> m_Triangles;
+		std::vector<uint32_t> m_Order;
+		std::vector<glm::vec3> m_Centers;
+		std::vector<Node> m_Nodes;
+	};
+
+	/**
+	 * builtin: the hierarchy each query had before the backends (the paths' nearest hits: ClosestBvh; the occlusion
+	 * rays: AnyBvh), each built the first time it is asked, so the results are exactly what they were.
+	 */
+	class BuiltinScene final : public UgcRays::Scene {
+	public:
+		explicit BuiltinScene(const UgcModel::Mesh& mesh) {
+			m_Mesh.positions = mesh.positions;
+			m_Mesh.indices = mesh.indices;
+		}
+
+		Hit Closest(const glm::vec3& origin, const glm::vec3& direction, uint32_t skip, float maxT) const override {
+			if (!m_Closest) m_Closest = std::make_unique<ClosestBvh>(m_Mesh);
+			return m_Closest->Closest(origin, direction, skip, maxT);
+		}
+
+		bool Occluded(const glm::vec3& origin, const glm::vec3& direction, float minT, float maxT) const override {
+			if (!m_Any) m_Any = std::make_unique<AnyBvh>(m_Mesh);
+			return m_Any->Hits(origin, direction, minT, maxT);
+		}
+
+	private:
+		UgcModel::Mesh m_Mesh; // positions and indices only
+		mutable std::unique_ptr<ClosestBvh> m_Closest;
+		mutable std::unique_ptr<AnyBvh> m_Any;
+	};
+
+	// Embree's device for the thread: one per thread, with no threads of its own (threads=1: the thread that commits
+	// a scene builds it), released when the thread ends
+	struct EmbreeDevice {
+		RTCDevice device{};
+
+		EmbreeDevice() {
+			device = rtcNewDevice("threads=1,set_affinity=0,verbose=0");
+			if (!device) throw std::runtime_error("Embree: no device (error " + std::to_string(rtcGetDeviceError(nullptr)) + ")");
+		}
+		~EmbreeDevice() { rtcReleaseDevice(device); }
+		EmbreeDevice(const EmbreeDevice&) = delete;
+		EmbreeDevice& operator=(const EmbreeDevice&) = delete;
+
+		static RTCDevice Get() {
+			thread_local EmbreeDevice instance;
+			return instance.device;
+		}
+	};
+
+	void EmbreeCheck(RTCDevice device, const char* what) {
+		const auto error = rtcGetDeviceError(device);
+		if (error != RTC_ERROR_NONE) throw std::runtime_error(std::string("Embree: ") + what + " failed (error " + std::to_string(error) + ")");
+	}
+
+	// A ray query's context, with the triangle the ray may not hit
+	struct SkipContext {
+		RTCRayQueryContext base;
+		uint32_t skip;
+	};
+
+	void SkipFilter(const RTCFilterFunctionNArguments* args) {
+		const auto* context = reinterpret_cast<const SkipContext*>(args->context);
+		for (unsigned int i = 0; i < args->N; i++) {
+			if (args->valid[i] != 0 && RTCHitN_primID(args->hit, args->N, i) == context->skip) args->valid[i] = 0;
+		}
+	}
+
+	RTCRay EmbreeRay(const glm::vec3& origin, const glm::vec3& direction, float minT, float maxT) {
+		RTCRay ray{};
+		ray.org_x = origin.x;
+		ray.org_y = origin.y;
+		ray.org_z = origin.z;
+		ray.dir_x = direction.x;
+		ray.dir_y = direction.y;
+		ray.dir_z = direction.z;
+		ray.tnear = minT;
+		ray.tfar = maxT;
+		ray.mask = 0xFFFFFFFFu;
+		ray.flags = 0;
+		ray.time = 0.0f;
+		return ray;
+	}
+
+	/**
+	 * embree: one triangle geometry, built at high quality and traced watertight (a ray through the edge two
+	 * triangles share hits one of them). The triangle a path leaves is skipped by a filter function.
+	 */
+	class EmbreeScene final : public UgcRays::Scene {
+	public:
+		explicit EmbreeScene(const UgcModel::Mesh& mesh) {
+			const auto device = EmbreeDevice::Get();
+			const size_t triangles = mesh.TriangleCount();
+			m_Scene = rtcNewScene(device);
+			EmbreeCheck(device, "rtcNewScene");
+			rtcSetSceneBuildQuality(m_Scene, RTC_BUILD_QUALITY_HIGH);
+			rtcSetSceneFlags(m_Scene, RTC_SCENE_FLAG_ROBUST | RTC_SCENE_FLAG_FILTER_FUNCTION_IN_ARGUMENTS);
+			if (triangles > 0 && !mesh.positions.empty()) {
+				const RTCGeometry geometry = rtcNewGeometry(device, RTC_GEOMETRY_TYPE_TRIANGLE);
+				// Copied into Embree's own buffers, which are padded for its vector loads (a std::vector isn't)
+				auto* vertices = static_cast<float*>(rtcSetNewGeometryBuffer(geometry, RTC_BUFFER_TYPE_VERTEX, 0, RTC_FORMAT_FLOAT3, sizeof(float) * 3, mesh.positions.size()));
+				auto* indices = static_cast<uint32_t*>(rtcSetNewGeometryBuffer(geometry, RTC_BUFFER_TYPE_INDEX, 0, RTC_FORMAT_UINT3, sizeof(uint32_t) * 3, triangles));
+				if (!vertices || !indices) {
+					rtcReleaseGeometry(geometry);
+					EmbreeCheck(device, "rtcSetNewGeometryBuffer");
+					throw std::runtime_error("Embree: no geometry buffers");
+				}
+				for (size_t i = 0; i < mesh.positions.size(); i++) {
+					vertices[i * 3] = mesh.positions[i].x;
+					vertices[i * 3 + 1] = mesh.positions[i].y;
+					vertices[i * 3 + 2] = mesh.positions[i].z;
+				}
+				std::memcpy(indices, mesh.indices.data(), sizeof(uint32_t) * 3 * triangles);
+				rtcSetGeometryEnableFilterFunctionFromArguments(geometry, true);
+				rtcCommitGeometry(geometry);
+				rtcAttachGeometry(m_Scene, geometry);
+				rtcReleaseGeometry(geometry);
+				m_Empty = false;
+			}
+			rtcCommitScene(m_Scene);
+			EmbreeCheck(device, "rtcCommitScene");
+		}
+
+		~EmbreeScene() override {
+			if (m_Scene) rtcReleaseScene(m_Scene);
+		}
+
+		EmbreeScene(const EmbreeScene&) = delete;
+		EmbreeScene& operator=(const EmbreeScene&) = delete;
+
+		Hit Closest(const glm::vec3& origin, const glm::vec3& direction, uint32_t skip, float maxT) const override {
+			Hit hit;
+			hit.t = maxT;
+			if (m_Empty) return hit;
+			RTCRayHit query{};
+			// Further than 0, as the builtin test
+			query.ray = EmbreeRay(origin, direction, std::numeric_limits<float>::min(), maxT);
+			query.hit.geomID = RTC_INVALID_GEOMETRY_ID;
+			query.hit.primID = RTC_INVALID_GEOMETRY_ID;
+			SkipContext context{};
+			rtcInitRayQueryContext(&context.base);
+			context.skip = skip;
+			RTCIntersectArguments arguments;
+			rtcInitIntersectArguments(&arguments);
+			arguments.context = &context.base;
+			if (skip != UgcRays::NONE) {
+				arguments.filter = SkipFilter;
+				arguments.flags = static_cast<RTCRayQueryFlags>(arguments.flags | RTC_RAY_QUERY_FLAG_INVOKE_ARGUMENT_FILTER);
+			}
+			rtcIntersect1(m_Scene, &query, &arguments);
+			if (query.hit.geomID == RTC_INVALID_GEOMETRY_ID) return hit;
+			hit.t = query.ray.tfar;
+			hit.triangle = query.hit.primID;
+			hit.u = query.hit.u;
+			hit.v = query.hit.v;
+			return hit;
+		}
+
+		bool Occluded(const glm::vec3& origin, const glm::vec3& direction, float minT, float maxT) const override {
+			if (m_Empty) return false;
+			// Embree counts hits at exactly tnear and tfar; the builtin test doesn't
+			auto ray = EmbreeRay(origin, direction, std::nextafter(minT, INF), std::nextafter(maxT, 0.0f));
+			rtcOccluded1(m_Scene, &ray, nullptr);
+			return ray.tfar < 0.0f;
+		}
+
+	private:
+		RTCScene m_Scene{};
+		bool m_Empty{ true };
+	};
+}
+
+namespace UgcRays {
+	std::string_view Name(eBackend backend) {
+		switch (backend) {
+		case eBackend::EMBREE: return "embree";
+		case eBackend::HIPRT: return "hiprt";
+		default: return "builtin";
+		}
+	}
+
+	std::optional<eBackend> Parse(std::string_view name) {
+		for (const auto backend : { eBackend::BUILTIN, eBackend::EMBREE, eBackend::HIPRT }) {
+			if (Name(backend) == name) return backend;
+		}
+		return std::nullopt;
+	}
+
+	bool Available(eBackend backend) {
+		return backend == eBackend::BUILTIN || backend == eBackend::EMBREE;
+	}
+
+	eBackend Resolve(eBackend wanted) {
+		return Available(wanted) ? wanted : eBackend::EMBREE;
+	}
+
+	std::unique_ptr<Scene> Make(eBackend backend, const UgcModel::Mesh& mesh) {
+		switch (Resolve(backend)) {
+		case eBackend::EMBREE: return std::make_unique<EmbreeScene>(mesh);
+		default: return std::make_unique<BuiltinScene>(mesh);
+		}
+	}
+}
diff --git a/dUgcServer/Render/UgcRays.h b/dUgcServer/Render/UgcRays.h
new file mode 100644
index 000000000..9e95f7f84
--- /dev/null
+++ b/dUgcServer/Render/UgcRays.h
@@ -0,0 +1,59 @@
+#pragma once
+
+#include <cstdint>
+#include <limits>
+#include <memory>
+#include <optional>
+#include <string_view>
+
+#include <glm/glm.hpp>
+
+#include "UgcModel.h"
+
+/**
+ * Rays against a mesh's triangles, for the hidden faces' paths (the nearest hit) and the ambient occlusion rays
+ * (whether anything is hit), by one of several backends (the ugc_ray_backend setting, or per job):
+ *   builtin: the UGC server's own bounding volume hierarchies (the ones it always had)
+ *   embree:  Intel's Embree 4 on the CPU, on the thread that asks (no threads of its own)
+ * A scene is built and traced on the thread that asks, so its time counts towards that thread's CPU time
+ * (UgcThrottle). Scenes aren't shared between threads. docs/UgcServer.md ("Processing options") has the details.
+ */
+namespace UgcRays {
+	constexpr uint32_t NONE = std::numeric_limits<uint32_t>::max();
+	constexpr float INF = std::numeric_limits<float>::infinity();
+
+	enum class eBackend : uint8_t { BUILTIN = 0, EMBREE, HIPRT };
+
+	// The setting's name of a backend (builtin, embree, hiprt)
+	std::string_view Name(eBackend backend);
+	// A backend by its name (case sensitive); nullopt for anything else
+	std::optional<eBackend> Parse(std::string_view name);
+	// Whether this build and machine can use the backend (builtin and embree always)
+	bool Available(eBackend backend);
+	// The backend that is used when `wanted` is asked for: itself, or embree when it isn't available
+	eBackend Resolve(eBackend wanted);
+
+	struct Hit {
+		float t{ INF };
+		uint32_t triangle{ NONE };
+		float u{}, v{}; // weights of the triangle's second and third vertex
+	};
+
+	class Scene {
+	public:
+		virtual ~Scene() = default;
+
+		// The nearest triangle along the ray (unit direction) further than 0 and before `maxT`, never `skip` (the
+		// triangle the ray leaves, as in Cycles); triangle NONE (and t = maxT) when there is none
+		virtual Hit Closest(const glm::vec3& origin, const glm::vec3& direction, uint32_t skip = NONE, float maxT = INF) const = 0;
+
+		// Whether the ray (unit direction) hits a triangle further than `minT` and nearer than `maxT`
+		virtual bool Occluded(const glm::vec3& origin, const glm::vec3& direction, float minT, float maxT) const = 0;
+	};
+
+	/**
+	 * The mesh's triangles as they are now (copied) in the backend Resolve(backend) picks. Throws std::runtime_error
+	 * when the backend fails.
+	 */
+	std::unique_ptr<Scene> Make(eBackend backend, const UgcModel::Mesh& mesh);
+}
diff --git a/dUgcServer/Render/UgcRender.cpp b/dUgcServer/Render/UgcRender.cpp
index 0aa4cd533..38a06901f 100644
--- a/dUgcServer/Render/UgcRender.cpp
+++ b/dUgcServer/Render/UgcRender.cpp
@@ -10,6 +10,7 @@
 
 #include "UgcIconPose.h"
 #include "UgcPalette.h"
+#include "UgcRays.h"
 #include "UgcThrottle.h"
 
 namespace {
@@ -45,140 +46,6 @@ namespace {
 	float ToLinear(float c) { return UgcPalette::SrgbToLinear(std::clamp(c, 0.0f, 1.0f)); }
 	float ToSrgb(float c) { return UgcPalette::LinearToSrgb(std::clamp(c, 0.0f, 1.0f)); }
 
-	// A bounding volume hierarchy over a mesh's triangles, for the occlusion rays
-	class Bvh {
-	public:
-		explicit Bvh(const UgcModel::Mesh& mesh) : m_Mesh(mesh) {
-			const size_t count = mesh.TriangleCount();
-			m_Order.resize(count);
-			std::iota(m_Order.begin(), m_Order.end(), 0u);
-			m_Centers.resize(count);
-			for (size_t t = 0; t < count; t++) m_Centers[t] = (Vertex(t, 0) + Vertex(t, 1) + Vertex(t, 2)) / 3.0f;
-			if (count > 0) Build(0, static_cast<uint32_t>(count));
-			Flatten();
-		}
-
-		// Whether a ray from `origin` along `direction` (unit) hits a triangle nearer than `maxDistance`
-		bool Hits(const glm::vec3& origin, const glm::vec3& direction, float maxDistance) const {
-			if (m_Nodes.empty()) return false;
-			const glm::vec3 inverse(1.0f / (std::abs(direction.x) > 1e-12f ? direction.x : 1e-12f), 1.0f / (std::abs(direction.y) > 1e-12f ? direction.y : 1e-12f),
-				1.0f / (std::abs(direction.z) > 1e-12f ? direction.z : 1e-12f));
-			uint32_t stack[64];
-			int top = 0;
-			stack[top++] = 0;
-			while (top > 0) {
-				const auto& node = m_Nodes[stack[--top]];
-				if (!BoxHit(node, origin, inverse, maxDistance)) continue;
-				if (node.count > 0) {
-					for (uint32_t i = node.first; i < node.first + node.count; i++) {
-						if (TriangleHit(m_Triangles[i], origin, direction, maxDistance)) return true;
-					}
-				} else if (top < 62) {
-					stack[top++] = node.first;
-					stack[top++] = node.first + 1;
-				}
-			}
-			return false;
-		}
-
-	private:
-		struct Node {
-			glm::vec3 min{};
-			glm::vec3 max{};
-			uint32_t first{}; // leaf: first triangle in m_Order; inner: the first of two children
-			uint32_t count{}; // triangles, 0 for inner nodes
-		};
-
-		glm::vec3 Vertex(size_t t, int k) const { return m_Mesh.positions[m_Mesh.indices[t * 3 + k]]; }
-
-		void Build(uint32_t first, uint32_t count) {
-			// Iterative, so deep trees don't use the stack
-			struct Task { uint32_t node, first, count; };
-			m_Nodes.push_back({});
-			std::vector<Task> tasks{ { 0, first, count } };
-			while (!tasks.empty()) {
-				const auto task = tasks.back();
-				tasks.pop_back();
-				Node node;
-				node.min = glm::vec3(INF);
-				node.max = glm::vec3(-INF);
-				glm::vec3 centerMin(INF), centerMax(-INF);
-				for (uint32_t i = task.first; i < task.first + task.count; i++) {
-					for (int k = 0; k < 3; k++) {
-						node.min = glm::min(node.min, Vertex(m_Order[i], k));
-						node.max = glm::max(node.max, Vertex(m_Order[i], k));
-					}
-					centerMin = glm::min(centerMin, m_Centers[m_Order[i]]);
-					centerMax = glm::max(centerMax, m_Centers[m_Order[i]]);
-				}
-				const auto extent = centerMax - centerMin;
-				const int axis = extent.x >= extent.y && extent.x >= extent.z ? 0 : extent.y >= extent.z ? 1 : 2;
-				if (task.count <= 4 || extent[axis] <= 0.0f) {
-					node.first = task.first;
-					node.count = task.count;
-					m_Nodes[task.node] = node;
-					continue;
-				}
-				const uint32_t half = task.count / 2;
-				auto* begin = m_Order.data() + task.first;
-				std::nth_element(begin, begin + half, begin + task.count, [&](uint32_t a, uint32_t b) { return m_Centers[a][axis] < m_Centers[b][axis]; });
-				node.first = static_cast<uint32_t>(m_Nodes.size());
-				node.count = 0;
-				m_Nodes[task.node] = node;
-				m_Nodes.push_back({});
-				m_Nodes.push_back({});
-				tasks.push_back({ node.first, task.first, half });
-				tasks.push_back({ node.first + 1, task.first + half, task.count - half });
-			}
-		}
-
-		static bool BoxHit(const Node& node, const glm::vec3& origin, const glm::vec3& inverse, float maxDistance) {
-			const auto t0 = (node.min - origin) * inverse;
-			const auto t1 = (node.max - origin) * inverse;
-			const auto near = glm::min(t0, t1), far = glm::max(t0, t1);
-			const float enter = std::max(std::max(near.x, near.y), std::max(near.z, 0.0f));
-			const float exit = std::min(std::min(far.x, far.y), std::min(far.z, maxDistance));
-			return enter <= exit;
-		}
-
-		struct Triangle {
-			glm::vec3 a, e1, e2;
-		};
-
-		// The triangles in leaf order, edges worked out once (the rays read them far more often than the tree is built)
-		void Flatten() {
-			m_Triangles.reserve(m_Order.size());
-			for (const auto t : m_Order) {
-				const auto a = Vertex(t, 0);
-				m_Triangles.push_back({ a, Vertex(t, 1) - a, Vertex(t, 2) - a });
-			}
-		}
-
-		static bool TriangleHit(const Triangle& triangle, const glm::vec3& origin, const glm::vec3& direction, float maxDistance) {
-			const auto& a = triangle.a;
-			const auto& e1 = triangle.e1;
-			const auto& e2 = triangle.e2;
-			const auto p = glm::cross(direction, e2);
-			const float det = glm::dot(e1, p);
-			if (std::abs(det) < 1e-12f) return false;
-			const float inv = 1.0f / det;
-			const auto s = origin - a;
-			const float u = glm::dot(s, p) * inv;
-			if (u < 0.0f || u > 1.0f) return false;
-			const auto q = glm::cross(s, e1);
-			const float v = glm::dot(direction, q) * inv;
-			if (v < 0.0f || u + v > 1.0f) return false;
-			const float distance = glm::dot(e2, q) * inv;
-			return distance > 1e-4f && distance < maxDistance;
-		}
-
-		const UgcModel::Mesh& m_Mesh;
-		std::vector<Triangle> m_Triangles;
-		std::vector<uint32_t> m_Order;
-		std::vector<glm::vec3> m_Centers;
-		std::vector<Node> m_Nodes;
-	};
-
 	float RadicalInverse(uint32_t bits) {
 		bits = (bits << 16u) | (bits >> 16u);
 		bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
@@ -216,10 +83,10 @@ namespace {
 }
 
 namespace UgcRender {
-	std::vector<float> AmbientOcclusion(const UgcModel::Mesh& mesh, const UgcModel::Mesh& occluders, float distance, int samples) {
+	std::vector<float> AmbientOcclusion(const UgcModel::Mesh& mesh, const UgcModel::Mesh& occluders, float distance, int samples, UgcRays::eBackend rays) {
 		std::vector<float> ao(mesh.positions.size(), 1.0f);
 		if (occluders.Empty() || samples <= 0 || distance <= 0.0f || mesh.normals.size() != mesh.positions.size()) return ao;
-		const Bvh bvh(occluders);
+		const auto scene = UgcRays::Make(rays, occluders);
 		const auto count = static_cast<uint32_t>(samples);
 		// Vertices at the same place facing the same way (bricks' shared corners) are worked out once
 		struct Key {
@@ -259,7 +126,7 @@ namespace UgcRender {
 				const float phi = 2.0f * 3.14159265f * std::fmod(RadicalInverse(i) + turn, 1.0f);
 				const float r = std::sqrt(u), z = std::sqrt(std::max(0.0f, 1.0f - u));
 				const auto direction = tangent * (r * std::cos(phi)) + bitangent * (r * std::sin(phi)) + normal * z;
-				if (!bvh.Hits(origin, direction, distance)) open++;
+				if (!scene->Occluded(origin, direction, 1e-4f, distance)) open++;
 			}
 			ao[v] = static_cast<float>(open) / static_cast<float>(count);
 			known.emplace(key, ao[v]);
@@ -270,7 +137,7 @@ namespace UgcRender {
 	std::vector<float> BakeAo(UgcModel::Model& model, const AoOptions& options) {
 		auto& opaque = model.opaque;
 		if (!options.enabled || opaque.Empty()) return {};
-		auto ao = AmbientOcclusion(opaque, opaque, options.distance, options.samples);
+		auto ao = AmbientOcclusion(opaque, opaque, options.distance, options.samples, options.rays);
 		const float strength = std::clamp(options.strength, 0.0f, 1.0f);
 		for (size_t v = 0; v < opaque.colors.size() && v < ao.size(); v++) {
 			glm::vec3 lit(1.0f - strength * (1.0f - ao[v]));
@@ -318,7 +185,7 @@ namespace UgcRender {
 		// Ambient occlusion darkens the world light (opaque bricks only, as they are what occludes)
 		std::vector<float> ao;
 		if (options.ao.enabled) {
-			ao = opaqueAo && opaqueAo->size() == model.opaque.positions.size() ? *opaqueAo : AmbientOcclusion(model.opaque, model.opaque, options.ao.distance, options.ao.samples);
+			ao = opaqueAo && opaqueAo->size() == model.opaque.positions.size() ? *opaqueAo : AmbientOcclusion(model.opaque, model.opaque, options.ao.distance, options.ao.samples, options.ao.rays);
 		}
 
 		// The sun's shadows: a depth map seen from the sun, looked up with a few taps for the sun's soft edge
diff --git a/dUgcServer/Render/UgcRender.h b/dUgcServer/Render/UgcRender.h
index 6dd6c4535..6e6e66f73 100644
--- a/dUgcServer/Render/UgcRender.h
+++ b/dUgcServer/Render/UgcRender.h
@@ -7,6 +7,7 @@
 
 #include "UgcGlitter.h"
 #include "UgcModel.h"
+#include "UgcRays.h"
 
 /**
  * A small software rasterizer (no GPU or display needed) for what the UGC server draws: the icons, and the occlusion
@@ -26,6 +27,7 @@ namespace UgcRender {
 		int samples{ 64 };          // rays per vertex (AO Samples)
 		float strength{ 1.0f };     // 0 leaves the colors, 1 is the full bake
 		float glowStrength{ 6.0f }; // what glowing colors add to the light (Glow Strength 3 x Glow Multiplier 2)
+		UgcRays::eBackend rays{};   // what traces the occlusion rays (ugc_ray_backend)
 	};
 
 	/**
@@ -74,7 +76,8 @@ namespace UgcRender {
 	 * normal, the same pattern every time) that leave without hitting a triangle of `occluders` within `distance`.
 	 * 1 is open, 0 fully hidden.
 	 */
-	std::vector<float> AmbientOcclusion(const UgcModel::Mesh& mesh, const UgcModel::Mesh& occluders, float distance, int samples);
+	std::vector<float> AmbientOcclusion(const UgcModel::Mesh& mesh, const UgcModel::Mesh& occluders, float distance, int samples,
+		UgcRays::eBackend rays = UgcRays::eBackend::BUILTIN);
 
 	/**
 	 * LU Toolbox's Bake Lighting with AO Only (its defaults): the opaque mesh's occlusion (transparent bricks are hidden
diff --git a/tests/dUgcTests/UgcTests.cpp b/tests/dUgcTests/UgcTests.cpp
index e08277a73..44a5f168f 100644
--- a/tests/dUgcTests/UgcTests.cpp
+++ b/tests/dUgcTests/UgcTests.cpp
@@ -22,6 +22,7 @@
 #include "UgcKeys.h"
 #include "UgcModular.h"
 #include "UgcPalette.h"
+#include "UgcRays.h"
 #include "UgcRender.h"
 #include "UgcStorage.h"
 #include "UgcThrottle.h"
@@ -1980,3 +1981,126 @@ TEST(UgcJobs, IconDrawnAgainChangesTheMakesTime) {
 	EXPECT_EQ(nlohmann::json::parse(*old)["ms"]["total"], 520);
 	EXPECT_FALSE(UgcJobs::WithIconTime("not json", 20.0, change));
 }
+
+namespace {
+	// A cluttered scene for the ray backends: the room with its doorway and box, a stack of boxes that touch and
+	// overlap, and a staircase of thin slabs
+	UgcModel::Mesh Clutter() {
+		auto mesh = Room(false);
+		for (int i = 0; i < 4; i++) AddBox(mesh, glm::vec3(-1.5f + i * 0.4f, -2.0f + i * 0.3f, -1.5f), glm::vec3(-1.0f + i * 0.4f, -1.6f + i * 0.3f, -0.8f));
+		for (int i = 0; i < 6; i++) AddBox(mesh, glm::vec3(0.5f, -2.0f + i * 0.2f, -1.8f + i * 0.25f), glm::vec3(1.8f, -1.95f + i * 0.2f, -1.4f + i * 0.25f));
+		return mesh;
+	}
+
+	// The backends other than builtin that this build and machine can use
+	std::vector<UgcRays::eBackend> OtherBackends() {
+		std::vector<UgcRays::eBackend> backends;
+		for (const auto backend : { UgcRays::eBackend::EMBREE, UgcRays::eBackend::HIPRT }) {
+			if (UgcRays::Available(backend)) backends.push_back(backend);
+		}
+		return backends;
+	}
+}
+
+TEST(UgcRays, NamesAndFallback) {
+	for (const auto backend : { UgcRays::eBackend::BUILTIN, UgcRays::eBackend::EMBREE, UgcRays::eBackend::HIPRT }) {
+		EXPECT_EQ(UgcRays::Parse(UgcRays::Name(backend)), backend);
+	}
+	EXPECT_FALSE(UgcRays::Parse("optix"));
+	EXPECT_TRUE(UgcRays::Available(UgcRays::eBackend::BUILTIN));
+	EXPECT_TRUE(UgcRays::Available(UgcRays::eBackend::EMBREE));
+	// A backend this machine can't use falls back to embree
+	EXPECT_EQ(UgcRays::Resolve(UgcRays::eBackend::HIPRT), UgcRays::Available(UgcRays::eBackend::HIPRT) ? UgcRays::eBackend::HIPRT : UgcRays::eBackend::EMBREE);
+	EXPECT_EQ(UgcRays::Resolve(UgcRays::eBackend::BUILTIN), UgcRays::eBackend::BUILTIN);
+	// An empty mesh is hit by nothing
+	for (const auto backend : { UgcRays::eBackend::BUILTIN, UgcRays::eBackend::EMBREE }) {
+		const auto scene = UgcRays::Make(backend, UgcModel::Mesh{});
+		EXPECT_EQ(scene->Closest(glm::vec3(0.0f), glm::vec3(0, 1, 0)).triangle, UgcRays::NONE);
+		EXPECT_FALSE(scene->Occluded(glm::vec3(0.0f), glm::vec3(0, 1, 0), 0.0f, 10.0f));
+	}
+}
+
+TEST(UgcRays, BackendsFindTheSameHits) {
+	// Rays in every direction from points around the clutter: the other backends find the same nearest triangle at
+	// the same distance, and agree on what blocks. A ray through an edge two triangles share may hit either, at the
+	// same distance.
+	const auto mesh = Clutter();
+	const auto builtin = UgcRays::Make(UgcRays::eBackend::BUILTIN, mesh);
+	for (const auto backend : OtherBackends()) {
+		const auto other = UgcRays::Make(backend, mesh);
+		uint64_t state = 12345;
+		const auto next = [&state]() {
+			state = state * 6364136223846793005ull + 1442695040888963407ull;
+			return static_cast<float>(state >> 40) / 16777216.0f;
+		};
+		size_t hits = 0, sameTriangle = 0;
+		for (int i = 0; i < 20000; i++) {
+			const glm::vec3 origin(next() * 3.8f - 1.9f, next() * 3.8f - 1.9f, next() * 3.8f - 1.9f);
+			const auto direction = glm::normalize(glm::vec3(next() - 0.5f, next() - 0.5f, next() - 0.5f) + glm::vec3(1e-4f));
+			const auto skip = static_cast<uint32_t>(next() * static_cast<float>(mesh.TriangleCount()));
+			const auto a = builtin->Closest(origin, direction, skip);
+			const auto b = other->Closest(origin, direction, skip);
+			ASSERT_EQ(a.triangle == UgcRays::NONE, b.triangle == UgcRays::NONE) << UgcRays::Name(backend) << " ray " << i;
+			if (a.triangle == UgcRays::NONE) continue;
+			hits++;
+			EXPECT_NEAR(a.t, b.t, 1e-4f * std::max(1.0f, a.t)) << UgcRays::Name(backend) << " ray " << i;
+			EXPECT_NE(b.triangle, skip);
+			if (a.triangle == b.triangle) {
+				sameTriangle++;
+				EXPECT_NEAR(a.u, b.u, 1e-3f);
+				EXPECT_NEAR(a.v, b.v, 1e-3f);
+			}
+			const float maxT = next() * 4.0f;
+			EXPECT_EQ(builtin->Occluded(origin, direction, 1e-4f, maxT), other->Occluded(origin, direction, 1e-4f, maxT)) << UgcRays::Name(backend) << " ray " << i;
+			// Limited: the nearest hit before maxT, or none
+			const auto limited = other->Closest(origin, direction, skip, maxT);
+			EXPECT_EQ(limited.triangle != UgcRays::NONE, b.t < maxT) << UgcRays::Name(backend) << " ray " << i;
+		}
+		EXPECT_GT(hits, 15000u);
+		EXPECT_GE(sameTriangle, hits - hits / 1000) << UgcRays::Name(backend);
+	}
+}
+
+TEST(UgcRays, OtherBackendsMakeTheSameModels) {
+	// The hidden faces and the occlusion with each backend. The paths bounce, so a hit found a rounding further away
+	// sends a path on from a slightly different point: the rare triangle decided by a path that only just gets out
+	// may go the other way. The small test model's files are the same.
+	const auto mesh = Clutter();
+	UgcHsr::Options hsr;
+	hsr.seed = 99;
+	const auto expected = UgcHsr::Visible(mesh, hsr);
+	const auto aoExpected = UgcRender::AmbientOcclusion(mesh, mesh, 2.0f, 64);
+	UgcBricks::BrickLibrary library(MakeRes(), 0);
+	const auto builtinModel = UgcJobs::ProcessModel(LOOKS_LXFML, library, SmallSettings(), 7);
+	ASSERT_TRUE(builtinModel.ok) << builtinModel.error;
+	for (const auto backend : OtherBackends()) {
+		hsr.rays = backend;
+		const auto visible = UgcHsr::Visible(mesh, hsr);
+		ASSERT_EQ(visible.size(), expected.size());
+		size_t differ = 0;
+		for (size_t t = 0; t < visible.size(); t++) differ += visible[t] != expected[t] ? 1 : 0;
+		EXPECT_LE(differ, visible.size() / 50) << UgcRays::Name(backend);
+		// The room's closed-off box is removed and what the doorway shows is kept, as with builtin
+		const auto closed = UgcHsr::Visible(Room(true), hsr);
+		for (size_t t = 0; t < 12; t++) EXPECT_FALSE(closed[t]) << UgcRays::Name(backend) << " " << t;
+		const auto open = UgcHsr::Visible(Room(false), hsr);
+		for (size_t t = 0; t < open.size(); t++) EXPECT_TRUE(open[t]) << UgcRays::Name(backend) << " " << t;
+
+		const auto ao = UgcRender::AmbientOcclusion(mesh, mesh, 2.0f, 64, backend);
+		ASSERT_EQ(ao.size(), aoExpected.size());
+		double total = 0.0;
+		for (size_t v = 0; v < ao.size(); v++) {
+			EXPECT_NEAR(ao[v], aoExpected[v], 0.05f) << UgcRays::Name(backend) << " vertex " << v;
+			total += std::abs(ao[v] - aoExpected[v]);
+		}
+		EXPECT_LT(total / static_cast<double>(ao.size()), 0.002) << UgcRays::Name(backend);
+
+		auto settings = SmallSettings();
+		settings.hsr.rays = backend;
+		settings.ao.rays = backend;
+		const auto made = UgcJobs::ProcessModel(LOOKS_LXFML, library, settings, 7);
+		ASSERT_TRUE(made.ok) << made.error;
+		EXPECT_EQ(made.files.at("model.nif.checksum"), builtinModel.files.at("model.nif.checksum")) << UgcRays::Name(backend);
+		EXPECT_EQ(made.files.at("icon.png"), builtinModel.files.at("icon.png")) << UgcRays::Name(backend);
+	}
+}