Home Reference Source

src/octree/world/WorldOctreeRaycaster.js

import { Box3, Line3, Ray, Vector3 } from "three";
import { layout } from "sparse-octree";
import { WorldOctantWrapper } from "./WorldOctantWrapper";

/**
 * A vector.
 *
 * @type {Vector3}
 * @private
 */

const v = new Vector3();

/**
 * A line.
 *
 * @type {Line3}
 * @private
 */

const l = new Line3();

/**
 * A box.
 *
 * @type {Box3}
 * @private
 */

const b = new Box3();

/**
 * A box.
 *
 * @type {Box3}
 * @private
 */

const d = new Box3();

/**
 * A ray.
 *
 * @type {Ray}
 * @private
 */

const r = new Ray();

/**
 * A lookup-table containing octant ids. Used to determine the exit plane from
 * an octant.
 *
 * @type {Uint8Array[]}
 * @private
 */

const octantTable = [

	new Uint8Array([4, 2, 1]),
	new Uint8Array([5, 3, 8]),
	new Uint8Array([6, 8, 3]),
	new Uint8Array([7, 8, 8]),
	new Uint8Array([8, 6, 5]),
	new Uint8Array([8, 7, 8]),
	new Uint8Array([8, 8, 7]),
	new Uint8Array([8, 8, 8])

];

/**
 * A byte that stores raycasting flags.
 *
 * @type {Number}
 * @private
 */

let flags = 0;

/**
 * Finds the entry plane of the first octant that a ray travels through.
 *
 * Determining the first octant requires knowing which of the t0s is the
 * largest. The tms of the other axes must also be compared against that
 * largest t0.
 *
 * @private
 * @param {Number} tx0 - Ray projection parameter.
 * @param {Number} ty0 - Ray projection parameter.
 * @param {Number} tz0 - Ray projection parameter.
 * @param {Number} txm - Ray projection parameter mean.
 * @param {Number} tym - Ray projection parameter mean.
 * @param {Number} tzm - Ray projection parameter mean.
 * @return {Number} The index of the first octant that the ray travels through.
 */

function findEntryOctant(tx0, ty0, tz0, txm, tym, tzm) {

	let entry = 0;

	// Find the entry plane.
	if(tx0 > ty0 && tx0 > tz0) {

		// YZ-plane.
		if(tym < tx0) {

			entry |= 2;

		}

		if(tzm < tx0) {

			entry |= 1;

		}

	} else if(ty0 > tz0) {

		// XZ-plane.
		if(txm < ty0) {

			entry |= 4;

		}

		if(tzm < ty0) {

			entry |= 1;

		}

	} else {

		// XY-plane.
		if(txm < tz0) {

			entry |= 4;

		}

		if(tym < tz0) {

			entry |= 2;

		}

	}

	return entry;

}

/**
 * Finds the next octant that intersects with the ray based on the exit plane of
 * the current one.
 *
 * @private
 * @param {Number} currentOctant - The index of the current octant.
 * @param {Number} tx1 - Ray projection parameter.
 * @param {Number} ty1 - Ray projection parameter.
 * @param {Number} tz1 - Ray projection parameter.
 * @return {Number} The index of the next octant that the ray travels through.
 */

function findNextOctant(currentOctant, tx1, ty1, tz1) {

	let min;
	let exit = 0;

	// Find the exit plane.
	if(tx1 < ty1) {

		min = tx1;
		exit = 0; // YZ-plane.

	} else {

		min = ty1;
		exit = 1; // XZ-plane.

	}

	if(tz1 < min) {

		exit = 2; // XY-plane.

	}

	return octantTable[currentOctant][exit];

}

/**
 * Recursively traverses the given octant to find (pseudo) leaf octants that
 * intersect with the given ray.
 *
 * @private
 * @param {WorldOctree} world - The world octree.
 * @param {WorldOctant} octant - The current octant.
 * @param {Number} keyX - The X-coordinate of the current octant key.
 * @param {Number} keyY - The Y-coordinate of the current octant key.
 * @param {Number} keyZ - The Z-coordinate of the current octant key.
 * @param {Number} lod - The current LOD.
 * @param {Number} tx0 - Ray projection parameter. Initial tx0 = (minX - rayOriginX) / rayDirectionX.
 * @param {Number} ty0 - Ray projection parameter. Initial ty0 = (minY - rayOriginY) / rayDirectionY.
 * @param {Number} tz0 - Ray projection parameter. Initial tz0 = (minZ - rayOriginZ) / rayDirectionZ.
 * @param {Number} tx1 - Ray projection parameter. Initial tx1 = (maxX - rayOriginX) / rayDirectionX.
 * @param {Number} ty1 - Ray projection parameter. Initial ty1 = (maxY - rayOriginY) / rayDirectionY.
 * @param {Number} tz1 - Ray projection parameter. Initial tz1 = (maxZ - rayOriginZ) / rayDirectionZ.
 * @param {WorldOctant[]} intersects - An array to be filled with the intersecting octants.
 */

function raycastOctant(world, octant, keyX, keyY, keyZ, lod, tx0, ty0, tz0, tx1, ty1, tz1, intersects) {

	let keyDesign, cellSize;
	let octantWrapper, grid;
	let children, offset;

	let currentOctant;
	let txm, tym, tzm;

	let i;

	if(tx1 >= 0.0 && ty1 >= 0.0 && tz1 >= 0.0) {

		keyDesign = world.getKeyDesign();

		if(lod === 0 || octant.isosurface !== null) {

			cellSize = world.getCellSize(lod);
			octantWrapper = new WorldOctantWrapper(octant);
			octantWrapper.id.set(lod, keyDesign.packKey(v.set(keyX, keyY, keyZ)));
			octantWrapper.min.copy(v).multiplyScalar(cellSize).add(world.min);
			octantWrapper.max.copy(octantWrapper.min).addScalar(cellSize);

			intersects.push(octantWrapper);

		} else if(octant.children > 0) {

			// Look at the next lower LOD.
			grid = world.getGrid(--lod);
			children = octant.children;

			// Translate the key coordinates to the next lower LOD.
			keyX <<= 1; keyY <<= 1; keyZ <<= 1;

			// Compute means.
			txm = 0.5 * (tx0 + tx1);
			tym = 0.5 * (ty0 + ty1);
			tzm = 0.5 * (tz0 + tz1);

			currentOctant = findEntryOctant(tx0, ty0, tz0, txm, tym, tzm);

			do {

				i = flags ^ currentOctant;

				switch(currentOctant) {

					case 0: {

						if((children & (1 << i)) !== 0) {

							offset = layout[i];
							v.set(keyX + offset[0], keyY + offset[1], keyZ + offset[2]);
							raycastOctant(world, grid.get(keyDesign.packKey(v)), v.x, v.y, v.z, lod, tx0, ty0, tz0, txm, tym, tzm, intersects);

						}

						currentOctant = findNextOctant(currentOctant, txm, tym, tzm);
						break;

					}

					case 1: {

						if((children & (1 << i)) !== 0) {

							offset = layout[i];
							v.set(keyX + offset[0], keyY + offset[1], keyZ + offset[2]);
							raycastOctant(world, grid.get(keyDesign.packKey(v)), v.x, v.y, v.z, lod, tx0, ty0, tzm, txm, tym, tz1, intersects);

						}

						currentOctant = findNextOctant(currentOctant, txm, tym, tz1);
						break;

					}

					case 2: {

						if((children & (1 << i)) !== 0) {

							offset = layout[i];
							v.set(keyX + offset[0], keyY + offset[1], keyZ + offset[2]);
							raycastOctant(world, grid.get(keyDesign.packKey(v)), v.x, v.y, v.z, lod, tx0, tym, tz0, txm, ty1, tzm, intersects);

						}

						currentOctant = findNextOctant(currentOctant, txm, ty1, tzm);
						break;

					}

					case 3: {

						if((children & (1 << i)) !== 0) {

							offset = layout[i];
							v.set(keyX + offset[0], keyY + offset[1], keyZ + offset[2]);
							raycastOctant(world, grid.get(keyDesign.packKey(v)), v.x, v.y, v.z, lod, tx0, tym, tzm, txm, ty1, tz1, intersects);

						}

						currentOctant = findNextOctant(currentOctant, txm, ty1, tz1);
						break;

					}

					case 4: {

						if((children & (1 << i)) !== 0) {

							offset = layout[i];
							v.set(keyX + offset[0], keyY + offset[1], keyZ + offset[2]);
							raycastOctant(world, grid.get(keyDesign.packKey(v)), v.x, v.y, v.z, lod, txm, ty0, tz0, tx1, tym, tzm, intersects);

						}

						currentOctant = findNextOctant(currentOctant, tx1, tym, tzm);
						break;

					}

					case 5: {

						if((children & (1 << i)) !== 0) {

							offset = layout[i];
							v.set(keyX + offset[0], keyY + offset[1], keyZ + offset[2]);
							raycastOctant(world, grid.get(keyDesign.packKey(v)), v.x, v.y, v.z, lod, txm, ty0, tzm, tx1, tym, tz1, intersects);

						}

						currentOctant = findNextOctant(currentOctant, tx1, tym, tz1);
						break;

					}

					case 6: {

						if((children & (1 << i)) !== 0) {

							offset = layout[i];
							v.set(keyX + offset[0], keyY + offset[1], keyZ + offset[2]);
							raycastOctant(world, grid.get(keyDesign.packKey(v)), v.x, v.y, v.z, lod, txm, tym, tz0, tx1, ty1, tzm, intersects);

						}

						currentOctant = findNextOctant(currentOctant, tx1, ty1, tzm);
						break;

					}

					case 7: {

						if((children & (1 << i)) !== 0) {

							offset = layout[i];
							v.set(keyX + offset[0], keyY + offset[1], keyZ + offset[2]);
							raycastOctant(world, grid.get(keyDesign.packKey(v)), v.x, v.y, v.z, lod, txm, tym, tzm, tx1, ty1, tz1, intersects);

						}

						// Far top right octant. No other octants can be reached from here.
						currentOctant = 8;
						break;

					}

				}

			} while(currentOctant < 8);

		}

	}

}

/**
 * Finds (pseudo) leaf octants in the given subtree that intersect with the
 * given ray.
 *
 * @private
 * @param {WorldOctree} world - The world octree.
 * @param {WorldOctantWrapper} subtree - A world octant, enriched with positional information.
 * @param {Vector3} keyCoordinates - The key coordinates of the octant.
 * @param {Ray} ray - A ray.
 * @param {WorldOctant[]} intersects - The intersecting octants. Sorted by distance, closest first
 */

function intersectSubtree(world, subtree, keyCoordinates, ray, intersects) {

	// Translate the octant extents to the scene origin.
	const min = b.min.set(0, 0, 0);
	const max = b.max.subVectors(subtree.max, subtree.min);

	const dimensions = subtree.getDimensions(d.min);
	const halfDimensions = d.max.copy(dimensions).multiplyScalar(0.5);

	const origin = r.origin.copy(ray.origin);
	const direction = r.direction.copy(ray.direction);

	let invDirX, invDirY, invDirZ;
	let tx0, tx1, ty0, ty1, tz0, tz1;

	// Translate the ray to the center of the octant.
	origin.sub(subtree.getCenter(v)).add(halfDimensions);

	// Reset all flags.
	flags = 0;

	// Handle rays with negative directions.
	if(direction.x < 0.0) {

		origin.x = dimensions.x - origin.x;
		direction.x = -direction.x;
		flags |= 4;

	}

	if(direction.y < 0.0) {

		origin.y = dimensions.y - origin.y;
		direction.y = -direction.y;
		flags |= 2;

	}

	if(direction.z < 0.0) {

		origin.z = dimensions.z - origin.z;
		direction.z = -direction.z;
		flags |= 1;

	}

	// Improve IEEE double stability.
	invDirX = 1.0 / direction.x;
	invDirY = 1.0 / direction.y;
	invDirZ = 1.0 / direction.z;

	// Project the ray to the octant's boundaries.
	tx0 = (min.x - origin.x) * invDirX;
	tx1 = (max.x - origin.x) * invDirX;
	ty0 = (min.y - origin.y) * invDirY;
	ty1 = (max.y - origin.y) * invDirY;
	tz0 = (min.z - origin.z) * invDirZ;
	tz1 = (max.z - origin.z) * invDirZ;

	// Find the intersecting children.
	raycastOctant(
		world, subtree.octant,
		keyCoordinates.x, keyCoordinates.y, keyCoordinates.z, world.getDepth(),
		tx0, ty0, tz0, tx1, ty1, tz1,
		intersects
	);

}

/**
 * A world octree raycaster.
 *
 * This raycaster is a specialised hybrid that uses a voxel traversal algorithm
 * to iterate over the octants of the highest LOD grid and an octree traversal
 * algorithm to raycast the identified subtrees.
 *
 * The voxel traversal implementation is a 3D supercover variant of the Digital
 * Differential Analyzer (DDA) line algorithm and is similar to the Bresenham
 * algorithm. The octree traversal algorithm relies on octant child existence
 * information to skip empty space and to avoid hashmap lookup misses.
 *
 * References:
 *
 *  "Voxel Traversal along a 3D Line"
 *  by D. Cohen (1994)
 *
 *  "An Efficient Parametric Algorithm for Octree Traversal"
 *  by J. Revelles et al. (2000)
 */

export class WorldOctreeRaycaster {

	/**
	 * Finds (pseudo) leaf octants that intersect with the given ray.
	 *
	 * @param {WorldOctree} world - A world octree.
	 * @param {Ray} ray - A ray.
	 * @param {Array} [intersects] - An optional target list to be filled with the intersecting octants.
	 * @return {WorldOctant[]} The intersecting octants. Sorted by distance, closest first.
	 */

	static intersectWorldOctree(world, ray, intersects = []) {

		const lod = world.getDepth();
		const grid = world.getGrid(lod);
		const cellSize = world.getCellSize(lod);
		const keyDesign = world.getKeyDesign();
		const octantWrapper = new WorldOctantWrapper();

		const keyCoordinates0 = l.start;
		const keyCoordinates1 = l.end;

		// Find the point at which the ray enters the world grid.
		const a = !world.containsPoint(r.copy(ray).origin) ?
			r.intersectBox(world, r.origin) :
			r.origin;

		let key, octant;
		let t, b, n;

		let dx, dy, dz;
		let ax, ay, az, bx, by, bz;
		let sx, sy, sz, exy, exz, ezy;

		octantWrapper.id.lod = lod;

		// Check if the ray hits the world octree.
		if(a !== null) {

			// Phase 1: Initialisation.

			// Find the ending point.
			t = cellSize << 1;
			b = r.at(t, v);

			// Calculate the starting and ending cell coordinates.
			world.calculateKeyCoordinates(a, lod, keyCoordinates0);
			world.calculateKeyCoordinates(b, lod, keyCoordinates1);

			// Calculate the key coordinate vector from start to end.
			dx = keyCoordinates1.x - keyCoordinates0.x;
			dy = keyCoordinates1.y - keyCoordinates0.y;
			dz = keyCoordinates1.z - keyCoordinates0.z;

			// Prepare step sizes and project the line onto the XY-, XZ- and ZY-plane.
			sx = Math.sign(dx); sy = Math.sign(dy); sz = Math.sign(dz);
			ax = Math.abs(dx); ay = Math.abs(dy); az = Math.abs(dz);
			bx = 2 * ax; by = 2 * ay; bz = 2 * az;
			exy = ay - ax; exz = az - ax; ezy = ay - az;

			// Phase 2: Incremental Traversal.
			for(n = ax + ay + az; n > 0; --n) {

				key = keyDesign.packKey(keyCoordinates0);

				// Check if this cell is populated.
				if(grid.has(key)) {

					octant = grid.get(key);

					// Setup a pseudo octree.
					octantWrapper.id.key = key;
					octantWrapper.octant = octant;
					octantWrapper.min.copy(keyCoordinates0);
					octantWrapper.min.multiplyScalar(cellSize);
					octantWrapper.min.add(world.min);
					octantWrapper.max.copy(octantWrapper.min).addScalar(cellSize);

					if(octant.isosurface === null) {

						// Raycast the subtree and collect intersecting children.
						intersectSubtree(world, octantWrapper, keyCoordinates0, ray, intersects);

					} else {

						// The octant contains a mesh. No need to look deeper.
						intersects.push(octantWrapper.clone());

					}

				}

				if(exy < 0) {

					if(exz < 0) {

						keyCoordinates0.x += sx;
						exy += by; exz += bz;

					} else {

						keyCoordinates0.z += sz;
						exz -= bx; ezy += by;

					}

				} else if(ezy < 0) {

					keyCoordinates0.z += sz;
					exz -= bx; ezy += by;

				} else {

					keyCoordinates0.y += sy;
					exy -= bx; ezy -= bz;

				}

			}

		}

		return intersects;

	}

}