The overlap check of “Fast Parallel Surface and Solid Voxelisation on GPUs” by Schwarz and Seidel in C# in Unity visualised on a single triangle

TriangleBB.cs

using System;
using System.Collections.Generic;
using UnityEngine;

[Flags]
public enum BitValues
{
	_0_ =  0,
	_1_ = 1 << 0,
	_2_ = 1 << 1,
	_4_ = 1 << 2,
	_8_ = 1 << 3,
	_16_ = 1 << 4,
	_32_ = 1 << 5,
	_64_ = 1 << 6,
	_128_ = 1 << 7,
	_256_ = 1 << 8,
	_512_ = 1 << 9,
	_1024_ = 1 << 10,
	_2048_ = 1 << 11,
	_4096_ = 1 << 12,
	_8192_ = 1 << 13,
	_16384_ = 1 << 14,
	_32768_ = 1 << 15,
	_65536_ = 1 << 16,
	_131072_ = 1 << 17,
	_262144_ = 1 << 18,
	_524288_ = 1 << 19,
	_1048576_ = 1 << 20,
	_2097152_ = 1 << 21,
	_4194304_ = 1 << 22,
	_8388608_ = 1 << 23,
	_16777216_ = 1 << 24,
	_33554432_ = 1 << 25,
	_67108864_ = 1 << 26,
	_134217728_= 1 << 27,
	_268435456_ = 1 << 28,
	_536870912_ = 1 << 29,
	_1073741824_ = 1 << 30,
	_2147483648_ = 1 << 31
}


public class TriangleBB : MonoBehaviour
{
	[Header("Triangle")]
	[SerializeField] private MeshFilter triangleFilter;
	[SerializeField] private Transform triangleTransform;

	[Header("Do not change manually!")]
	public Vector3 v0;
	public Vector3 v1;
	public Vector3 v2;
	public Vector3 n;
	public Vector3 e0 => v1 - v0;
	public Vector3 ce0 => -e0 / 2 + v1;
	public Vector3 e1 => v2 - v1;
	public Vector3 ce1 => -e1 / 2 + v2;
	public Vector3 e2 => v0 - v2;
	public Vector3 ce2 => -e2 / 2 + v0;
	private Vector3 center => 1 / 3f * ( v0 + v1 + v2 );

	[Header("Grid")]
	[SerializeField] private Vector3Int grid = new (2, 2, 2);
	[SerializeField] private uint cellSize = 1;   // Is deltaP (dp)

	[Header("Debug Overlapped")]
	[SerializeField] private int[] _voxels;
	[SerializeField] private BitValues _voxels0ToBits;
	[SerializeField] private List<uint> _indices;

	private void OnDrawGizmos()
	{
		_voxels = new int[Math.Max(grid.x * grid.y * grid.z / 32, 1)];
		_voxels0ToBits = BitValues._0_;
		_indices = new List<uint>();

		// Triangle data from mesh
		n = Vector3.Cross(e0, e1).normalized;
		Matrix4x4 localToWorld = triangleTransform.localToWorldMatrix;
		var mesh = triangleFilter.sharedMesh;
		v0 = localToWorld.MultiplyPoint3x4(mesh.vertices[0]);
		v1 = localToWorld.MultiplyPoint3x4(mesh.vertices[1]);
		v2 = localToWorld.MultiplyPoint3x4(mesh.vertices[2]);
		// Gizmos.color = Color.green;
		// Gizmos.DrawSphere(v0, 0.05f);
		// Gizmos.DrawSphere(v1, 0.05f);
		// Gizmos.DrawSphere(v2, 0.05f);

		// Triangle world bounding box
		var min = Vector3.Min(v0, Vector3.Min(v1, v2));
		var max = Vector3.Max(v0, Vector3.Max(v1, v2));
		Gizmos.color = Color.yellow;
		Gizmos.DrawSphere(min, 0.05f);
		Gizmos.DrawSphere(max, 0.05f);
		Gizmos.DrawWireCube(min + ( max - min ) / 2, max - min);
	
		// Grid
		var dp = new Vector3(cellSize, cellSize, cellSize);
		Gizmos.color = Color.blue;
		for (int x = 0; x < grid.x; x++)
		for (int y = 0; y < grid.y; y++)
		for (int z = 0; z < grid.z; z++)
		{
			var voxelCenter = new Vector3(x * cellSize + cellSize / 2f, y * cellSize + cellSize / 2f, z * cellSize + cellSize / 2f);
			Gizmos.DrawWireCube(voxelCenter, dp);
		}
	
		// Triangle bounding box in voxel space (all voxels overlapped by triangle world bounding box)
		Vector3Int voxelBBMin = Vector3Int.FloorToInt(min / cellSize) * (int)cellSize;
		Vector3Int voxelBBMax = Vector3Int.CeilToInt(max / cellSize) * (int)cellSize;
		Gizmos.color = Color.red;
		Gizmos.DrawSphere(voxelBBMin, 0.05f);
		Gizmos.DrawSphere(voxelBBMax, 0.05f);
		// Color the voxels inside the triangle's world bounding box that pass the equations' tests red, otherwise black.
		for (int z = voxelBBMin.z; z < voxelBBMax.z; z += (int)cellSize)
		for (int y = voxelBBMin.y; y < voxelBBMax.y; y += (int)cellSize)
		for (int x = voxelBBMin.x; x < voxelBBMax.x; x += (int)cellSize)
		{
			var voxelP = new Vector3(x, y, z);
			var voxelCenter = new Vector3(x + cellSize / 2f, y + cellSize / 2f, z + cellSize / 2f);
			bool overlapped = CalculateEquation1(voxelP, dp) && CalculateEquation2And3(voxelP, dp);
			Gizmos.color = overlapped ? Color.red : Color.black;
			Gizmos.DrawWireCube(voxelCenter, dp);
			if (overlapped)
			{
				int index = x + grid.x * (y + grid.y * z);
				int shift = index % 32;
				// I think index needs to be divided by 32 because each one contains 32 voxels.
				_voxels[index / 32] |= 1 << shift;
				_voxels0ToBits |= (BitValues)(1 << shift);
				_indices.Add((uint)index);
			}
		}
	}

	public bool CalculateEquation1(Vector3 p, Vector3 dp)
	{
		var c = GetCriticalPoint(dp);
		var np = Vector3.Dot(n, p);
		var d1 = Vector3.Dot(n, c - v0);
		var d2 = Vector3.Dot(n, ( dp - c ) - v0);
		var result = (np + d1) * (np + d2);
		var trianglePlaneOverlapsBox = result <= 0;
		return trianglePlaneOverlapsBox;
	}

	// Returns true if condition of equation 3 is met.
	public bool CalculateEquation2And3(Vector3 p, Vector3 dp)
	{
		// Equation 2
		var neXY = GetEdgeNormalsXY();
		var neYZ = GetEdgeNormalsYZ();
		var neZX = GetEdgeNormalsZX();
		var dXY = GetEdgeFunctionValuesXY(dp, neXY);
		var dYZ = GetEdgeFunctionValuesYZ(dp, neYZ);
		var dZX = GetEdgeFunctionValuesZX(dp, neZX);
	
		// Equation 3
	
		// XY
		var pXY = new Vector2(p.x, p.y);
		if (Vector2.Dot(neXY.e0, pXY) + dXY.e0 < 0) return false;
		if (Vector2.Dot(neXY.e1, pXY) + dXY.e1 < 0) return false;
		if (Vector2.Dot(neXY.e2, pXY) + dXY.e2 < 0) return false;
	
		// YZ
		var pYZ = new Vector2(p.y, p.z);
		if (Vector2.Dot(neYZ.e0, pYZ) + dYZ.e0 < 0) return false;
		if (Vector2.Dot(neYZ.e1, pYZ) + dYZ.e1 < 0) return false;
		if (Vector2.Dot(neYZ.e2, pYZ) + dYZ.e2 < 0) return false;
	
		// ZX
		var pZX = new Vector2(p.z, p.x);
		if (Vector2.Dot(neZX.e0, pZX) + dZX.e0 < 0) return false;
		if (Vector2.Dot(neZX.e1, pZX) + dZX.e1 < 0) return false;
		if (Vector2.Dot(neZX.e2, pZX) + dZX.e2 < 0) return false;

		return true;
	}

	public Vector3 GetCriticalPoint(Vector3 dp)
	{
		return new Vector3(n.x > 0 ? dp.x : 0, n.y > 0 ? dp.y : 0, n.z > 0 ? dp.z : 0);
	}

	public (Vector2 e0, Vector2 e1, Vector2 e2) GetEdgeNormalsXY()
	{
		var zMultiplier = n.z < 0.0f ? -1 : 1;
		return
		(
			new Vector2(-e0.y, e0.x) * zMultiplier,
			new Vector2(-e1.y, e1.x) * zMultiplier,
			new Vector2(-e2.y, e2.x) * zMultiplier
		);
	}

	public (Vector2 e0, Vector2 e1, Vector2 e2) GetEdgeNormalsYZ()
	{
		var xMultiplier = n.x < 0.0f ? -1 : 1;
		return
		(
			new Vector2(-e0.z, e0.y) * xMultiplier,
			new Vector2(-e1.z, e1.y) * xMultiplier,
			new Vector2(-e2.z, e2.y) * xMultiplier
		);
	}

	public (Vector2 e0, Vector2 e1, Vector2 e2) GetEdgeNormalsZX()
	{
		var yMultiplier = n.y < 0.0f ? -1 : 1;
		return
		(
			new Vector2(-e0.x, e0.z) * yMultiplier,
			new Vector2(-e1.x, e1.z) * yMultiplier,
			new Vector2(-e2.x, e2.z) * yMultiplier
		);
	}

	public (float e0, float e1, float e2) GetEdgeFunctionValuesXY(Vector3 dp, (Vector2 e0, Vector2 e1, Vector2 e2) neXY)
	{
		return
		(
			-Vector2.Dot(neXY.e0, new Vector2(v0.x, v0.y)) + Mathf.Max(0.0f, dp.x * neXY.e0.x) + Mathf.Max(0.0f, dp.y * neXY.e0.y), // e0
			-Vector2.Dot(neXY.e1, new Vector2(v1.x, v1.y)) + Mathf.Max(0.0f, dp.x * neXY.e1.x) + Mathf.Max(0.0f, dp.y * neXY.e1.y), // e1
			-Vector2.Dot(neXY.e2, new Vector2(v2.x, v2.y)) + Mathf.Max(0.0f, dp.x * neXY.e2.x) + Mathf.Max(0.0f, dp.y * neXY.e2.y) // e2
		);
	}


	public (float e0, float e1, float e2) GetEdgeFunctionValuesYZ(Vector3 dp, (Vector2 e0, Vector2 e1, Vector2 e2) neYZ)
	{
		return
		(
			-Vector2.Dot(neYZ.e0, new Vector2(v0.y, v0.z)) + Mathf.Max(0.0f, dp.y * neYZ.e0.x) + Mathf.Max(0.0f, dp.z * neYZ.e0.y), // e0
			-Vector2.Dot(neYZ.e1, new Vector2(v1.y, v1.z)) + Mathf.Max(0.0f, dp.y * neYZ.e1.x) + Mathf.Max(0.0f, dp.z * neYZ.e1.y), // e1
			-Vector2.Dot(neYZ.e2, new Vector2(v2.y, v2.z)) + Mathf.Max(0.0f, dp.y * neYZ.e2.x) + Mathf.Max(0.0f, dp.z * neYZ.e2.y) // e2
		);
	}

	public (float e0, float e1, float e2) GetEdgeFunctionValuesZX(Vector3 dp, (Vector2 e0, Vector2 e1, Vector2 e2) neZX)
	{
		return
		(
			-Vector2.Dot(neZX.e0, new Vector2(v0.z, v0.x)) + Mathf.Max(0.0f, dp.z * neZX.e0.x) + Mathf.Max(0.0f, dp.x * neZX.e0.y), // e0
			-Vector2.Dot(neZX.e1, new Vector2(v1.z, v1.x)) + Mathf.Max(0.0f, dp.z * neZX.e1.x) + Mathf.Max(0.0f, dp.x * neZX.e1.y), // e1
			-Vector2.Dot(neZX.e2, new Vector2(v2.z, v2.x)) + Mathf.Max(0.0f, dp.z * neZX.e2.x) + Mathf.Max(0.0f, dp.x * neZX.e2.y) // e2
		);
	}
}