entity bvh with exact old tie order and allocation-free entity/element intersectors, entity benchmark harness and perf test scenes

src/main/java/eu/mhsl/minecraft/pixelpics/assets/model/Element.java

@@ -19,6 +19,10 @@ public final class Element {
     public final double rotAngleRad; // radians
     public final boolean rescale;
 
+    // cos/sin of the element rotation (+angle: local -> block) and of its inverse (-angle: block -> local),
+    // precomputed for the intersector's hot path.
+    public final double rotCos, rotSin, invRotCos, invRotSin;
+
     // Blockstate rotation: the model-space element and the 90° steps (x first, then y) that turn it into this
     // world-space one. Intersectors hit the unrotated element with the inverse-rotated ray, so textures stay
     // bound to the model (vanilla); under uvlock its faces carry a UV transform keeping them world-aligned.
@@ -59,6 +63,10 @@ public final class Element {
         this.rotAxis = rotAxis;
         this.rotAngleRad = rotAngleRad;
         this.rescale = rescale;
+        this.rotCos = Math.cos(rotAngleRad);
+        this.rotSin = Math.sin(rotAngleRad);
+        this.invRotCos = Math.cos(-rotAngleRad);
+        this.invRotSin = Math.sin(-rotAngleRad);
     }
 
     /** Model space -> world: vanilla blockstate rotation (x steps, then y steps; clockwise about the centre). */

src/main/java/eu/mhsl/minecraft/pixelpics/render/entity/EntityIntersector.java

@@ -10,7 +10,9 @@ import org.bukkit.util.Vector;
 /**
  * Intersects a world-space ray with a single {@link EntityCube} (oriented box). The ray is mapped
  * into the cube's local frame, slab-tested, and the entry face is sampled. Fully transparent texels
- * are treated as holes (alpha cutout); partially transparent ones keep their alpha for the tracer to blend. The returned {@code t} is a world-space distance.
+ * are treated as holes (alpha cutout); partially transparent ones keep their alpha for the tracer to blend.
+ * The {@code t} is a world-space distance. The test writes into a caller-owned {@link Hit} and allocates
+ * nothing; {@link #toFaceHit} builds the result object only for the hit the caller finally keeps.
  */
 public final class EntityIntersector {
 
@@ -18,21 +20,48 @@ public final class EntityIntersector {
 
     private EntityIntersector() {}
 
+    /** Mutable per-thread result/scratch of {@link #intersect}; never shared between threads. */
+    public static final class Hit {
+        public double t;
+        public int color;
+        int axis;
+        boolean neg;
+        final double[] o = new double[3], d = new double[3];
+    }
+
     /**
      * Nearest drawn face of the cube with {@code t > minT} (entry face first, else the far face); {@code minT}
-     * lets the tracer continue behind a translucent hit.
+     * lets the tracer continue behind a translucent hit. Cubes whose entry lies beyond {@code maxT} are
+     * rejected early (no texture lookup); a hit at exactly {@code maxT} is still reported. On a hit the face
+     * is left in {@code out}.
      */
-    public static FaceHit intersect(
-            EntityCube cube, double ox, double oy, double oz, double dx, double dy, double dz, double minT) {
-        double[] o = cube.toLocal.apply(ox, oy, oz);
-        double[] d = cube.toLocal.applyLinear(dx, dy, dz);
+    public static boolean intersect(
+            EntityCube cube,
+            double ox,
+            double oy,
+            double oz,
+            double dx,
+            double dy,
+            double dz,
+            double minT,
+            double maxT,
+            Hit out) {
+        double[] r = cube.toLocal.r, tr = cube.toLocal.t;
+        double[] o = out.o, d = out.d;
+        o[0] = r[0] * ox + r[1] * oy + r[2] * oz + tr[0];
+        o[1] = r[3] * ox + r[4] * oy + r[5] * oz + tr[1];
+        o[2] = r[6] * ox + r[7] * oy + r[8] * oz + tr[2];
+        d[0] = r[0] * dx + r[1] * dy + r[2] * dz;
+        d[1] = r[3] * dx + r[4] * dy + r[5] * dz;
+        d[2] = r[6] * dx + r[7] * dy + r[8] * dz;
+        double[] from = cube.from, to = cube.to;
 
         double tmin = Double.NEGATIVE_INFINITY, tmax = Double.POSITIVE_INFINITY;
         int inAxis = -1, outAxis = -1;
         boolean inNeg = false, outNeg = false;
         for (int a = 0; a < 3; a++) {
             if (Math.abs(d[a]) < EPS) {
-                if (o[a] < cube.from[a] - EPS || o[a] > cube.to[a] + EPS) return null;
+                if (o[a] < from[a] - EPS || o[a] > to[a] + EPS) return false;
                 continue;
             }
             double inv = 1.0 / d[a];
@@ -40,8 +69,8 @@ public final class EntityIntersector {
             // max face. Decided by the direction, not by comparing t1/t2: for a zero-thickness box
             // (wing membranes, fins) both are equal and the comparison can't tell the sides apart.
             boolean forward = inv > 0;
-            double tNear = ((forward ? cube.from[a] : cube.to[a]) - o[a]) * inv;
-            double tFar = ((forward ? cube.to[a] : cube.from[a]) - o[a]) * inv;
+            double tNear = ((forward ? from[a] : to[a]) - o[a]) * inv;
+            double tFar = ((forward ? to[a] : from[a]) - o[a]) * inv;
             if (tNear > tmin) {
                 tmin = tNear;
                 inAxis = a;
@@ -52,58 +81,54 @@ public final class EntityIntersector {
                 outAxis = a;
                 outNeg = !forward;
             }
-            if (tmin > tmax + EPS) return null;
+            if (tmin > tmax + EPS) return false;
         }
         double lo = Math.max(EPS, minT);
-        if (inAxis < 0 || tmax <= lo) return null;
+        if (inAxis < 0 || tmax <= lo || tmin > maxT) return false;
 
         // Vanilla draws entity models without back-face culling: a ray looking through a transparent
         // texel of the entry face still sees the inside of the far face (and both sides of a plane).
-        if (tmin > lo) {
-            FaceHit hit = sample(cube, o, d, tmin, inAxis, inNeg, ox, oy, oz, dx, dy, dz);
-            if (hit != null) return hit;
-        }
-        return outAxis < 0 ? null : sample(cube, o, d, tmax, outAxis, outNeg, ox, oy, oz, dx, dy, dz);
+        if (tmin > lo && sample(cube, o, d, tmin, inAxis, inNeg, out)) return true;
+        return outAxis >= 0 && sample(cube, o, d, tmax, outAxis, outNeg, out);
     }
 
-    private static FaceHit sample(
-            EntityCube cube,
-            double[] o,
-            double[] d,
-            double t,
-            int axis,
-            boolean neg,
-            double ox,
-            double oy,
-            double oz,
-            double dx,
-            double dy,
-            double dz) {
-        double px = o[0] + d[0] * t, py = o[1] + d[1] * t, pz = o[2] + d[2] * t;
-        Direction dir =
-                switch (axis) {
-                    case 0 -> neg ? Direction.WEST : Direction.EAST;
-                    case 1 -> neg ? Direction.DOWN : Direction.UP;
-                    default -> neg ? Direction.NORTH : Direction.SOUTH;
-                };
+    private static boolean sample(EntityCube cube, double[] o, double[] d, double t, int axis, boolean neg, Hit out) {
+        Direction dir = direction(axis, neg);
         Face face = cube.faces[dir.ordinal()];
-        if (face == null) return null;
-
+        if (face == null) return false;
+        double px = o[0] + d[0] * t, py = o[1] + d[1] * t, pz = o[2] + d[2] * t;
         double fx = frac(px, cube.from[0], cube.to[0]);
         double fy = frac(py, cube.from[1], cube.to[1]);
         double fz = frac(pz, cube.from[2], cube.to[2]);
         // Vanilla face-UV orientation, shared by block-model faces and entity box-UV (see BoxUv).
-        double[] st = ElementIntersector.faceUv(dir, fx, fy, fz);
-        double s = st[0], tt = st[1];
-        int color = face.sample(s, tt);
-        if (!ColorUtil.opaque(color)) return null;
+        int color = face.sample(ElementIntersector.faceS(dir, fx, fy, fz), ElementIntersector.faceT(dir, fy, fz));
+        if (!ColorUtil.opaque(color)) return false;
+        out.t = t;
+        out.color = color;
+        out.axis = axis;
+        out.neg = neg;
+        return true;
+    }
 
+    /** The result object for a hit of {@code cube} recorded by {@link #intersect} (world point + normal). */
+    public static FaceHit toFaceHit(
+            EntityCube cube, Hit h, double ox, double oy, double oz, double dx, double dy, double dz) {
+        double t = h.t;
+        Direction dir = direction(h.axis, h.neg);
         Vector world = new Vector(ox + dx * t, oy + dy * t, oz + dz * t);
         double[] n = cube.toWorld.applyLinear(dir.nx, dir.ny, dir.nz);
         Vector normal = new Vector(n[0], n[1], n[2]).normalize();
         // A far/back face is seen from inside: shade it as facing the viewer.
         if (normal.getX() * dx + normal.getY() * dy + normal.getZ() * dz > 0) normal.multiply(-1);
-        return new FaceHit(t, world, normal, color, -1, cube.additive);
+        return new FaceHit(t, world, normal, h.color, -1, cube.additive);
+    }
+
+    private static Direction direction(int axis, boolean neg) {
+        return switch (axis) {
+            case 0 -> neg ? Direction.WEST : Direction.EAST;
+            case 1 -> neg ? Direction.DOWN : Direction.UP;
+            default -> neg ? Direction.NORTH : Direction.SOUTH;
+        };
     }
 
     private static double frac(double v, double lo, double hi) {

src/main/java/eu/mhsl/minecraft/pixelpics/render/entity/EntityScene.java

@@ -7,17 +7,63 @@ import java.util.ArrayList;
 import java.util.List;
 
 /**
- * All baked entities (mobs + block-entities) for one render. Provides the nearest hit along a ray,
- * using per-entity and per-cube AABB broad-phase culling. Block-entities are baked into the same list
- * and tested identically, so the raytracer needs no special case. Immutable after construction → safe
- * for the parallel tracer.
+ * All baked entities (mobs + block-entities + decorations + objects) for one render. Provides the nearest hit
+ * along a ray. Every entity cube goes into one bounding volume hierarchy over the cubes' world AABBs (binned
+ * SAH, flat arrays), traversed front to back with the current best hit as the far bound, so a ray only tests
+ * the cubes near it instead of every entity — farms with hundreds of mobs/items cost about as much as a few.
+ * Near-ties (coplanar faces of different cubes) are decided by replaying the former linear scan over the few
+ * candidates, so images stay pixel-identical. Immutable after construction; the per-ray scratch is per thread
+ * → safe for the parallel tracer.
  */
 public final class EntityScene {
 
     private static final double EPS = 1e-7;
-    private final List<RenderedEntity> entities;
+    private static final int LEAF_SIZE = 4;
+    private static final int SAH_BINS = 16;
+    private static final int MAX_SAH_DEPTH = 40;
     private static volatile boolean warned;
 
+    private final List<RenderedEntity> entities;
+
+    // Cubes in BVH leaf order, their bake-order rank (tie-break) and world AABBs (minX..maxZ, 6 per cube).
+    private final EntityCube[] cubes;
+    private final int[] cubeRank;
+    private final double[] cubeBox;
+
+    // BVH nodes: box (6 per node); inner nodes have left child = node + 1 and right child = nodeRight,
+    // leaves (nodeCount > 0) cover cubes[nodeStart .. nodeStart + nodeCount).
+    private double[] nodeBox;
+    private int[] nodeRight, nodeStart, nodeCount;
+    private int nodes;
+    private final int maxDepth;
+
+    // Hits within this distance of the nearest one are candidates for the exact tie-break (pass 2 of
+    // nearestHit); far above double rounding at render distances, far below anything visible.
+    private static final double TIE_SLACK = 1e-6;
+
+    /** Per-thread traversal scratch: stack and the candidate hits (reused, grown on demand). */
+    private static final class Traversal {
+        int[] stack = new int[64];
+        double[] stackT = new double[64];
+        EntityIntersector.Hit[] cand = new EntityIntersector.Hit[0];
+        int[] candCube = new int[8];
+        int[] order = new int[8];
+
+        EntityIntersector.Hit candidate(int k) {
+            if (k == cand.length) {
+                int len = Math.max(8, k * 2);
+                EntityIntersector.Hit[] c = java.util.Arrays.copyOf(cand, len);
+                for (int i = k; i < len; i++) c[i] = new EntityIntersector.Hit();
+                cand = c;
+                candCube = java.util.Arrays.copyOf(candCube, len);
+                order = java.util.Arrays.copyOf(order, len);
+            }
+            return cand[k];
+        }
+    }
+
+    private static final ThreadLocal<Traversal> TRAVERSAL = ThreadLocal.withInitial(Traversal::new);
+
     public EntityScene(
             List<EntityState> states,
             CemBaker baker,
@@ -32,6 +78,17 @@ public final class EntityScene {
         addAll(blockEntities, beBaker);
         addAll(decorations, decoBaker);
         addAll(objects, objectBaker);
+
+        int n = 0;
+        for (RenderedEntity e : entities) n += e.cubes.size();
+        EntityCube[] all = new EntityCube[n];
+        int k = 0;
+        for (RenderedEntity e : entities) for (EntityCube c : e.cubes) all[k++] = c;
+
+        cubes = new EntityCube[n];
+        cubeRank = new int[n];
+        cubeBox = new double[n * 6];
+        maxDepth = build(all);
     }
 
     /** Bakes each state and keeps the ones with geometry. A null baker (path disabled) is skipped. */
@@ -61,52 +118,365 @@ public final class EntityScene {
     /** Nearest entity hit with {@code minT < t < maxT}, or null. */
     public FaceHit nearestHit(
             double ox, double oy, double oz, double dx, double dy, double dz, double minT, double maxT) {
-        FaceHit best = null;
-        double bestT = maxT;
-        for (RenderedEntity e : entities) {
-            if (!rayAabb(e.aabbMin, e.aabbMax, ox, oy, oz, dx, dy, dz, minT, bestT)) continue;
-            for (EntityCube cube : e.cubes) {
-                if (!rayAabb(cube.aabbMin, cube.aabbMax, ox, oy, oz, dx, dy, dz, minT, bestT)) continue;
-                FaceHit hit = EntityIntersector.intersect(cube, ox, oy, oz, dx, dy, dz, minT);
-                if (hit != null && hit.t() > EPS && hit.t() > minT && hit.t() < bestT) {
-                    best = hit;
-                    bestT = hit.t();
+        if (nodes == 0) return null;
+        Traversal tr = TRAVERSAL.get();
+        if (tr.stack.length < maxDepth + 2) {
+            tr.stack = new int[maxDepth + 2];
+            tr.stackT = new double[maxDepth + 2];
+        }
+        int[] stack = tr.stack;
+        double[] stackT = tr.stackT;
+
+        // Per-axis slab setup, the same arithmetic as a per-box test (inv = 1/d; |d| < EPS = parallel).
+        boolean px = Math.abs(dx) < EPS, py = Math.abs(dy) < EPS, pz = Math.abs(dz) < EPS;
+        double ix = px ? 0 : 1.0 / dx, iy = py ? 0 : 1.0 / dy, iz = pz ? 0 : 1.0 / dz;
+        double t0 = Math.max(0, minT);
+
+        // Pass 1 (BVH, front to back): every cube hit within TIE_SLACK of the nearest one becomes a candidate.
+        // The far bound stays 2·TIE_SLACK behind the nearest hit so rounding in the AABB tests can't drop one.
+        double nearest = Double.POSITIVE_INFINITY, bound = maxT;
+        int n = 0;
+        double[] nb = nodeBox, cb = cubeBox;
+        int sp = 0;
+        double rootT = entry(nb, 0, ox, oy, oz, ix, iy, iz, px, py, pz, t0, bound);
+        if (rootT != rootT) return null;
+        stack[sp] = 0;
+        stackT[sp++] = rootT;
+        while (sp > 0) {
+            int node = stack[--sp];
+            if (stackT[sp] > bound) continue; // the nearest hit got closer since this node was queued
+            int count = nodeCount[node];
+            if (count > 0) {
+                for (int i = nodeStart[node], end = i + count; i < end; i++) {
+                    double ct = entry(cb, i, ox, oy, oz, ix, iy, iz, px, py, pz, t0, bound);
+                    if (ct != ct) continue; // NaN: AABB missed
+                    EntityIntersector.Hit hit = tr.candidate(n);
+                    if (!EntityIntersector.intersect(cubes[i], ox, oy, oz, dx, dy, dz, minT, bound, hit)) continue;
+                    double t = hit.t;
+                    if (!(t < maxT) || t > nearest + TIE_SLACK) continue;
+                    tr.candCube[n++] = i;
+                    if (t < nearest) {
+                        nearest = t;
+                        bound = Math.min(maxT, t + 2 * TIE_SLACK);
+                    }
+                }
+            } else {
+                int l = node + 1, r = nodeRight[node];
+                double tl = entry(nb, l, ox, oy, oz, ix, iy, iz, px, py, pz, t0, bound);
+                double trr = entry(nb, r, ox, oy, oz, ix, iy, iz, px, py, pz, t0, bound);
+                boolean hl = tl == tl, hr = trr == trr;
+                if (hl && hr) {
+                    // Near child on top of the stack: visited first, so the far one is usually culled.
+                    if (tl <= trr) {
+                        stack[sp] = r;
+                        stackT[sp++] = trr;
+                        stack[sp] = l;
+                        stackT[sp++] = tl;
+                    } else {
+                        stack[sp] = l;
+                        stackT[sp++] = tl;
+                        stack[sp] = r;
+                        stackT[sp++] = trr;
+                    }
+                } else if (hl) {
+                    stack[sp] = l;
+                    stackT[sp++] = tl;
+                } else if (hr) {
+                    stack[sp] = r;
+                    stackT[sp++] = trr;
                 }
             }
         }
-        return best;
+        if (n == 0) return null;
+
+        // Pass 2: replay the former linear scan over the near candidates in bake order — a cube counts only if
+        // its AABB passes with the best t so far as far bound and it is strictly nearer. Coplanar faces of
+        // different cubes (equal t up to rounding) thus resolve exactly as before: pixel-identical images.
+        int[] candCube = tr.candCube;
+        int[] order = tr.order;
+        int m = 0;
+        for (int k = 0; k < n; k++) {
+            if (tr.cand[k].t > nearest + TIE_SLACK) continue;
+            // Insertion by bake rank (a handful of candidates at most).
+            int rank = cubeRank[candCube[k]], j = m++;
+            while (j > 0 && cubeRank[candCube[order[j - 1]]] > rank) {
+                order[j] = order[j - 1];
+                j--;
+            }
+            order[j] = k;
+        }
+        double bestT = maxT;
+        int best = -1;
+        for (int j = 0; j < m; j++) {
+            int k = order[j];
+            double t = tr.cand[k].t;
+            double e = entry(cb, candCube[k], ox, oy, oz, ix, iy, iz, px, py, pz, t0, bestT);
+            if (e == e && t < bestT) {
+                bestT = t;
+                best = k;
+            }
+        }
+        return EntityIntersector.toFaceHit(cubes[candCube[best]], tr.cand[best], ox, oy, oz, dx, dy, dz);
     }
 
-    private static boolean rayAabb(
-            double[] min,
-            double[] max,
+    /**
+     * Slab test of box {@code i} of {@code boxes} (6 doubles per box) against the ray within {@code [t0, tMax]};
+     * returns the entry t, or NaN on a miss. Same comparisons as the former per-entity/per-cube AABB check
+     * (a box touching the interval at {@code tMax} still counts).
+     */
+    private static double entry(
+            double[] boxes,
+            int i,
             double ox,
             double oy,
             double oz,
-            double dx,
-            double dy,
-            double dz,
-            double minT,
-            double maxT) {
-        double tmin = Math.max(0, minT), tmax = maxT;
-        double[] o = {ox, oy, oz}, d = {dx, dy, dz};
-        for (int a = 0; a < 3; a++) {
-            if (Math.abs(d[a]) < EPS) {
-                if (o[a] < min[a] || o[a] > max[a]) return false;
-            } else {
-                double inv = 1.0 / d[a];
-                double t1 = (min[a] - o[a]) * inv;
-                double t2 = (max[a] - o[a]) * inv;
-                if (t1 > t2) {
-                    double tmp = t1;
-                    t1 = t2;
-                    t2 = tmp;
-                }
-                if (t1 > tmin) tmin = t1;
-                if (t2 < tmax) tmax = t2;
-                if (tmin > tmax) return false;
+            double ix,
+            double iy,
+            double iz,
+            boolean px,
+            boolean py,
+            boolean pz,
+            double t0,
+            double tMax) {
+        int b = i * 6;
+        double tmin = t0, tmax = tMax;
+        if (px) {
+            if (ox < boxes[b] || ox > boxes[b + 3]) return Double.NaN;
+        } else {
+            double t1 = (boxes[b] - ox) * ix, t2 = (boxes[b + 3] - ox) * ix;
+            if (t1 > t2) {
+                double tmp = t1;
+                t1 = t2;
+                t2 = tmp;
+            }
+            if (t1 > tmin) tmin = t1;
+            if (t2 < tmax) tmax = t2;
+            if (tmin > tmax) return Double.NaN;
+        }
+        if (py) {
+            if (oy < boxes[b + 1] || oy > boxes[b + 4]) return Double.NaN;
+        } else {
+            double t1 = (boxes[b + 1] - oy) * iy, t2 = (boxes[b + 4] - oy) * iy;
+            if (t1 > t2) {
+                double tmp = t1;
+                t1 = t2;
+                t2 = tmp;
+            }
+            if (t1 > tmin) tmin = t1;
+            if (t2 < tmax) tmax = t2;
+            if (tmin > tmax) return Double.NaN;
+        }
+        if (pz) {
+            if (oz < boxes[b + 2] || oz > boxes[b + 5]) return Double.NaN;
+        } else {
+            double t1 = (boxes[b + 2] - oz) * iz, t2 = (boxes[b + 5] - oz) * iz;
+            if (t1 > t2) {
+                double tmp = t1;
+                t1 = t2;
+                t2 = tmp;
+            }
+            if (t1 > tmin) tmin = t1;
+            if (t2 < tmax) tmax = t2;
+            if (tmin > tmax) return Double.NaN;
+        }
+        return tmin;
+    }
+
+    // ---- BVH construction (binned SAH over the cube AABB centroids) ----
+
+    private double[] centroid; // 3 per original cube
+    private double[] srcBox; // 6 per original cube
+    private int[] perm; // original cube indices, partitioned in place
+    private int filled; // cubes written in leaf order
+
+    private int build(EntityCube[] all) {
+        int n = all.length;
+        if (n == 0) {
+            nodeBox = new double[0];
+            nodeRight = nodeStart = nodeCount = new int[0];
+            return 0;
+        }
+        srcBox = new double[n * 6];
+        centroid = new double[n * 3];
+        perm = new int[n];
+        for (int i = 0; i < n; i++) {
+            EntityCube c = all[i];
+            for (int a = 0; a < 3; a++) {
+                srcBox[i * 6 + a] = c.aabbMin[a];
+                srcBox[i * 6 + 3 + a] = c.aabbMax[a];
+                centroid[i * 3 + a] = (c.aabbMin[a] + c.aabbMax[a]) * 0.5;
+            }
+            perm[i] = i;
+        }
+        int cap = 2 * n;
+        nodeBox = new double[cap * 6];
+        nodeRight = new int[cap];
+        nodeStart = new int[cap];
+        nodeCount = new int[cap];
+        int depth = buildNode(all, 0, n, 0);
+        srcBox = centroid = null;
+        perm = null;
+        return depth;
+    }
+
+    /** Builds the subtree over perm[lo, hi); returns its depth. */
+    private int buildNode(EntityCube[] all, int lo, int hi, int depth) {
+        int node = nodes++;
+        int nb = node * 6;
+        nodeBox[nb] = nodeBox[nb + 1] = nodeBox[nb + 2] = Double.MAX_VALUE;
+        nodeBox[nb + 3] = nodeBox[nb + 4] = nodeBox[nb + 5] = -Double.MAX_VALUE;
+        double[] cmin = {Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE};
+        double[] cmax = {-Double.MAX_VALUE, -Double.MAX_VALUE, -Double.MAX_VALUE};
+        for (int i = lo; i < hi; i++) {
+            int c = perm[i];
+            for (int a = 0; a < 3; a++) {
+                nodeBox[nb + a] = Math.min(nodeBox[nb + a], srcBox[c * 6 + a]);
+                nodeBox[nb + 3 + a] = Math.max(nodeBox[nb + 3 + a], srcBox[c * 6 + 3 + a]);
+                cmin[a] = Math.min(cmin[a], centroid[c * 3 + a]);
+                cmax[a] = Math.max(cmax[a], centroid[c * 3 + a]);
             }
         }
-        return true;
+        int count = hi - lo;
+        int axis = 0;
+        for (int a = 1; a < 3; a++) if (cmax[a] - cmin[a] > cmax[axis] - cmin[axis]) axis = a;
+        double extent = cmax[axis] - cmin[axis];
+        int mid;
+        if (count <= LEAF_SIZE || !(extent > 1e-9)) {
+            mid = -1;
+        } else if (depth >= MAX_SAH_DEPTH) {
+            // Pathological centroid distributions can make SAH splits very lopsided: median keeps it bounded.
+            sortByCentroid(lo, hi, axis);
+            mid = lo + count / 2;
+        } else {
+            mid = sahSplit(nb, lo, hi, axis, cmin[axis], extent);
+        }
+        if (mid < 0) {
+            // Leaf: cubes in bake order inside the leaf (ascending original index).
+            java.util.Arrays.sort(perm, lo, hi);
+            nodeStart[node] = filled;
+            nodeCount[node] = count;
+            for (int i = lo; i < hi; i++) {
+                int c = perm[i];
+                cubes[filled] = all[c];
+                cubeRank[filled] = c;
+                System.arraycopy(srcBox, c * 6, cubeBox, filled * 6, 6);
+                filled++;
+            }
+            return depth;
+        }
+        nodeCount[node] = 0;
+        int dl = buildNode(all, lo, mid, depth + 1);
+        nodeRight[node] = nodes;
+        int dr = buildNode(all, mid, hi, depth + 1);
+        return Math.max(dl, dr);
+    }
+
+    /**
+     * Partitions perm[lo, hi) along {@code axis} at the cheapest of {@link #SAH_BINS}-1 bin boundaries
+     * (surface area heuristic); falls back to the median when every boundary leaves one side empty. Returns
+     * the split index, or -1 when a leaf is cheaper.
+     */
+    private int sahSplit(int nb, int lo, int hi, int axis, double cmin, double extent) {
+        int[] binCount = new int[SAH_BINS];
+        double[] binBox = new double[SAH_BINS * 6];
+        for (int b = 0; b < SAH_BINS; b++) {
+            for (int a = 0; a < 3; a++) {
+                binBox[b * 6 + a] = Double.MAX_VALUE;
+                binBox[b * 6 + 3 + a] = -Double.MAX_VALUE;
+            }
+        }
+        double scale = SAH_BINS / extent;
+        for (int i = lo; i < hi; i++) {
+            int c = perm[i];
+            int b = bin(centroid[c * 3 + axis], cmin, scale);
+            binCount[b]++;
+            for (int a = 0; a < 3; a++) {
+                binBox[b * 6 + a] = Math.min(binBox[b * 6 + a], srcBox[c * 6 + a]);
+                binBox[b * 6 + 3 + a] = Math.max(binBox[b * 6 + 3 + a], srcBox[c * 6 + 3 + a]);
+            }
+        }
+        // Sweep: cost(split after bin s) = area(left)·nLeft + area(right)·nRight.
+        double[] leftArea = new double[SAH_BINS];
+        int[] leftCount = new int[SAH_BINS];
+        double[] acc = {
+            Double.MAX_VALUE,
+            Double.MAX_VALUE,
+            Double.MAX_VALUE,
+            -Double.MAX_VALUE,
+            -Double.MAX_VALUE,
+            -Double.MAX_VALUE
+        };
+        int cnt = 0;
+        for (int b = 0; b < SAH_BINS - 1; b++) {
+            grow(acc, binBox, b);
+            cnt += binCount[b];
+            leftArea[b] = area(acc);
+            leftCount[b] = cnt;
+        }
+        acc[0] = acc[1] = acc[2] = Double.MAX_VALUE;
+        acc[3] = acc[4] = acc[5] = -Double.MAX_VALUE;
+        cnt = 0;
+        double bestCost = Double.POSITIVE_INFINITY;
+        int bestSplit = -1;
+        for (int b = SAH_BINS - 1; b >= 1; b--) {
+            grow(acc, binBox, b);
+            cnt += binCount[b];
+            int left = leftCount[b - 1];
+            if (left == 0 || cnt == 0) continue;
+            double cost = leftArea[b - 1] * left + area(acc) * cnt;
+            if (cost < bestCost) {
+                bestCost = cost;
+                bestSplit = b;
+            }
+        }
+        int count = hi - lo;
+        if (bestSplit < 0) {
+            // All centroids in one bin (should not happen with extent > 0): median split.
+            sortByCentroid(lo, hi, axis);
+            return lo + count / 2;
+        }
+        // A leaf costs area(node)·count (one cube test per cube), a split one node visit more.
+        double nodeArea = area(nodeBox, nb);
+        if (count <= 16 && bestCost + nodeArea >= nodeArea * count) return -1;
+        int i = lo, j = hi - 1;
+        while (i <= j) {
+            if (bin(centroid[perm[i] * 3 + axis], cmin, scale) < bestSplit) {
+                i++;
+            } else {
+                int tmp = perm[i];
+                perm[i] = perm[j];
+                perm[j--] = tmp;
+            }
+        }
+        return i;
+    }
+
+    private static int bin(double c, double cmin, double scale) {
+        int b = (int) ((c - cmin) * scale);
+        return b < 0 ? 0 : Math.min(b, SAH_BINS - 1);
+    }
+
+    private void sortByCentroid(int lo, int hi, int axis) {
+        Integer[] tmp = new Integer[hi - lo];
+        for (int i = lo; i < hi; i++) tmp[i - lo] = perm[i];
+        java.util.Arrays.sort(tmp, (a, b) -> Double.compare(centroid[a * 3 + axis], centroid[b * 3 + axis]));
+        for (int i = lo; i < hi; i++) perm[i] = tmp[i - lo];
+    }
+
+    private static void grow(double[] acc, double[] boxes, int b) {
+        for (int a = 0; a < 3; a++) {
+            acc[a] = Math.min(acc[a], boxes[b * 6 + a]);
+            acc[3 + a] = Math.max(acc[3 + a], boxes[b * 6 + 3 + a]);
+        }
+    }
+
+    private static double area(double[] box) {
+        return area(box, 0);
+    }
+
+    private static double area(double[] box, int o) {
+        double x = box[o + 3] - box[o], y = box[o + 4] - box[o + 1], z = box[o + 5] - box[o + 2];
+        if (!(x >= 0) || !(y >= 0) || !(z >= 0)) return 0;
+        return 2 * (x * y + y * z + z * x);
     }
 }

src/main/java/eu/mhsl/minecraft/pixelpics/render/raytrace/ElementIntersector.java

@@ -5,7 +5,8 @@ import eu.mhsl.minecraft.pixelpics.assets.model.Element;
 import eu.mhsl.minecraft.pixelpics.assets.model.Face;
 import eu.mhsl.minecraft.pixelpics.assets.model.ResolvedModel;
 import eu.mhsl.minecraft.pixelpics.render.util.ColorUtil;
-import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Comparator;
 import java.util.List;
 import org.bukkit.util.Vector;
 
@@ -22,6 +23,50 @@ public final class ElementIntersector {
 
     private ElementIntersector() {}
 
+    /**
+     * Per-thread scratch of one {@link #intersect} call: the candidate faces as parallel primitive arrays, so
+     * the hot path allocates nothing per element test (the tracer runs rays in parallel, hence per thread).
+     */
+    private static final class Scratch {
+        int n;
+        Element[] element = new Element[16];
+        double[] t = new double[16];
+        Direction[] dir = new Direction[16];
+        double[] s = new double[16];
+        double[] t2 = new double[16];
+        double[] nx = new double[16], ny = new double[16], nz = new double[16];
+        int[] order = new int[16];
+        // Candidate indices to sort; pre-boxed (cached Integers) so the sort allocates nothing.
+        Integer[] idx = new Integer[16];
+        final double[] o = new double[3], d = new double[3], v = new double[3];
+        // Result of the last slab test.
+        double hitT, hitS, hitT2, hitNx, hitNy, hitNz;
+        Direction hitDir;
+
+        // Same ordering as before (depth; coplanar within 1e-4: later elements first) over candidate indices.
+        // Sorted with the same TimSort as the former List<Candidate>.sort, so ties resolve identically.
+        final Comparator<Integer> byDepth = (a, b) -> {
+            if (Math.abs(t[a] - t[b]) > 1e-4) return Double.compare(t[a], t[b]);
+            return Integer.compare(order[b], order[a]);
+        };
+
+        void grow() {
+            int len = element.length * 2;
+            element = Arrays.copyOf(element, len);
+            t = Arrays.copyOf(t, len);
+            dir = Arrays.copyOf(dir, len);
+            s = Arrays.copyOf(s, len);
+            t2 = Arrays.copyOf(t2, len);
+            nx = Arrays.copyOf(nx, len);
+            ny = Arrays.copyOf(ny, len);
+            nz = Arrays.copyOf(nz, len);
+            order = Arrays.copyOf(order, len);
+            idx = new Integer[len];
+        }
+    }
+
+    private static final ThreadLocal<Scratch> SCRATCH = ThreadLocal.withInitial(Scratch::new);
+
     /**
      * @param ox,oy,oz ray origin in block-local coordinates (world origin minus block min corner)
      * @param dx,dy,dz ray direction (need not be normalized)
@@ -40,72 +85,88 @@ public final class ElementIntersector {
             int by,
             int bz,
             double minT) {
-        List<Candidate> candidates = new ArrayList<>(model.elements.size());
-        for (int i = 0; i < model.elements.size(); i++) {
-            Element element = model.elements.get(i);
-            Candidate c;
+        List<Element> elements = model.elements;
+        int count = elements.size();
+        if (count == 0) return null;
+        Scratch sc = SCRATCH.get();
+        sc.n = 0;
+        double[] v = sc.v;
+        for (int i = 0; i < count; i++) {
+            Element element = elements.get(i);
+            boolean hit;
+            Element e;
             if (element.unrotated != null) {
                 // Blockstate rotation: hit the model-space element with the inverse-rotated ray (same t: rigid
                 // transform), so faces and UVs stay bound to the model like in vanilla (uvlock: see Face#uvLock).
-                Element m = element.unrotated;
-                double[] mo = Element.unrotate(ox, oy, oz, element.xSteps, element.ySteps, true);
-                double[] md = Element.unrotate(dx, dy, dz, element.xSteps, element.ySteps, false);
-                c = m.isAxisAligned()
-                        ? intersectAabb(m, mo[0], mo[1], mo[2], md[0], md[1], md[2])
-                        : intersectObb(m, mo[0], mo[1], mo[2], md[0], md[1], md[2]);
-                if (c != null) {
-                    double[] n = Element.rotate(
-                            c.normal[0], c.normal[1], c.normal[2], element.xSteps, element.ySteps, false);
-                    c = new Candidate(c.element, c.t, c.dir, c.s, c.t2, n, 0);
+                e = element.unrotated;
+                unrotate(v, ox, oy, oz, element.xSteps, element.ySteps, true);
+                double mox = v[0], moy = v[1], moz = v[2];
+                unrotate(v, dx, dy, dz, element.xSteps, element.ySteps, false);
+                hit = e.isAxisAligned()
+                        ? slab(sc, mox, moy, moz, v[0], v[1], v[2], e.from, e.to, null)
+                        : intersectObb(sc, e, mox, moy, moz, v[0], v[1], v[2]);
+                if (hit) {
+                    rotate(v, sc.hitNx, sc.hitNy, sc.hitNz, element.xSteps, element.ySteps, false);
+                    sc.hitNx = v[0];
+                    sc.hitNy = v[1];
+                    sc.hitNz = v[2];
                 }
             } else {
-                c = element.isAxisAligned()
-                        ? intersectAabb(element, ox, oy, oz, dx, dy, dz)
-                        : intersectObb(element, ox, oy, oz, dx, dy, dz);
+                e = element;
+                hit = e.isAxisAligned()
+                        ? slab(sc, ox, oy, oz, dx, dy, dz, e.from, e.to, null)
+                        : intersectObb(sc, e, ox, oy, oz, dx, dy, dz);
             }
-            if (c != null && c.t > minT) candidates.add(new Candidate(c.element, c.t, c.dir, c.s, c.t2, c.normal, i));
+            if (!hit || !(sc.hitT > minT)) continue;
+            int k = sc.n;
+            if (k == sc.element.length) sc.grow();
+            sc.element[k] = e;
+            sc.t[k] = sc.hitT;
+            sc.dir[k] = sc.hitDir;
+            sc.s[k] = sc.hitS;
+            sc.t2[k] = sc.hitT2;
+            sc.nx[k] = sc.hitNx;
+            sc.ny[k] = sc.hitNy;
+            sc.nz[k] = sc.hitNz;
+            sc.order[k] = i;
+            sc.n = k + 1;
         }
-        if (candidates.isEmpty()) return null;
+        int n = sc.n;
+        if (n == 0) return null;
+        Integer[] idx = sc.idx;
+        for (int k = 0; k < n; k++) idx[k] = k; // autoboxing hits the Integer cache for k < 128
         // Sort by depth; for coplanar faces (equal t) render later elements first, matching vanilla's
         // draw order so overlays (e.g. the tinted grass side overlay) sit on top of the base face.
-        candidates.sort((a, b) -> {
-            if (Math.abs(a.t - b.t) > 1e-4) return Double.compare(a.t, b.t);
-            return Integer.compare(b.order, a.order);
-        });
+        if (n > 1) Arrays.sort(idx, 0, n, sc.byDepth);
 
-        for (Candidate c : candidates) {
-            Face face = c.element.faces[c.dir.ordinal()];
+        for (int j = 0; j < n; j++) {
+            int k = idx[j];
+            Face face = sc.element[k].faces[sc.dir[k].ordinal()];
             if (face == null) continue;
-            int color = face.sample(c.s, c.t2);
+            int color = face.sample(sc.s[k], sc.t2[k]);
             if (!ColorUtil.opaque(color)) continue;
 
-            Vector world = new Vector(bx + ox + dx * c.t, by + oy + dy * c.t, bz + oz + dz * c.t);
-            Vector normal = new Vector(c.normal[0], c.normal[1], c.normal[2]);
-            return new FaceHit(c.t, world, normal, color, face.tintIndex, false);
+            double t = sc.t[k];
+            Vector world = new Vector(bx + ox + dx * t, by + oy + dy * t, bz + oz + dz * t);
+            Vector normal = new Vector(sc.nx[k], sc.ny[k], sc.nz[k]);
+            return new FaceHit(t, world, normal, color, face.tintIndex, false);
         }
         return null;
     }
 
-    private record Candidate(
-            Element element, double t, Direction dir, double s, double t2, double[] normal, int order) {}
-
-    private static Candidate intersectAabb(
-            Element e, double ox, double oy, double oz, double dx, double dy, double dz) {
-        return slab(e, ox, oy, oz, dx, dy, dz, e.from, e.to, null);
-    }
-
     /** Rotated element: transform the ray into the element's local (unrotated) frame, then slab-test. */
-    private static Candidate intersectObb(Element e, double ox, double oy, double oz, double dx, double dy, double dz) {
-        double[] o = rotate(ox - e.rotOrigin[0], oy - e.rotOrigin[1], oz - e.rotOrigin[2], e.rotAxis, -e.rotAngleRad);
-        o[0] += e.rotOrigin[0];
-        o[1] += e.rotOrigin[1];
-        o[2] += e.rotOrigin[2];
-        double[] d = rotate(dx, dy, dz, e.rotAxis, -e.rotAngleRad);
-        return slab(e, o[0], o[1], o[2], d[0], d[1], d[2], e.from, e.to, e);
+    private static boolean intersectObb(
+            Scratch sc, Element e, double ox, double oy, double oz, double dx, double dy, double dz) {
+        double[] v = sc.v;
+        rotate(v, ox - e.rotOrigin[0], oy - e.rotOrigin[1], oz - e.rotOrigin[2], e.rotAxis, e.invRotCos, e.invRotSin);
+        double lox = v[0] + e.rotOrigin[0], loy = v[1] + e.rotOrigin[1], loz = v[2] + e.rotOrigin[2];
+        rotate(v, dx, dy, dz, e.rotAxis, e.invRotCos, e.invRotSin);
+        return slab(sc, lox, loy, loz, v[0], v[1], v[2], e.from, e.to, e);
     }
 
-    private static Candidate slab(
-            Element e,
+    /** Slab test; on a hit the entry (or, from inside, exit) face is left in the scratch's hit fields. */
+    private static boolean slab(
+            Scratch sc,
             double ox,
             double oy,
             double oz,
@@ -120,11 +181,16 @@ public final class ElementIntersector {
         int axis = -1;
         boolean negFace = false; // entered through the low-coordinate face
 
-        double[] o = {ox, oy, oz};
-        double[] d = {dx, dy, dz};
+        double[] o = sc.o, d = sc.d;
+        o[0] = ox;
+        o[1] = oy;
+        o[2] = oz;
+        d[0] = dx;
+        d[1] = dy;
+        d[2] = dz;
         for (int a = 0; a < 3; a++) {
             if (Math.abs(d[a]) < EPS) {
-                if (o[a] < from[a] - EPS || o[a] > to[a] + EPS) return null;
+                if (o[a] < from[a] - EPS || o[a] > to[a] + EPS) return false;
                 continue;
             }
             double inv = 1.0 / d[a];
@@ -144,15 +210,15 @@ public final class ElementIntersector {
                 negFace = neg;
             }
             if (t2 < tmax) tmax = t2;
-            if (tmin > tmax) return null;
+            if (tmin > tmax) return false;
         }
-        if (axis < 0) return null;
+        if (axis < 0) return false;
 
         double tEntry = tmin;
         if (tEntry < EPS) {
             // origin inside the box (e.g. camera within a block): use exit point instead
             tEntry = tmax;
-            if (tEntry < EPS) return null;
+            if (tEntry < EPS) return false;
         }
 
         double px = o[0] + d[0] * tEntry;
@@ -160,35 +226,51 @@ public final class ElementIntersector {
         double pz = o[2] + d[2] * tEntry;
 
         Direction dir = faceFor(axis, negFace);
-        double[] normal = {dir.nx, dir.ny, dir.nz};
         if (obb != null) {
             // rotate the normal back into block space
-            normal = rotate(normal[0], normal[1], normal[2], obb.rotAxis, obb.rotAngleRad);
+            double[] v = sc.v;
+            rotate(v, dir.nx, dir.ny, dir.nz, obb.rotAxis, obb.rotCos, obb.rotSin);
+            sc.hitNx = v[0];
+            sc.hitNy = v[1];
+            sc.hitNz = v[2];
+        } else {
+            sc.hitNx = dir.nx;
+            sc.hitNy = dir.ny;
+            sc.hitNz = dir.nz;
         }
 
         double fracX = frac(px, from[0], to[0]);
         double fracY = frac(py, from[1], to[1]);
         double fracZ = frac(pz, from[2], to[2]);
 
-        double[] st = faceUv(dir, fracX, fracY, fracZ);
-        double s = st[0], t = st[1];
-        return new Candidate(e, tEntry, dir, s, t, normal, 0);
+        sc.hitT = tEntry;
+        sc.hitDir = dir;
+        sc.hitS = faceS(dir, fracX, fracY, fracZ);
+        sc.hitT2 = faceT(dir, fracY, fracZ);
+        return true;
     }
 
     /**
-     * Face-local (s, t) in 0..1 for a hit at block fractions (fx, fy, fz), matching vanilla's
-     * {@code BlockElement.uvsByFace}: a face's UV rect {@code [u1,v1,u2,v2]} is lerped with these, so north/east
+     * Face-local horizontal coordinate s in 0..1 for a hit at block fractions (fx, fy, fz), matching vanilla's
+     * {@code BlockElement.uvsByFace}: a face's UV rect {@code [u1,v1,u2,v2]} is lerped with (s, t), so north/east
      * faces run their horizontal axis towards -x/-z (left-to-right seen from outside) and down runs v towards
      * north. Shared with the entity intersector (entity box-UV uses the same orientation).
      */
-    public static double[] faceUv(Direction dir, double fx, double fy, double fz) {
+    public static double faceS(Direction dir, double fx, double fy, double fz) {
         return switch (dir) {
-            case UP -> new double[] {fx, fz};
-            case DOWN -> new double[] {fx, 1 - fz};
-            case NORTH -> new double[] {1 - fx, 1 - fy};
-            case SOUTH -> new double[] {fx, 1 - fy};
-            case WEST -> new double[] {fz, 1 - fy};
-            case EAST -> new double[] {1 - fz, 1 - fy};
+            case UP, DOWN, SOUTH -> fx;
+            case NORTH -> 1 - fx;
+            case WEST -> fz;
+            case EAST -> 1 - fz;
+        };
+    }
+
+    /** Face-local vertical coordinate t in 0..1 (see {@link #faceS}); t does not depend on fx. */
+    public static double faceT(Direction dir, double fy, double fz) {
+        return switch (dir) {
+            case UP -> fz;
+            case DOWN -> 1 - fz;
+            default -> 1 - fy;
         };
     }
 
@@ -207,14 +289,49 @@ public final class ElementIntersector {
         return f < 0 ? 0 : Math.min(f, 1);
     }
 
-    /** Rotate (x,y,z) around the given axis (0=x,1=y,2=z) by angle radians. */
-    private static double[] rotate(double x, double y, double z, int axis, double angle) {
-        double cos = Math.cos(angle);
-        double sin = Math.sin(angle);
-        return switch (axis) {
-            case 0 -> new double[] {x, y * cos - z * sin, y * sin + z * cos};
-            case 1 -> new double[] {x * cos + z * sin, y, -x * sin + z * cos};
-            default -> new double[] {x * cos - y * sin, x * sin + y * cos, z};
-        };
+    /** Rotate (x,y,z) around the given axis (0=x,1=y,2=z) by the angle with the given cos/sin, into out. */
+    private static void rotate(double[] out, double x, double y, double z, int axis, double cos, double sin) {
+        switch (axis) {
+            case 0 -> {
+                out[0] = x;
+                out[1] = y * cos - z * sin;
+                out[2] = y * sin + z * cos;
+            }
+            case 1 -> {
+                out[0] = x * cos + z * sin;
+                out[1] = y;
+                out[2] = -x * sin + z * cos;
+            }
+            default -> {
+                out[0] = x * cos - y * sin;
+                out[1] = x * sin + y * cos;
+                out[2] = z;
+            }
+        }
+    }
+
+    /** {@link Element#rotate} into out (same arithmetic, no allocation). */
+    private static void rotate(double[] out, double x, double y, double z, int xSteps, int ySteps, boolean point) {
+        double c = point ? 0.5 : 0;
+        double dx = x - c, dy = y - c, dz = z - c;
+        for (int i = 0; i < xSteps; i++) {
+            double ny = dz, nz = -dy; // (dy, dz) -> (dz, -dy)
+            dy = ny;
+            dz = nz;
+        }
+        for (int i = 0; i < ySteps; i++) {
+            double nx = -dz, nz = dx; // (dx, dz) -> (-dz, dx)
+            dx = nx;
+            dz = nz;
+        }
+        out[0] = dx + c;
+        out[1] = dy + c;
+        out[2] = dz + c;
+    }
+
+    /** {@link Element#unrotate} into out (same arithmetic, no allocation). */
+    private static void unrotate(double[] out, double x, double y, double z, int xSteps, int ySteps, boolean point) {
+        rotate(out, x, y, z, 0, (4 - ySteps) % 4, point);
+        rotate(out, out[0], out[1], out[2], (4 - xSteps) % 4, 0, point);
     }
 }