DarkflameServer/dUgcServer/Render/UgcHsr.cpp

#include "UgcHsr.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
#include "UgcRays.h"
#include "UgcRender.h"
#include "UgcThrottle.h"

namespace {
	constexpr float INF = std::numeric_limits<float>::infinity();
	constexpr float PI = 3.14159265358979f;
	// Blender's default Glossy bounces (LU Toolbox leaves it)
	constexpr int GLOSSY_BOUNCES = 4;
	// Points on one triangle at most (very big triangles get their points further apart)
	constexpr size_t MAX_POINTS = 4096;

	// PCG32 (O'Neill), seeded per path so every path is the same whatever order they're traced in
	class Random {
	public:
		explicit Random(uint64_t seed) {
			m_State = 0;
			Next();
			m_State += seed;
			Next();
		}

		uint32_t Next() {
			const uint64_t old = m_State;
			m_State = old * 6364136223846793005ull + 1442695040888963407ull;
			const auto xorshifted = static_cast<uint32_t>(((old >> 18u) ^ old) >> 27u);
			const auto rot = static_cast<uint32_t>(old >> 59u);
			return (xorshifted >> rot) | (xorshifted << ((32u - rot) & 31u));
		}

		// In [0, 1)
		float Float() { return static_cast<float>(Next() >> 8) * (1.0f / 16777216.0f); }

	private:
		uint64_t m_State;
	};

	uint64_t Mix(uint64_t h) {
		h ^= h >> 33;
		h *= 0xff51afd7ed558ccdull;
		h ^= h >> 33;
		h *= 0xc4ceb9fe1a85ec53ull;
		h ^= h >> 33;
		return h;
	}

	uint64_t PathSeed(uint64_t seed, uint64_t triangle, uint64_t point, uint64_t sample) {
		return Mix(Mix(Mix(seed ^ 0x9E3779B97F4A7C15ull) ^ triangle) ^ (point << 16 | sample));
	}

	// Directions around a unit normal (Duff et al., "Building an Orthonormal Basis, Revisited")
	void Basis(const glm::vec3& n, glm::vec3& t, glm::vec3& b) {
		const float sign = std::copysign(1.0f, n.z);
		const float a = -1.0f / (sign + n.z);
		const float c = n.x * n.y * a;
		t = glm::vec3(1.0f + sign * n.x * n.x * a, sign * c, -sign * n.x);
		b = glm::vec3(c, sign + n.y * n.y * a, -n.y);
	}

	glm::vec3 CosineHemisphere(const glm::vec3& n, float u, float v) {
		glm::vec3 t, b;
		Basis(n, t, b);
		const float r = std::sqrt(u), phi = 2.0f * PI * v;
		return t * (r * std::cos(phi)) + b * (r * std::sin(phi)) + n * std::sqrt(std::max(0.0f, 1.0f - u));
	}

	/**
	 * Cycles' ensure_valid_reflection, which the Principled BSDF applies to its normal: N turned towards Ng just
	 * enough that the mirror reflection of I (towards where the light goes) stays above the surface.
	 */
	glm::vec3 EnsureValidReflection(const glm::vec3& ng, const glm::vec3& i, const glm::vec3& n) {
		const auto r = 2.0f * glm::dot(n, i) * n - i;
		const float threshold = std::min(0.9f * glm::dot(ng, i), 0.01f);
		if (glm::dot(ng, r) >= threshold) return n;
		const float ndotng = glm::dot(n, ng);
		const auto x = glm::normalize(n - ndotng * ng);
		const float ix = glm::dot(i, x), iz = glm::dot(i, ng);
		const float ix2 = ix * ix, iz2 = iz * iz;
		const float a = ix2 + iz2;
		const float b = std::sqrt(std::max(ix2 * (a - threshold * threshold), 0.0f));
		const float c = iz * threshold + a;
		const float fac = 0.5f / a;
		const float n1z2 = fac * (b + c), n2z2 = fac * (-b + c);
		bool valid1 = n1z2 > 1e-5f && n1z2 <= 1.0f + 1e-5f;
		bool valid2 = n2z2 > 1e-5f && n2z2 <= 1.0f + 1e-5f;
		const auto root = [](float v) { return std::sqrt(std::max(v, 0.0f)); };
		glm::vec2 chosen;
		if (valid1 && valid2) {
			const glm::vec2 n1(root(1.0f - n1z2), root(n1z2)), n2(root(1.0f - n2z2), root(n2z2));
			const float r1 = 2.0f * (n1.x * ix + n1.y * iz) * n1.y - iz;
			const float r2 = 2.0f * (n2.x * ix + n2.y * iz) * n2.y - iz;
			valid1 = r1 >= 1e-5f;
			valid2 = r2 >= 1e-5f;
			if (valid1 && valid2) chosen = r1 < r2 ? n1 : n2;
			else chosen = r1 > r2 ? n1 : n2;
		} else if (valid1 || valid2) {
			const float z2 = valid1 ? n1z2 : n2z2;
			chosen = glm::vec2(root(1.0f - z2), root(z2));
		} else {
			return ng;
		}
		return chosen.x * x + chosen.y * ng;
	}

	float SchlickFresnel(float u) {
		const float m = std::clamp(1.0f - u, 0.0f, 1.0f);
		const float m2 = m * m;
		return m2 * m2 * m;
	}

	float FresnelDielectricCos(float cosi, float eta) {
		const float c = std::abs(cosi);
		float g = eta * eta - 1.0f + c * c;
		if (!(g > 0.0f)) return 1.0f;
		g = std::sqrt(g);
		const float a = (g - c) / (g + c);
		const float b = (c * (g + c) - 1.0f) / (c * (g - c) + 1.0f);
		return 0.5f * a * a * (1.0f + b * b);
	}

	/**
	 * The bake material: a new material's Principled BSDF (Blender 3.1: base color 0.8, specular 0.5, roughness 0.5,
	 * GGX) as Cycles 3.1 samples it, a diffuse closure (Burley, with retro-reflection) and a specular one (GGX with a
	 * dielectric Fresnel, IOR 1.5) picked by their sample weights. Only the directions a path takes and its throughput
	 * (for the Russian roulette) matter here.
	 */
	class Principled {
	public:
		static constexpr float BASE = 0.8f;      // Base Color
		static constexpr float ROUGHNESS = 0.5f; // Roughness
		static constexpr float IOR = 1.5f;       // (2 / (1 - sqrt(0.08 * Specular 0.5))) - 1
		static constexpr float CSPEC0 = 0.04f;   // Specular 0.5 * 0.08

		// n: the material's normal (after EnsureValidReflection), i: towards where the path came from, ng: the
		// geometric normal on the same side, sn: the shading normal. `first`: the bake point (Cycles' defensive
		// sampling there); minRayPdf: the smallest bounce pdf so far (Cycles' Filter Glossy 1 blurs the specular)
		Principled(const glm::vec3& n, const glm::vec3& i, const glm::vec3& ng, const glm::vec3& sn, bool first, float minRayPdf)
			: m_N(n), m_I(i), m_Ng(ng), m_Sn(sn) {
			m_F0 = FresnelDielectricCos(1.0f, IOR);
			m_WeightDiffuse = BASE;
			m_WeightSpecular = Fresnel(m_I, m_N); // the closure's weight 1 times its Fresnel color's average
			if (first) {
				const float sum = m_WeightDiffuse + m_WeightSpecular;
				m_WeightDiffuse = std::max(m_WeightDiffuse, 0.125f * sum);
				m_WeightSpecular = std::max(m_WeightSpecular, 0.125f * sum);
			}
			m_Alpha = std::clamp(ROUGHNESS * ROUGHNESS, 0.0f, 1.0f);
			if (minRayPdf < 1.0f) m_Alpha = std::max(std::sqrt(1.0f - minRayPdf) * 0.5f, m_Alpha);
		}

		struct Sample {
			glm::vec3 direction{};
			float pdf{};        // of the mixture; 0: nothing sampled, the path ends
			float throughput{}; // eval / pdf
			bool glossy{};
		};

		Sample Draw(Random& random) const {
			Sample sample;
			const float pick = random.Float() * (m_WeightDiffuse + m_WeightSpecular);
			const float u = random.Float(), v = random.Float();
			float pdfDiffuse = 0.0f, pdfSpecular = 0.0f;
			float eval = 0.0f;
			if (pick < m_WeightDiffuse) {
				sample.direction = CosineHemisphere(m_N, u, v);
				if (!(glm::dot(m_Ng, sample.direction) > 0.0f)) return sample;
				eval = EvalDiffuse(sample.direction, pdfDiffuse);
				if (!(pdfDiffuse > 0.0f) || !(eval > 0.0f)) return sample;
				if (!Transmission(sample.direction)) eval += EvalSpecular(sample.direction, pdfSpecular);
			} else {
				sample.glossy = true;
				if (!SampleSpecular(u, v, sample.direction, eval, pdfSpecular)) return sample;
				if (!Transmission(sample.direction)) eval += EvalDiffuse(sample.direction, pdfDiffuse);
			}
			sample.pdf = (pdfDiffuse * m_WeightDiffuse + pdfSpecular * m_WeightSpecular) / (m_WeightDiffuse + m_WeightSpecular);
			sample.throughput = sample.pdf > 0.0f ? eval / sample.pdf : 0.0f;
			return sample;
		}

	private:
		// Cycles decides reflection or transmission by the shading normal; both closures only reflect
		bool Transmission(const glm::vec3& l) const { return glm::dot(m_Sn, l) < 0.0f; }

		float Fresnel(const glm::vec3& l, const glm::vec3& h) const {
			const float fh = (FresnelDielectricCos(glm::dot(l, h), IOR) - m_F0) / (1.0f - m_F0);
			return CSPEC0 * (1.0f - fh) + fh;
		}

		// bsdf_principled_diffuse_eval_reflect times the closure weight
		float EvalDiffuse(const glm::vec3& l, float& pdf) const {
			const float nl = glm::dot(m_N, l);
			if (!(nl > 0.0f)) {
				pdf = 0.0f;
				return 0.0f;
			}
			pdf = nl / PI;
			const float nv = glm::dot(m_N, m_I);
			const float fv = SchlickFresnel(nv), fl = SchlickFresnel(nl);
			float f = (1.0f - 0.5f * fv) * (1.0f - 0.5f * fl);
			const float rr = ROUGHNESS * (glm::dot(l, m_I) + 1.0f);
			f += rr * (fl + fv + fl * fv * (rr - 1.0f));
			return BASE * nl / PI * f;
		}

		// bsdf_microfacet_ggx_eval_reflect (isotropic, Fresnel) times the closure weight 1
		float EvalSpecular(const glm::vec3& l, float& pdf) const {
			pdf = 0.0f;
			const float cosNO = glm::dot(m_N, m_I), cosNI = glm::dot(m_N, l);
			if (!(cosNI > 0.0f && cosNO > 0.0f) || m_Alpha * m_Alpha <= 1e-7f) return 0.0f;
			const auto m = glm::normalize(l + m_I);
			const float alpha2 = m_Alpha * m_Alpha;
			const float cosThetaM = glm::dot(m_N, m);
			const float cosThetaM2 = cosThetaM * cosThetaM;
			const float tanThetaM2 = (1.0f - cosThetaM2) / cosThetaM2;
			const float d = alpha2 / (PI * cosThetaM2 * cosThetaM2 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
			const float g1o = 2.0f / (1.0f + std::sqrt(std::max(0.0f, 1.0f + alpha2 * (1.0f - cosNO * cosNO) / (cosNO * cosNO))));
			const float g1i = 2.0f / (1.0f + std::sqrt(std::max(0.0f, 1.0f + alpha2 * (1.0f - cosNI * cosNI) / (cosNI * cosNI))));
			const float common = d * 0.25f / cosNO;
			pdf = g1o * common;
			return Fresnel(l, m) * g1o * g1i * common;
		}

		// bsdf_microfacet_ggx_sample (isotropic): a visible normal (Heitz and d'Eon), the view reflected in it
		bool SampleSpecular(float randu, float randv, glm::vec3& l, float& eval, float& pdf) const {
			const float cosNO = glm::dot(m_N, m_I);
			if (!(cosNO > 0.0f)) return false;
			glm::vec3 x, y;
			Basis(m_N, x, y);
			// Stretch the view, sample the slopes, rotate and unstretch
			auto local = glm::normalize(glm::vec3(m_Alpha * glm::dot(x, m_I), m_Alpha * glm::dot(y, m_I), cosNO));
			float costheta = 1.0f, sintheta = 0.0f, cosphi = 1.0f, sinphi = 0.0f;
			if (local.z < 0.99999f) {
				costheta = local.z;
				sintheta = std::sqrt(std::max(0.0f, 1.0f - costheta * costheta));
				cosphi = local.x / sintheta;
				sinphi = local.y / sintheta;
			}
			float slopeX, slopeY, g1o;
			if (costheta >= 0.99999f) {
				const float r = std::sqrt(randu / (1.0f - randu));
				const float phi = 2.0f * PI * randv;
				slopeX = r * std::cos(phi);
				slopeY = r * std::sin(phi);
				g1o = 1.0f;
			} else {
				const float tanThetaI = sintheta / costheta;
				const float g1Inv = 0.5f * (1.0f + std::sqrt(std::max(0.0f, 1.0f + tanThetaI * tanThetaI)));
				g1o = 1.0f / g1Inv;
				const float a = 2.0f * randu * g1Inv - 1.0f;
				const float aa = a * a;
				const float tmp = 1.0f / (aa - 1.0f);
				const float b = tanThetaI, bb = b * b;
				const float dd = std::sqrt(std::max(0.0f, bb * (tmp * tmp) - (aa - bb) * tmp));
				const float x1 = b * tmp - dd, x2 = b * tmp + dd;
				slopeX = (a < 0.0f || x2 * tanThetaI > 1.0f) ? x1 : x2;
				float s;
				if (randv > 0.5f) {
					s = 1.0f;
					randv = 2.0f * (randv - 0.5f);
				} else {
					s = -1.0f;
					randv = 2.0f * (0.5f - randv);
				}
				const float z = (randv * (randv * (randv * 0.27385f - 0.73369f) + 0.46341f)) / (randv * (randv * (randv * 0.093073f + 0.309420f) - 1.0f) + 0.597999f);
				slopeY = s * z * std::sqrt(std::max(0.0f, 1.0f + slopeX * slopeX));
			}
			const float rotated = cosphi * slopeX - sinphi * slopeY;
			slopeY = sinphi * slopeX + cosphi * slopeY;
			slopeX = rotated * m_Alpha;
			slopeY *= m_Alpha;
			const auto localM = glm::normalize(glm::vec3(-slopeX, -slopeY, 1.0f));
			const auto m = x * localM.x + y * localM.y + m_N * localM.z;
			const float cosMO = glm::dot(m, m_I);
			if (!(cosMO > 0.0f)) return false;
			l = 2.0f * cosMO * m - m_I;
			if (!(glm::dot(m_Ng, l) > 0.0f)) return false;
			const float alpha2 = m_Alpha * m_Alpha;
			const float cosThetaM2 = localM.z * localM.z;
			const float tanThetaM2 = 1.0f / cosThetaM2 - 1.0f;
			const float d = alpha2 / (PI * cosThetaM2 * cosThetaM2 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
			const float cosNI = glm::dot(m_N, l);
			const float g1i = 2.0f / (1.0f + std::sqrt(std::max(0.0f, 1.0f + alpha2 * (1.0f - cosNI * cosNI) / (cosNI * cosNI))));
			const float common = g1o * d * 0.25f / cosNO;
			pdf = common;
			eval = g1i * common * Fresnel(l, m);
			return pdf > 0.0f && eval > 0.0f;
		}

		glm::vec3 m_N, m_I, m_Ng, m_Sn;
		float m_F0{};
		float m_WeightDiffuse{}, m_WeightSpecular{};
		float m_Alpha{};
	};

	// 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);
		float enter = 0.0f, exit = maxT;
		for (int axis = 0; axis < 3; axis++) {
			if (std::abs(direction[axis]) < 1e-20f) {
				if (origin[axis] < MIN[axis] || origin[axis] > MAX[axis]) return INF;
				continue;
			}
			const float inv = 1.0f / direction[axis];
			float t0 = (MIN[axis] - origin[axis]) * inv, t1 = (MAX[axis] - origin[axis]) * inv;
			if (t0 > t1) std::swap(t0, t1);
			enter = std::max(enter, t0);
			exit = std::min(exit, t1);
			if (enter > exit) return INF;
		}
		return enter;
	}

	class Tracer {
	public:
		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();
		}

		// The triangle's normal from its winding (unit; zero when it has no area)
		glm::vec3 FaceNormal(uint32_t t) const {
			const auto& a = m_Mesh.positions[m_Mesh.indices[t * 3]];
			const auto n = glm::cross(m_Mesh.positions[m_Mesh.indices[t * 3 + 1]] - a, m_Mesh.positions[m_Mesh.indices[t * 3 + 2]] - a);
			const float length = glm::length(n);
			return length > 0.0f ? n / length : glm::vec3(0.0f);
		}

		// The vertex normals at barycentric w (Cycles' smooth normal: the face's own when they sum to nothing)
		glm::vec3 ShadingNormal(uint32_t t, const glm::vec3& w, const glm::vec3& faceNormal) const {
			if (!m_Smooth) return faceNormal;
			const auto n = m_Mesh.normals[m_Mesh.indices[t * 3]] * w.x + m_Mesh.normals[m_Mesh.indices[t * 3 + 1]] * w.y + m_Mesh.normals[m_Mesh.indices[t * 3 + 2]] * w.z;
			const float length = glm::length(n);
			return length > 0.0f ? n / length : faceNormal;
		}

		glm::vec3 Position(uint32_t t, const glm::vec3& w) const {
			return m_Mesh.positions[m_Mesh.indices[t * 3]] * w.x + m_Mesh.positions[m_Mesh.indices[t * 3 + 1]] * w.y + m_Mesh.positions[m_Mesh.indices[t * 3 + 2]] * w.z;
		}

		// A path between its rays
		struct Path {
			glm::vec3 ng{}, n{}, incoming{}, p{};
			glm::vec3 direction{}; // of the ray it waits for
			uint32_t self{};       // the triangle it is on (the ray never hits it)
			float throughput{ 1.0f };
			float minRayPdf{ INF };
			int glossy{};
			int bounce{};
			bool sampleGlossy{};
			Random random;
		};

		enum class eStep : uint8_t { RAY, ESCAPED, ENDED };

		// A path from point w (barycentric) of triangle t, its random numbers from `seed`
		Path Start(uint32_t t, const glm::vec3& w, uint64_t seed) const {
			Path path{ .random = Random(seed) };
			path.ng = FaceNormal(t);
			// The bake looks at the point along its smooth normal (the ray comes from there): that side is lit, and
			// when the winding faces the other way Cycles treats it as the back of the face (both normals turned)
			path.n = ShadingNormal(t, w, path.ng);
			path.incoming = path.n;
			if (glm::dot(path.ng, path.incoming) < 0.0f) {
				path.ng = -path.ng;
				path.n = -path.n;
			}
			path.p = Position(t, w);
			path.self = t;
			return path;
		}

		// The direction the path goes on in: RAY (from path.p along path.direction), or ENDED
		eStep Scatter(Path& path) const {
			// The material's normal: the Principled BSDF keeps its mirror reflection above the surface
			const Principled material(EnsureValidReflection(path.ng, path.incoming, path.n), path.incoming, path.ng, path.n, path.bounce == 0, path.minRayPdf);
			const auto sample = material.Draw(path.random);
			// Nothing sampled (below the surface): the path ends
			if (!(sample.pdf > 0.0f) || !(sample.throughput > 0.0f)) return eStep::ENDED;
			path.direction = sample.direction;
			path.throughput *= sample.throughput;
			path.minRayPdf = std::min(path.minRayPdf, sample.pdf);
			path.sampleGlossy = sample.glossy;
			return eStep::RAY;
		}

		// What the ray hit: ESCAPED (nothing: the sky), ENDED, or RAY (it bounced: Scatter again)
		eStep Bounce(Path& path, const UgcRays::Hit& hit) const {
			const auto& direction = path.direction;
			if (m_Options.groundPlane && GroundHit(path.p, direction, hit.t) < hit.t) return eStep::ENDED;
			if (hit.triangle == UgcRays::NONE) return eStep::ESCAPED;
			// Past the bounce limits the next surface doesn't scatter (Cycles: Max Bounces, Glossy 4)
			if (path.bounce + 1 > m_Options.bounces) return eStep::ENDED;
			if (path.sampleGlossy && ++path.glossy > GLOSSY_BOUNCES) return eStep::ENDED;
			if (path.bounce + 1 >= 2) {
				const float probability = std::min(std::sqrt(path.throughput), 1.0f);
				if (path.random.Float() >= probability) return eStep::ENDED;
				path.throughput /= probability;
			}
			path.self = hit.triangle;
			path.ng = FaceNormal(path.self);
			path.n = ShadingNormal(path.self, glm::vec3(1.0f - hit.u - hit.v, hit.u, hit.v), path.ng);
			path.incoming = -direction;
			// Hit from behind: the back is lit (both normals turned towards the ray)
			if (glm::dot(path.ng, path.incoming) < 0.0f) {
				path.ng = -path.ng;
				path.n = -path.n;
			}
			path.p = path.p + direction * hit.t;
			path.bounce++;
			return eStep::RAY;
		}

		/**
		 * One path from point w (barycentric) of triangle t, as a Cycles diffuse bake traces it: whether it reaches the
		 * sky. The path bounces off the model in the directions the bake material draws until a bounce ray hits
		 * nothing, which is the sky. The sky isn't sampled directly: Cycles doesn't sample a world of one flat color as
		 * a light, so only rays that happen to leave count. Up to `bounces` bounces (4 of them glossy), and from the
		 * second bounce on the path may end early (Cycles' Russian roulette: it goes on with probability
		 * sqrt(throughput)).
		 */
		bool Escapes(uint32_t t, const glm::vec3& w, uint64_t seed) const {
			auto path = Start(t, w, seed);
			while (Scatter(path) == eStep::RAY) {
				const auto step = Bounce(path, m_Rays->Closest(path.p, path.direction, path.self));
				if (step != eStep::RAY) return step == eStep::ESCAPED;
			}
			return false;
		}

		const UgcRays::Scene& Rays() const { return *m_Rays; }

	private:
		const UgcModel::Mesh& m_Mesh;
		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{};
	};
}

namespace UgcHsr {
	std::string_view Name(eMethod method) {
		return method == eMethod::FAST ? "fast" : "toolbox";
	}

	std::optional<eMethod> Parse(std::string_view name) {
		for (const auto method : { eMethod::TOOLBOX, eMethod::FAST }) {
			if (Name(method) == name) return method;
		}
		return std::nullopt;
	}

	std::vector<glm::vec3> SamplePoints(const glm::vec3& a, const glm::vec3& b, const glm::vec3& c, float spacing, size_t minimum) {
		// Too small to lay out: the centre and one point towards each corner (the centres of its four halved-side
		// triangles)
		const std::vector<glm::vec3> fallback{ { 1.0f / 3, 1.0f / 3, 1.0f / 3 }, { 2.0f / 3, 1.0f / 6, 1.0f / 6 }, { 1.0f / 6, 2.0f / 3, 1.0f / 6 },
			{ 1.0f / 6, 1.0f / 6, 2.0f / 3 } };
		const std::array<glm::vec3, 3> corners{ a, b, c };
		// The longest side from corner `first` to the next; the third is the apex
		int first = 0;
		float longest = -1.0f;
		for (int i = 0; i < 3; i++) {
			const float length = glm::length(corners[(i + 1) % 3] - corners[i]);
			if (length > longest) {
				longest = length;
				first = i;
			}
		}
		const float area = 0.5f * glm::length(glm::cross(b - a, c - a));
		if (!(longest > 0.0f) || !(area > 0.0f) || !(spacing > 0.0f) || !std::isfinite(area)) return fallback;
		minimum = std::min(minimum, MAX_POINTS);
		// Very big triangles: points further apart, at most MAX_POINTS
		spacing = std::max(spacing, std::sqrt(area / static_cast<float>(MAX_POINTS)));
		const float height = 2.0f * area / longest;
		std::vector<glm::vec3> points;
		for (int attempt = 0; attempt < 16; attempt++) {
			points.clear();
			const int along = std::max(1, static_cast<int>(std::ceil(longest / spacing)));
			const int rows = std::max(1, static_cast<int>(std::ceil(height / spacing)));
			for (int row = 0; row < rows; row++) {
				// v: from the longest side (0) to the apex (1); the row is (1 - v) of the side's length
				const float v = (row + 0.5f) / rows;
				const int count = std::max(1, static_cast<int>(std::lround(along * (1.0f - v))));
				for (int j = 0; j < count; j++) {
					const float u = (j + 0.5f) / count;
					glm::vec3 w{};
					w[first] = (1.0f - u) * (1.0f - v);
					w[(first + 1) % 3] = u * (1.0f - v);
					w[(first + 2) % 3] = v;
					points.push_back(w);
				}
			}
			if (points.size() >= minimum) break;
			// Fewer than the minimum: closer together
			spacing *= 0.95f * std::sqrt(static_cast<float>(points.size()) / static_cast<float>(minimum));
		}
		return points.size() < fallback.size() ? fallback : points;
	}

	std::vector<bool> Visible(const UgcModel::Mesh& mesh, const Options& options, uint64_t* pointCount, uint64_t* pathCount) {
		const size_t triangles = mesh.TriangleCount();
		std::vector<bool> visible(triangles, true);
		if (triangles == 0) return visible;
		const Tracer tracer(mesh, options);
		const int samples = std::max(options.samples, 1);
		uint64_t points = 0, paths = 0, sinceCheckpoint = 0;
		if (tracer.Rays().PrefersBatches() || options.sideBySide) {
			// Side by side (a GPU): the triangles in groups of about GROUP_PATHS paths; in each round every point of every
			// triangle of the group not seen yet starts a path, and the round's paths are traced together a bounce at a
			// time. Each path is the one traced one by one below (the same random numbers), and a triangle is kept when
			// any of its paths escapes, so the triangles decided are the same; only paths the one by one tracing would
			// have skipped after one escaped are traced too.
			constexpr size_t GROUP_PATHS = 1u << 18;
			uint32_t next = 0;
			while (next < triangles) {
				std::vector<uint32_t> group;
				std::vector<std::vector<glm::vec3>> groupWeights;
				size_t groupPaths = 0;
				for (; next < triangles && groupPaths < GROUP_PATHS; next++) {
					// A triangle without area draws nothing (LU Toolbox's bake leaves it dark too): removed
					if (tracer.FaceNormal(next) == glm::vec3(0.0f)) {
						visible[next] = false;
						continue;
					}
					const auto& a = mesh.positions[mesh.indices[next * 3]];
					const auto& b = mesh.positions[mesh.indices[next * 3 + 1]];
					const auto& c = mesh.positions[mesh.indices[next * 3 + 2]];
					groupWeights.push_back(SamplePoints(a, b, c, options.spacing, static_cast<size_t>(std::max(options.minPoints, 1))));
					group.push_back(next);
					groupPaths += groupWeights.back().size();
					points += groupWeights.back().size();
				}
				std::vector<uint8_t> escaped(group.size(), 0);
				for (int sample = 0; sample < samples; sample++) {
					std::vector<Tracer::Path> live;
					std::vector<uint32_t> owner; // the path's triangle in the group
					for (uint32_t g = 0; g < group.size(); g++) {
						if (escaped[g]) continue;
						for (size_t i = 0; i < groupWeights[g].size(); i++) {
							auto path = tracer.Start(group[g], groupWeights[g][i], PathSeed(options.seed, group[g], i, static_cast<uint64_t>(sample)));
							paths++;
							if (tracer.Scatter(path) != Tracer::eStep::RAY) continue;
							live.push_back(std::move(path));
							owner.push_back(g);
						}
					}
					if (live.empty()) break;
					std::vector<UgcRays::Ray> rays;
					std::vector<UgcRays::Hit> hits;
					while (!live.empty()) {
						UgcThrottle::Checkpoint();
						rays.resize(live.size());
						hits.resize(live.size());
						for (size_t k = 0; k < live.size(); k++) rays[k] = UgcRays::Ray{ live[k].p, 0.0f, live[k].direction, INF, live[k].self };
						tracer.Rays().Closest(rays.data(), hits.data(), live.size());
						size_t kept = 0;
						for (size_t k = 0; k < live.size(); k++) {
							if (escaped[owner[k]]) continue;
							const auto step = tracer.Bounce(live[k], hits[k]);
							if (step == Tracer::eStep::ESCAPED) escaped[owner[k]] = 1;
							if (step != Tracer::eStep::RAY || tracer.Scatter(live[k]) != Tracer::eStep::RAY) continue;
							if (kept != k) {
								live[kept] = std::move(live[k]);
								owner[kept] = owner[k];
							}
							kept++;
						}
						live.erase(live.begin() + static_cast<std::ptrdiff_t>(kept), live.end());
						owner.resize(kept);
					}
				}
				for (size_t g = 0; g < group.size(); g++) visible[group[g]] = escaped[g] != 0;
			}
			if (pointCount) *pointCount = points;
			if (pathCount) *pathCount = paths;
			return visible;
		}
		for (uint32_t t = 0; t < triangles; 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]];
			// A triangle without area draws nothing (LU Toolbox's bake leaves it dark too): removed
			if (tracer.FaceNormal(t) == glm::vec3(0.0f)) {
				visible[t] = false;
				continue;
			}
			const auto weights = SamplePoints(a, b, c, options.spacing, static_cast<size_t>(std::max(options.minPoints, 1)));
			points += weights.size();
			bool escaped = false;
			// A path from each point, then another from each, ...: a triangle that's seen is usually known at once
			for (int sample = 0; sample < samples && !escaped; sample++) {
				for (size_t i = 0; i < weights.size(); i++) {
					paths++;
					if (tracer.Escapes(t, weights[i], PathSeed(options.seed, t, i, static_cast<uint64_t>(sample)))) {
						escaped = true;
						break;
					}
				}
				sinceCheckpoint += weights.size();
				if (sinceCheckpoint >= 256) {
					UgcThrottle::Checkpoint();
					sinceCheckpoint = 0;
				}
			}
			visible[t] = escaped;
		}
		if (pointCount) *pointCount = points;
		if (pathCount) *pathCount = paths;
		return visible;
	}

	Result RemoveHiddenFaces(UgcModel::Model& model, const Options& options) {
		Result result;
		auto& opaque = model.opaque;
		result.trianglesBefore = opaque.TriangleCount() + model.transparent.TriangleCount();
		if (opaque.Empty() || !options.enabled) return result;
		result.kept = options.method == eMethod::FAST ? UgcRender::VisibleFromAround(model, options.fastResolution, options.groundPlane) :
			Visible(opaque, options, &result.points, &result.paths);
		for (const bool kept : result.kept) result.trianglesRemoved += kept ? 0 : 1;
		UgcModel::KeepTriangles(opaque, result.kept);
		return result;
	}
}