Initial 3D blocks support

libs/cpml/LICENSE.md

@@ -0,0 +1,60 @@
+# Licenses
+
+CPML is Copyright (c) 2016 Colby Klein <shakesoda@gmail.com>.
+
+CPML is Copyright (c) 2016 Landon Manning <lmanning17@gmail.com>.
+
+Code in vec3.lua is derived from hump.vector. (c) 2010-2013 Matthias Richter. MIT.
+
+Portions of mat4.lua are from LuaMatrix, (c) 2010 Michael Lutz. MIT.
+
+Code in simplex.lua is (c) 2011 Stefan Gustavson. MIT.
+
+Code in bound2.lua and bound3.lua are (c) 2018 Andi McClure. MIT.
+
+Code in quat.lua is from Andrew Stacey and covered under the CC0 license.
+
+Code in octree.lua is derived from UnityOctree. (c) 2014 Nition. BSD-2-Clause.
+
+# The MIT License (MIT)
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+# The BSD License (BSD-2-Clause)
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this
+  list of conditions and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright notice,
+  this list of conditions and the following disclaimer in the documentation
+  and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

libs/cpml/_private_precond.lua

@@ -0,0 +1,17 @@
+-- Preconditions for cpml functions.
+local precond = {}
+
+
+function precond.typeof(t, expected, msg)
+	if type(t) ~= expected then
+		error(("%s: %s (<%s> expected)"):format(msg, type(t), expected), 3)
+	end
+end
+
+function precond.assert(cond, msg, ...)
+	if not cond then
+		error(msg:format(...), 3)
+	end
+end
+
+return precond

libs/cpml/_private_utils.lua

@@ -0,0 +1,16 @@
+-- Functions exported by utils.lua but needed by vec2 or vec3 (which utils.lua requires)
+
+local private = {}
+local floor   = math.floor
+local ceil    = math.ceil
+
+function private.round(value, precision)
+	if precision then return private.round(value / precision) * precision end
+	return value >= 0 and floor(value+0.5) or ceil(value-0.5)
+end
+
+function private.is_nan(a)
+	return a ~= a
+end
+
+return private

libs/cpml/bound2.lua

@@ -0,0 +1,199 @@
+--- A 2 component bounding box.
+-- @module bound2
+
+local modules = (...):gsub('%.[^%.]+$', '') .. "."
+local vec2    = require(modules .. "vec2")
+
+local bound2    = {}
+local bound2_mt = {}
+
+-- Private constructor.
+local function new(min, max)
+	return setmetatable({
+		min=min, -- min: vec2, minimum value for each component 
+		max=max, -- max: vec2, maximum value for each component 
+	}, bound2_mt)
+end
+
+-- Do the check to see if JIT is enabled. If so use the optimized FFI structs.
+local status, ffi
+if type(jit) == "table" and jit.status() then
+	status, ffi = pcall(require, "ffi")
+	if status then
+		ffi.cdef "typedef struct { cpml_vec2 min, max; } cpml_bound2;"
+		new = ffi.typeof("cpml_bound2")
+	end
+end
+
+bound2.zero = new(vec2.zero, vec2.zero)
+
+--- The public constructor.
+-- @param min Can be of two types: </br>
+-- vec2 min, minimum value for each component
+-- nil Create bound at single point 0,0
+-- @tparam vec2 max, maximum value for each component
+-- @treturn bound2 out
+function bound2.new(min, max)
+	if min and max then
+		return new(min:clone(), max:clone())
+	elseif min or max then
+		error("Unexpected nil argument to bound2.new")
+	else
+		return new(vec2.zero, vec2.zero)
+	end
+end
+
+--- Clone a bound.
+-- @tparam bound2 a bound to be cloned
+-- @treturn bound2 out
+function bound2.clone(a)
+	return new(a.min, a.max)
+end
+
+--- Construct a bound covering one or two points 
+-- @tparam vec2 a Any vector
+-- @tparam vec2 b Any second vector (optional)
+-- @treturn vec2 Minimum bound containing the given points
+function bound2.at(a, b) -- "bounded by". b may be nil
+	if b then
+		return bound2.new(a,b):check()
+	else
+		return bound2.zero:with_center(a)
+	end
+end
+
+--- Extend bound to include point
+-- @tparam bound2 a bound
+-- @tparam vec2 point to include
+-- @treturn bound2 Bound covering current min, current max and new point
+function bound2.extend(a, center)
+	return bound2.new(a.min:component_min(center), a.max:component_max(center))
+end
+
+--- Extend bound to entirety of other bound
+-- @tparam bound2 a bound
+-- @tparam bound2 bound to cover
+-- @treturn bound2 Bound covering current min and max of each bound in the pair
+function bound2.extend_bound(a, b)
+	return a:extend(b.min):extend(b.max)
+end
+
+--- Get size of bounding box as a vector 
+-- @tparam bound2 a bound
+-- @treturn vec2 Vector spanning min to max points
+function bound2.size(a)
+	return a.max - a.min
+end
+
+--- Resize bounding box from minimum corner
+-- @tparam bound2 a bound
+-- @tparam vec2 new size
+-- @treturn bound2 resized bound
+function bound2.with_size(a, size)
+	return bound2.new(a.min, a.min + size)
+end
+
+--- Get half-size of bounding box as a vector. A more correct term for this is probably "apothem"
+-- @tparam bound2 a bound
+-- @treturn vec2 Vector spanning center to max point
+function bound2.radius(a)
+	return a:size()/2
+end
+
+--- Get center of bounding box
+-- @tparam bound2 a bound
+-- @treturn bound2 Point in center of bound
+function bound2.center(a)
+	return (a.min + a.max)/2
+end
+
+--- Move bounding box to new center
+-- @tparam bound2 a bound
+-- @tparam vec2 new center
+-- @treturn bound2 Bound with same size as input but different center
+function bound2.with_center(a, center)
+	return bound2.offset(a, center - a:center())
+end
+
+--- Resize bounding box from center
+-- @tparam bound2 a bound
+-- @tparam vec2 new size
+-- @treturn bound2 resized bound
+function bound2.with_size_centered(a, size)
+	local center = a:center()
+	local rad = size/2
+	return bound2.new(center - rad, center + rad)
+end
+
+--- Convert possibly-invalid bounding box to valid one
+-- @tparam bound2 a bound
+-- @treturn bound2 bound with all components corrected for min-max property
+function bound2.check(a)
+	if a.min.x > a.max.x or a.min.y > a.max.y then
+		return bound2.new(vec2.component_min(a.min, a.max), vec2.component_max(a.min, a.max))
+	end
+	return a
+end
+
+--- Shrink bounding box with fixed margin
+-- @tparam bound2 a bound
+-- @tparam vec2 a margin
+-- @treturn bound2 bound with margin subtracted from all edges. May not be valid, consider calling check()
+function bound2.inset(a, v)
+	return bound2.new(a.min + v, a.max - v)
+end
+
+--- Expand bounding box with fixed margin
+-- @tparam bound2 a bound
+-- @tparam vec2 a margin
+-- @treturn bound2 bound with margin added to all edges. May not be valid, consider calling check()
+function bound2.outset(a, v)
+	return bound2.new(a.min - v, a.max + v)
+end
+
+--- Offset bounding box
+-- @tparam bound2 a bound
+-- @tparam vec2 offset
+-- @treturn bound2 bound with same size, but position moved by offset
+function bound2.offset(a, v)
+	return bound2.new(a.min + v, a.max + v)
+end
+
+--- Test if point in bound
+-- @tparam bound2 a bound
+-- @tparam vec2 point to test
+-- @treturn boolean true if point in bounding box
+function bound2.contains(a, v)
+	return a.min.x <= v.x and a.min.y <= v.y
+	   and a.max.x >= v.x and a.max.y >= v.y
+end
+
+-- Round all components of all vectors to nearest int (or other precision).
+-- @tparam vec3 a bound to round.
+-- @tparam precision Digits after the decimal (round number if unspecified)
+-- @treturn vec3 Rounded bound
+function bound2.round(a, precision)
+	return bound2.new(a.min:round(precision), a.max:round(precision))
+end
+
+--- Return a formatted string.
+-- @tparam bound2 a bound to be turned into a string
+-- @treturn string formatted
+function bound2.to_string(a)
+	return string.format("(%s-%s)", a.min, a.max)
+end
+
+bound2_mt.__index    = bound2
+bound2_mt.__tostring = bound2.to_string
+
+function bound2_mt.__call(_, a, b)
+	return bound2.new(a, b)
+end
+
+if status then
+	xpcall(function() -- Allow this to silently fail; assume failure means someone messed with package.loaded
+		ffi.metatype(new, bound2_mt)
+	end, function() end)
+end
+
+return setmetatable({}, bound2_mt)

libs/cpml/bound3.lua

@@ -0,0 +1,199 @@
+--- A 3-component axis-aligned bounding box.
+-- @module bound3
+
+local modules = (...):gsub('%.[^%.]+$', '') .. "."
+local vec3    = require(modules .. "vec3")
+
+local bound3    = {}
+local bound3_mt = {}
+
+-- Private constructor.
+local function new(min, max)
+	return setmetatable({
+		min=min, -- min: vec3, minimum value for each component 
+		max=max  -- max: vec3, maximum value for each component
+	}, bound3_mt)
+end
+
+-- Do the check to see if JIT is enabled. If so use the optimized FFI structs.
+local status, ffi
+if type(jit) == "table" and jit.status() then
+	status, ffi = pcall(require, "ffi")
+	if status then
+		ffi.cdef "typedef struct { cpml_vec3 min, max; } cpml_bound3;"
+		new = ffi.typeof("cpml_bound3")
+	end
+end
+
+bound3.zero = new(vec3.zero, vec3.zero)
+
+--- The public constructor.
+-- @param min Can be of two types: </br>
+-- vec3 min, minimum value for each component
+-- nil Create bound at single point 0,0,0
+-- @tparam vec3 max, maximum value for each component
+-- @treturn bound3 out
+function bound3.new(min, max)
+	if min and max then
+		return new(min:clone(), max:clone())
+	elseif min or max then
+		error("Unexpected nil argument to bound3.new")
+	else
+		return new(vec3.zero, vec3.zero)
+	end
+end
+
+--- Clone a bound.
+-- @tparam bound3 a bound to be cloned
+-- @treturn bound3 out
+function bound3.clone(a)
+	return new(a.min, a.max)
+end
+
+--- Construct a bound covering one or two points 
+-- @tparam vec3 a Any vector
+-- @tparam vec3 b Any second vector (optional)
+-- @treturn vec3 Minimum bound containing the given points
+function bound3.at(a, b) -- "bounded by". b may be nil
+	if b then
+		return bound3.new(a,b):check()
+	else
+		return bound3.zero:with_center(a)
+	end
+end
+
+--- Extend bound to include point
+-- @tparam bound3 a bound
+-- @tparam vec3 point to include
+-- @treturn bound3 Bound covering current min, current max and new point
+function bound3.extend(a, center)
+	return bound3.new(a.min:component_min(center), a.max:component_max(center))
+end
+
+--- Extend bound to entirety of other bound
+-- @tparam bound3 a bound
+-- @tparam bound3 bound to cover
+-- @treturn bound3 Bound covering current min and max of each bound in the pair
+function bound3.extend_bound(a, b)
+	return a:extend(b.min):extend(b.max)
+end
+
+--- Get size of bounding box as a vector 
+-- @tparam bound3 a bound
+-- @treturn vec3 Vector spanning min to max points
+function bound3.size(a)
+	return a.max - a.min
+end
+
+--- Resize bounding box from minimum corner
+-- @tparam bound3 a bound
+-- @tparam vec3 new size
+-- @treturn bound3 resized bound
+function bound3.with_size(a, size)
+	return bound3.new(a.min, a.min + size)
+end
+
+--- Get half-size of bounding box as a vector. A more correct term for this is probably "apothem"
+-- @tparam bound3 a bound
+-- @treturn vec3 Vector spanning center to max point
+function bound3.radius(a)
+	return a:size()/2
+end
+
+--- Get center of bounding box
+-- @tparam bound3 a bound
+-- @treturn bound3 Point in center of bound
+function bound3.center(a)
+	return (a.min + a.max)/2
+end
+
+--- Move bounding box to new center
+-- @tparam bound3 a bound
+-- @tparam vec3 new center
+-- @treturn bound3 Bound with same size as input but different center
+function bound3.with_center(a, center)
+	return bound3.offset(a, center - a:center())
+end
+
+--- Resize bounding box from center
+-- @tparam bound3 a bound
+-- @tparam vec3 new size
+-- @treturn bound3 resized bound
+function bound3.with_size_centered(a, size)
+	local center = a:center()
+	local rad = size/2
+	return bound3.new(center - rad, center + rad)
+end
+
+--- Convert possibly-invalid bounding box to valid one
+-- @tparam bound3 a bound
+-- @treturn bound3 bound with all components corrected for min-max property
+function bound3.check(a)
+	if a.min.x > a.max.x or a.min.y > a.max.y or a.min.z > a.max.z then
+		return bound3.new(vec3.component_min(a.min, a.max), vec3.component_max(a.min, a.max))
+	end
+	return a
+end
+
+--- Shrink bounding box with fixed margin
+-- @tparam bound3 a bound
+-- @tparam vec3 a margin
+-- @treturn bound3 bound with margin subtracted from all edges. May not be valid, consider calling check()
+function bound3.inset(a, v)
+	return bound3.new(a.min + v, a.max - v)
+end
+
+--- Expand bounding box with fixed margin
+-- @tparam bound3 a bound
+-- @tparam vec3 a margin
+-- @treturn bound3 bound with margin added to all edges. May not be valid, consider calling check()
+function bound3.outset(a, v)
+	return bound3.new(a.min - v, a.max + v)
+end
+
+--- Offset bounding box
+-- @tparam bound3 a bound
+-- @tparam vec3 offset
+-- @treturn bound3 bound with same size, but position moved by offset
+function bound3.offset(a, v)
+	return bound3.new(a.min + v, a.max + v)
+end
+
+--- Test if point in bound
+-- @tparam bound3 a bound
+-- @tparam vec3 point to test
+-- @treturn boolean true if point in bounding box
+function bound3.contains(a, v)
+	return a.min.x <= v.x and a.min.y <= v.y and a.min.z <= v.z
+	   and a.max.x >= v.x and a.max.y >= v.y and a.max.z >= v.z
+end
+
+-- Round all components of all vectors to nearest int (or other precision).
+-- @tparam vec3 a bound to round.
+-- @tparam precision Digits after the decimal (round number if unspecified)
+-- @treturn vec3 Rounded bound
+function bound3.round(a, precision)
+	return bound3.new(a.min:round(precision), a.max:round(precision))
+end
+
+--- Return a formatted string.
+-- @tparam bound3 a bound to be turned into a string
+-- @treturn string formatted
+function bound3.to_string(a)
+	return string.format("(%s-%s)", a.min, a.max)
+end
+
+bound3_mt.__index    = bound3
+bound3_mt.__tostring = bound3.to_string
+
+function bound3_mt.__call(_, a, b)
+	return bound3.new(a, b)
+end
+
+if status then
+	xpcall(function() -- Allow this to silently fail; assume failure means someone messed with package.loaded
+		ffi.metatype(new, bound3_mt)
+	end, function() end)
+end
+
+return setmetatable({}, bound3_mt)

libs/cpml/bvh.lua

@@ -0,0 +1,557 @@
+-- https://github.com/benraziel/bvh-tree
+
+--- BVH Tree
+-- @module bvh
+
+local modules   = (...):gsub('%.[^%.]+$', '') .. "."
+local intersect = require(modules .. "intersect")
+local vec3      = require(modules .. "vec3")
+local EPSILON   = 1e-6
+local BVH       = {}
+local BVHNode   = {}
+local Node
+
+BVH.__index     = BVH
+BVHNode.__index = BVHNode
+
+local function new(triangles, maxTrianglesPerNode)
+	local tree = setmetatable({}, BVH)
+	local trianglesArray = {}
+
+	for _, triangle in ipairs(triangles) do
+		local p1 = triangle[1]
+		local p2 = triangle[2]
+		local p3 = triangle[3]
+
+		table.insert(trianglesArray, p1.x or p1[1])
+		table.insert(trianglesArray, p1.y or p1[2])
+		table.insert(trianglesArray, p1.z or p1[3])
+
+		table.insert(trianglesArray, p2.x or p2[1])
+		table.insert(trianglesArray, p2.y or p2[2])
+		table.insert(trianglesArray, p2.z or p2[3])
+
+		table.insert(trianglesArray, p3.x or p3[1])
+		table.insert(trianglesArray, p3.y or p3[2])
+		table.insert(trianglesArray, p3.z or p3[3])
+	end
+
+	tree._trianglesArray      = trianglesArray
+	tree._maxTrianglesPerNode = maxTrianglesPerNode or 10
+	tree._bboxArray           = tree.calcBoundingBoxes(trianglesArray)
+
+	-- clone a helper array
+	tree._bboxHelper = {}
+	for _, bbox in ipairs(tree._bboxArray) do
+		table.insert(tree._bboxHelper, bbox)
+	end
+
+	-- create the root node, add all the triangles to it
+	local triangleCount = #triangles
+	local extents = tree:calcExtents(1, triangleCount, EPSILON)
+	tree._rootNode = Node(extents[1], extents[2], 1, triangleCount, 1)
+
+	tree._nodes_to_split = { tree._rootNode }
+	while #tree._nodes_to_split > 0 do
+		local node = table.remove(tree._nodes_to_split)
+		tree:splitNode(node)
+	end
+	return tree
+end
+
+function BVH:intersectAABB(aabb)
+	local nodesToIntersect             = { self._rootNode }
+	local trianglesInIntersectingNodes = {} -- a list of nodes that intersect the ray (according to their bounding box)
+	local intersectingTriangles        = {}
+
+	-- go over the BVH tree, and extract the list of triangles that lie in nodes that intersect the box.
+	-- note: these triangles may not intersect the box themselves
+	while #nodesToIntersect > 0 do
+		local node = table.remove(nodesToIntersect)
+
+		local node_aabb = {
+			min = node._extentsMin,
+			max = node._extentsMax
+		}
+
+		if intersect.aabb_aabb(aabb, node_aabb) then
+			if node._node0 then
+				table.insert(nodesToIntersect, node._node0)
+			end
+
+			if node._node1 then
+				table.insert(nodesToIntersect, node._node1)
+			end
+
+			for i=node._startIndex, node._endIndex do
+				table.insert(trianglesInIntersectingNodes, self._bboxArray[1+(i-1)*7])
+			end
+		end
+	end
+
+	-- insert all node triangles, don't bother being more specific yet.
+	local triangle = { vec3(), vec3(), vec3() }
+
+	for i=1, #trianglesInIntersectingNodes do
+		local triIndex = trianglesInIntersectingNodes[i]
+
+		-- print(triIndex, #self._trianglesArray)
+		triangle[1].x = self._trianglesArray[1+(triIndex-1)*9]
+		triangle[1].y = self._trianglesArray[1+(triIndex-1)*9+1]
+		triangle[1].z = self._trianglesArray[1+(triIndex-1)*9+2]
+		triangle[2].x = self._trianglesArray[1+(triIndex-1)*9+3]
+		triangle[2].y = self._trianglesArray[1+(triIndex-1)*9+4]
+		triangle[2].z = self._trianglesArray[1+(triIndex-1)*9+5]
+		triangle[3].x = self._trianglesArray[1+(triIndex-1)*9+6]
+		triangle[3].y = self._trianglesArray[1+(triIndex-1)*9+7]
+		triangle[3].z = self._trianglesArray[1+(triIndex-1)*9+8]
+
+		table.insert(intersectingTriangles, {
+			triangle             = { triangle[1]:clone(), triangle[2]:clone(), triangle[3]:clone() },
+			triangleIndex        = triIndex
+		})
+	end
+
+	return intersectingTriangles
+end
+
+function BVH:intersectRay(rayOrigin, rayDirection, backfaceCulling)
+	local nodesToIntersect             = { self._rootNode }
+	local trianglesInIntersectingNodes = {} -- a list of nodes that intersect the ray (according to their bounding box)
+	local intersectingTriangles        = {}
+
+	local invRayDirection = vec3(
+		1 / rayDirection.x,
+		1 / rayDirection.y,
+		1 / rayDirection.z
+	)
+
+	-- go over the BVH tree, and extract the list of triangles that lie in nodes that intersect the ray.
+	-- note: these triangles may not intersect the ray themselves
+	while #nodesToIntersect > 0 do
+		local node = table.remove(nodesToIntersect)
+
+		if BVH.intersectNodeBox(rayOrigin, invRayDirection, node) then
+			if node._node0 then
+				table.insert(nodesToIntersect, node._node0)
+			end
+
+			if node._node1 then
+				table.insert(nodesToIntersect, node._node1)
+			end
+
+			for i=node._startIndex, node._endIndex do
+				table.insert(trianglesInIntersectingNodes, self._bboxArray[1+(i-1)*7])
+			end
+		end
+	end
+
+	-- go over the list of candidate triangles, and check each of them using ray triangle intersection
+	local triangle = { vec3(), vec3(), vec3() }
+	local ray      = {
+		position  = vec3(rayOrigin.x,    rayOrigin.y,    rayOrigin.z),
+		direction = vec3(rayDirection.x, rayDirection.y, rayDirection.z)
+	}
+
+	for i=1, #trianglesInIntersectingNodes do
+		local triIndex = trianglesInIntersectingNodes[i]
+
+		-- print(triIndex, #self._trianglesArray)
+		triangle[1].x = self._trianglesArray[1+(triIndex-1)*9]
+		triangle[1].y = self._trianglesArray[1+(triIndex-1)*9+1]
+		triangle[1].z = self._trianglesArray[1+(triIndex-1)*9+2]
+		triangle[2].x = self._trianglesArray[1+(triIndex-1)*9+3]
+		triangle[2].y = self._trianglesArray[1+(triIndex-1)*9+4]
+		triangle[2].z = self._trianglesArray[1+(triIndex-1)*9+5]
+		triangle[3].x = self._trianglesArray[1+(triIndex-1)*9+6]
+		triangle[3].y = self._trianglesArray[1+(triIndex-1)*9+7]
+		triangle[3].z = self._trianglesArray[1+(triIndex-1)*9+8]
+
+		local intersectionPoint, intersectionDistance = intersect.ray_triangle(ray, triangle, backfaceCulling)
+
+		if intersectionPoint then
+			table.insert(intersectingTriangles, {
+				triangle             = { triangle[1]:clone(), triangle[2]:clone(), triangle[3]:clone() },
+				triangleIndex        = triIndex,
+				intersectionPoint    = intersectionPoint,
+				intersectionDistance = intersectionDistance
+			})
+		end
+	end
+
+	return intersectingTriangles
+end
+
+function BVH.calcBoundingBoxes(trianglesArray)
+	local p1x, p1y, p1z
+	local p2x, p2y, p2z
+	local p3x, p3y, p3z
+	local minX, minY, minZ
+	local maxX, maxY, maxZ
+
+	local bboxArray = {}
+
+	for i=1, #trianglesArray / 9 do
+		p1x = trianglesArray[1+(i-1)*9]
+		p1y = trianglesArray[1+(i-1)*9+1]
+		p1z = trianglesArray[1+(i-1)*9+2]
+		p2x = trianglesArray[1+(i-1)*9+3]
+		p2y = trianglesArray[1+(i-1)*9+4]
+		p2z = trianglesArray[1+(i-1)*9+5]
+		p3x = trianglesArray[1+(i-1)*9+6]
+		p3y = trianglesArray[1+(i-1)*9+7]
+		p3z = trianglesArray[1+(i-1)*9+8]
+
+		minX = math.min(p1x, p2x, p3x)
+		minY = math.min(p1y, p2y, p3y)
+		minZ = math.min(p1z, p2z, p3z)
+		maxX = math.max(p1x, p2x, p3x)
+		maxY = math.max(p1y, p2y, p3y)
+		maxZ = math.max(p1z, p2z, p3z)
+
+		BVH.setBox(bboxArray, i, i, minX, minY, minZ, maxX, maxY, maxZ)
+	end
+
+	return bboxArray
+end
+
+function BVH:calcExtents(startIndex, endIndex, expandBy)
+	expandBy = expandBy or 0
+
+	if startIndex > endIndex then
+		return { vec3(), vec3() }
+	end
+
+	local minX =  math.huge
+	local minY =  math.huge
+	local minZ =  math.huge
+	local maxX = -math.huge
+	local maxY = -math.huge
+	local maxZ = -math.huge
+
+	for i=startIndex, endIndex do
+		minX = math.min(self._bboxArray[1+(i-1)*7+1], minX)
+		minY = math.min(self._bboxArray[1+(i-1)*7+2], minY)
+		minZ = math.min(self._bboxArray[1+(i-1)*7+3], minZ)
+		maxX = math.max(self._bboxArray[1+(i-1)*7+4], maxX)
+		maxY = math.max(self._bboxArray[1+(i-1)*7+5], maxY)
+		maxZ = math.max(self._bboxArray[1+(i-1)*7+6], maxZ)
+	end
+
+	return {
+		vec3(minX - expandBy, minY - expandBy, minZ - expandBy),
+		vec3(maxX + expandBy, maxY + expandBy, maxZ + expandBy)
+	}
+end
+
+function BVH:splitNode(node)
+	local num_elements = node:elementCount()
+	if (num_elements <= self._maxTrianglesPerNode) or (num_elements <= 0) then
+		return
+	end
+
+	local startIndex = node._startIndex
+	local endIndex   = node._endIndex
+
+	local leftNode  = { {},{},{} }
+	local rightNode = { {},{},{} }
+	local extentCenters = { node:centerX(), node:centerY(), node:centerZ() }
+
+	local extentsLength = {
+		node._extentsMax.x - node._extentsMin.x,
+		node._extentsMax.y - node._extentsMin.y,
+		node._extentsMax.z - node._extentsMin.z
+	}
+
+	local objectCenter = {}
+	for i=startIndex, endIndex do
+		objectCenter[1] = (self._bboxArray[1+(i-1)*7+1] + self._bboxArray[1+(i-1)*7+4]) * 0.5 -- center = (min + max) / 2
+		objectCenter[2] = (self._bboxArray[1+(i-1)*7+2] + self._bboxArray[1+(i-1)*7+5]) * 0.5 -- center = (min + max) / 2
+		objectCenter[3] = (self._bboxArray[1+(i-1)*7+3] + self._bboxArray[1+(i-1)*7+6]) * 0.5 -- center = (min + max) / 2
+
+		for j=1, 3 do
+			if objectCenter[j] < extentCenters[j] then
+				table.insert(leftNode[j], i)
+			else
+				table.insert(rightNode[j], i)
+			end
+		end
+	end
+
+	-- check if we couldn't split the node by any of the axes (x, y or z). halt
+	-- here, dont try to split any more (cause it will always fail, and we'll
+	-- enter an infinite loop
+	local splitFailed = {
+		#leftNode[1] == 0 or #rightNode[1] == 0,
+		#leftNode[2] == 0 or #rightNode[2] == 0,
+		#leftNode[3] == 0 or #rightNode[3] == 0
+	}
+
+	if splitFailed[1] and splitFailed[2] and splitFailed[3] then
+		return
+	end
+
+	-- choose the longest split axis. if we can't split by it, choose next best one.
+	local splitOrder = { 1, 2, 3 }
+	table.sort(splitOrder, function(a, b)
+		return extentsLength[a] > extentsLength[b]
+	end)
+
+	local leftElements
+	local rightElements
+
+	for i=1, 3 do
+		local candidateIndex = splitOrder[i]
+		if not splitFailed[candidateIndex] then
+			leftElements  = leftNode[candidateIndex]
+			rightElements = rightNode[candidateIndex]
+			break
+		end
+	end
+
+	-- sort the elements in range (startIndex, endIndex) according to which node they should be at
+	local node0Start = startIndex
+	local node1Start = node0Start + #leftElements
+	local node0End = node1Start - 1
+	local node1End = endIndex
+	local currElement
+
+	local helperPos = node._startIndex
+	local concatenatedElements = {}
+
+	for _, element in ipairs(leftElements) do
+		table.insert(concatenatedElements, element)
+	end
+
+	for _, element in ipairs(rightElements) do
+		table.insert(concatenatedElements, element)
+	end
+
+	-- print(#leftElements, #rightElements, #concatenatedElements)
+
+	for i=1, #concatenatedElements do
+		currElement = concatenatedElements[i]
+		BVH.copyBox(self._bboxArray, currElement, self._bboxHelper, helperPos)
+		helperPos = helperPos + 1
+	end
+
+	-- copy results back to main array
+	for i=1+(node._startIndex-1)*7, node._endIndex*7 do
+		self._bboxArray[i] = self._bboxHelper[i]
+	end
+
+	-- create 2 new nodes for the node we just split, and add links to them from the parent node
+	local node0Extents = self:calcExtents(node0Start, node0End, EPSILON)
+	local node1Extents = self:calcExtents(node1Start, node1End, EPSILON)
+
+	local node0 = Node(node0Extents[1], node0Extents[2], node0Start, node0End, node._level + 1)
+	local node1 = Node(node1Extents[1], node1Extents[2], node1Start, node1End, node._level + 1)
+
+	node._node0 = node0
+	node._node1 = node1
+	node:clearShapes()
+
+	-- add new nodes to the split queue
+	table.insert(self._nodes_to_split, node0)
+	table.insert(self._nodes_to_split, node1)
+end
+
+function BVH._calcTValues(minVal, maxVal, rayOriginCoord, invdir)
+	local res = { min=0, max=0 }
+
+	if invdir >= 0 then
+		res.min = ( minVal - rayOriginCoord ) * invdir
+		res.max = ( maxVal - rayOriginCoord ) * invdir
+	else
+		res.min = ( maxVal - rayOriginCoord ) * invdir
+		res.max = ( minVal - rayOriginCoord ) * invdir
+	end
+
+	return res
+end
+
+function BVH.intersectNodeBox(rayOrigin, invRayDirection, node)
+	local t  = BVH._calcTValues(node._extentsMin.x, node._extentsMax.x, rayOrigin.x, invRayDirection.x)
+	local ty = BVH._calcTValues(node._extentsMin.y, node._extentsMax.y, rayOrigin.y, invRayDirection.y)
+
+	if t.min > ty.max or ty.min > t.max then
+		return false
+	end
+
+	-- These lines also handle the case where tmin or tmax is NaN
+	-- (result of 0 * Infinity). x !== x returns true if x is NaN
+	if ty.min > t.min or t.min ~= t.min then
+		t.min = ty.min
+	end
+
+	if ty.max < t.max or t.max ~= t.max then
+		t.max = ty.max
+	end
+
+	local tz = BVH._calcTValues(node._extentsMin.z, node._extentsMax.z, rayOrigin.z, invRayDirection.z)
+
+	if t.min > tz.max or tz.min > t.max then
+		return false
+	end
+
+	if tz.min > t.min or t.min ~= t.min then
+		t.min = tz.min
+	end
+
+	if tz.max < t.max or t.max ~= t.max then
+		t.max = tz.max
+	end
+
+	--return point closest to the ray (positive side)
+	if t.max < 0 then
+		return false
+	end
+
+	return true
+end
+
+function BVH.setBox(bboxArray, pos, triangleId, minX, minY, minZ, maxX, maxY, maxZ)
+	bboxArray[1+(pos-1)*7]   = triangleId
+	bboxArray[1+(pos-1)*7+1] = minX
+	bboxArray[1+(pos-1)*7+2] = minY
+	bboxArray[1+(pos-1)*7+3] = minZ
+	bboxArray[1+(pos-1)*7+4] = maxX
+	bboxArray[1+(pos-1)*7+5] = maxY
+	bboxArray[1+(pos-1)*7+6] = maxZ
+end
+
+function BVH.copyBox(sourceArray, sourcePos, destArray, destPos)
+	destArray[1+(destPos-1)*7]   = sourceArray[1+(sourcePos-1)*7]
+	destArray[1+(destPos-1)*7+1] = sourceArray[1+(sourcePos-1)*7+1]
+	destArray[1+(destPos-1)*7+2] = sourceArray[1+(sourcePos-1)*7+2]
+	destArray[1+(destPos-1)*7+3] = sourceArray[1+(sourcePos-1)*7+3]
+	destArray[1+(destPos-1)*7+4] = sourceArray[1+(sourcePos-1)*7+4]
+	destArray[1+(destPos-1)*7+5] = sourceArray[1+(sourcePos-1)*7+5]
+	destArray[1+(destPos-1)*7+6] = sourceArray[1+(sourcePos-1)*7+6]
+end
+
+function BVH.getBox(bboxArray, pos, outputBox)
+	outputBox.triangleId = bboxArray[1+(pos-1)*7]
+	outputBox.minX       = bboxArray[1+(pos-1)*7+1]
+	outputBox.minY       = bboxArray[1+(pos-1)*7+2]
+	outputBox.minZ       = bboxArray[1+(pos-1)*7+3]
+	outputBox.maxX       = bboxArray[1+(pos-1)*7+4]
+	outputBox.maxY       = bboxArray[1+(pos-1)*7+5]
+	outputBox.maxZ       = bboxArray[1+(pos-1)*7+6]
+end
+
+local function new_node(extentsMin, extentsMax, startIndex, endIndex, level)
+	return setmetatable({
+		_extentsMin = extentsMin,
+		_extentsMax = extentsMax,
+		_startIndex = startIndex,
+		_endIndex   = endIndex,
+		_level      = level
+		--_node0    = nil
+		--_node1    = nil
+	}, BVHNode)
+end
+
+function BVHNode:elementCount()
+	return (self._endIndex + 1) - self._startIndex
+end
+
+function BVHNode:centerX()
+	return (self._extentsMin.x + self._extentsMax.x) * 0.5
+end
+
+function BVHNode:centerY()
+	return (self._extentsMin.y + self._extentsMax.y) * 0.5
+end
+
+function BVHNode:centerZ()
+	return (self._extentsMin.z + self._extentsMax.z) * 0.5
+end
+
+function BVHNode:clearShapes()
+	self._startIndex =  0
+	self._endIndex   = -1
+end
+
+function BVHNode.ngSphereRadius(extentsMin, extentsMax)
+	local centerX = (extentsMin.x + extentsMax.x) * 0.5
+	local centerY = (extentsMin.y + extentsMax.y) * 0.5
+	local centerZ = (extentsMin.z + extentsMax.z) * 0.5
+
+	local extentsMinDistSqr =
+		(centerX - extentsMin.x) * (centerX - extentsMin.x) +
+		(centerY - extentsMin.y) * (centerY - extentsMin.y) +
+		(centerZ - extentsMin.z) * (centerZ - extentsMin.z)
+
+	local extentsMaxDistSqr =
+		(centerX - extentsMax.x) * (centerX - extentsMax.x) +
+		(centerY - extentsMax.y) * (centerY - extentsMax.y) +
+		(centerZ - extentsMax.z) * (centerZ - extentsMax.z)
+
+	return math.sqrt(math.max(extentsMinDistSqr, extentsMaxDistSqr))
+end
+
+--[[
+
+--- Draws node boundaries visually for debugging.
+-- @param cube Cube model to draw
+-- @param depth Used for recurcive calls to self method
+function OctreeNode:draw_bounds(cube, depth)
+	depth = depth or 0
+	local tint = depth / 7 -- Will eventually get values > 1. Color rounds to 1 automatically
+
+	love.graphics.setColor(tint * 255, 0, (1 - tint) * 255)
+	local m = mat4()
+		:translate(self.center)
+		:scale(vec3(self.adjLength, self.adjLength, self.adjLength))
+
+	love.graphics.updateMatrix("transform", m)
+	love.graphics.setWireframe(true)
+	love.graphics.draw(cube)
+	love.graphics.setWireframe(false)
+
+	for _, child in ipairs(self.children) do
+		child:draw_bounds(cube, depth + 1)
+	end
+
+	love.graphics.setColor(255, 255, 255)
+end
+
+--- Draws the bounds of all objects in the tree visually for debugging.
+-- @param cube Cube model to draw
+-- @param filter a function returning true or false to determine visibility.
+function OctreeNode:draw_objects(cube, filter)
+	local tint = self.baseLength / 20
+	love.graphics.setColor(0, (1 - tint) * 255, tint * 255, 63)
+
+	for _, object in ipairs(self.objects) do
+		if filter and filter(object.data) or not filter then
+			local m = mat4()
+				:translate(object.bounds.center)
+				:scale(object.bounds.size)
+
+			love.graphics.updateMatrix("transform", m)
+			love.graphics.draw(cube)
+		end
+	end
+
+	for _, child in ipairs(self.children) do
+		child:draw_objects(cube, filter)
+	end
+
+	love.graphics.setColor(255, 255, 255)
+end
+
+--]]
+
+Node = setmetatable({
+	new = new_node
+}, {
+	__call = function(_, ...) return new_node(...) end
+})
+
+return setmetatable({
+	new = new
+}, {
+	__call = function(_, ...) return new(...) end
+})

libs/cpml/color.lua

@@ -0,0 +1,400 @@
+--- Color utilities
+-- @module color
+
+local modules  = (...):gsub('%.[^%.]+$', '') .. "."
+local constants = require(modules .. "constants")
+local utils    = require(modules .. "utils")
+local precond  = require(modules .. "_private_precond")
+local color    = {}
+local color_mt = {}
+
+local function new(r, g, b, a)
+	local c = { r, g, b, a }
+	c._c = c
+	return setmetatable(c, color_mt)
+end
+
+-- HSV utilities (adapted from http://www.cs.rit.edu/~ncs/color/t_convert.html)
+-- hsv_to_color(hsv)
+-- Converts a set of HSV values to a color. hsv is a table.
+-- See also: hsv(h, s, v)
+local function hsv_to_color(hsv)
+	local i
+	local f, q, p, t
+	local h, s, v
+	local a = hsv[4] or 1
+	s = hsv[2]
+	v = hsv[3]
+
+	if s == 0 then
+		return new(v, v, v, a)
+	end
+
+	h = hsv[1] * 6 -- sector 0 to 5
+
+	i = math.floor(h)
+	f = h - i -- factorial part of h
+	p = v * (1-s)
+	q = v * (1-s*f)
+	t = v * (1-s*(1-f))
+
+	if     i == 0 then return new(v, t, p, a)
+	elseif i == 1 then return new(q, v, p, a)
+	elseif i == 2 then return new(p, v, t, a)
+	elseif i == 3 then return new(p, q, v, a)
+	elseif i == 4 then return new(t, p, v, a)
+	else               return new(v, p, q, a)
+	end
+end
+
+-- color_to_hsv(c)
+-- Takes in a normal color and returns a table with the HSV values.
+local function color_to_hsv(c)
+	local r = c[1]
+	local g = c[2]
+	local b = c[3]
+	local a = c[4] or 1
+	local h, s, v
+
+	local min = math.min(r, g, b)
+	local max = math.max(r, g, b)
+	v = max
+
+	local delta = max - min
+
+	-- black, nothing else is really possible here.
+	if min == 0 and max == 0 then
+		return { 0, 0, 0, a }
+	end
+
+	if max ~= 0 then
+		s = delta / max
+	else
+		-- r = g = b = 0 s = 0, v is undefined
+		s = 0
+		h = -1
+		return { h, s, v, 1 }
+	end
+
+	-- Prevent division by zero.
+	if delta == 0 then
+		delta = constants.DBL_EPSILON
+	end
+
+	if r == max then
+		h = ( g - b ) / delta     -- yellow/magenta
+	elseif g == max then
+		h = 2 + ( b - r ) / delta -- cyan/yellow
+	else
+		h = 4 + ( r - g ) / delta -- magenta/cyan
+	end
+
+	h = h / 6 -- normalize from segment 0..5
+
+	if h < 0 then
+		h = h + 1
+	end
+
+	return { h, s, v, a }
+end
+
+--- The public constructor.
+-- @param x Can be of three types: </br>
+-- number red component 0-1
+-- table {r, g, b, a}
+-- nil for {0,0,0,0}
+-- @tparam number g Green component 0-1
+-- @tparam number b Blue component 0-1
+-- @tparam number a Alpha component 0-1
+-- @treturn color out
+function color.new(r, g, b, a)
+	-- number, number, number, number
+	if r and g and b and a then
+		precond.typeof(r, "number", "new: Wrong argument type for r")
+		precond.typeof(g, "number", "new: Wrong argument type for g")
+		precond.typeof(b, "number", "new: Wrong argument type for b")
+		precond.typeof(a, "number", "new: Wrong argument type for a")
+
+		return new(r, g, b, a)
+
+	-- {r, g, b, a}
+	elseif type(r) == "table" then
+		local rr, gg, bb, aa = r[1], r[2], r[3], r[4]
+		precond.typeof(rr, "number", "new: Wrong argument type for r")
+		precond.typeof(gg, "number", "new: Wrong argument type for g")
+		precond.typeof(bb, "number", "new: Wrong argument type for b")
+		precond.typeof(aa, "number", "new: Wrong argument type for a")
+
+		return new(rr, gg, bb, aa)
+	end
+
+	return new(0, 0, 0, 0)
+end
+
+--- Convert hue,saturation,value table to color object.
+-- @tparam table hsva {hue 0-1, saturation 0-1, value 0-1, alpha 0-1}
+-- @treturn color out
+color.hsv_to_color_table = hsv_to_color
+
+--- Convert color to hue,saturation,value table
+-- @tparam color in
+-- @treturn table hsva {hue 0-1, saturation 0-1, value 0-1, alpha 0-1}
+color.color_to_hsv_table = color_to_hsv
+
+--- Convert hue,saturation,value to color object.
+-- @tparam number h hue 0-1
+-- @tparam number s saturation 0-1
+-- @tparam number v value 0-1
+-- @treturn color out
+function color.from_hsv(h, s, v)
+	return hsv_to_color { h, s, v }
+end
+
+--- Convert hue,saturation,value to color object.
+-- @tparam number h hue 0-1
+-- @tparam number s saturation 0-1
+-- @tparam number v value 0-1
+-- @tparam number a alpha 0-1
+-- @treturn color out
+function color.from_hsva(h, s, v, a)
+	return hsv_to_color { h, s, v, a }
+end
+
+--- Invert a color.
+-- @tparam color to invert
+-- @treturn color out
+function color.invert(c)
+	return new(1 - c[1], 1 - c[2], 1 - c[3], c[4])
+end
+
+--- Lighten a color by a component-wise fixed amount (alpha unchanged)
+-- @tparam color to lighten
+-- @tparam number amount to increase each component by, 0-1 scale
+-- @treturn color out
+function color.lighten(c, v)
+	return new(
+		utils.clamp(c[1] + v, 0, 1),
+		utils.clamp(c[2] + v, 0, 1),
+		utils.clamp(c[3] + v, 0, 1),
+		c[4]
+	)
+end
+
+--- Interpolate between two colors.
+-- @tparam color at start
+-- @tparam color at end
+-- @tparam number s in 0-1 progress between the two colors
+-- @treturn color out
+function color.lerp(a, b, s)
+	return a + s * (b - a)
+end
+
+--- Unpack a color into individual components in 0-1.
+-- @tparam color to unpack
+-- @treturn number r in 0-1
+-- @treturn number g in 0-1
+-- @treturn number b in 0-1
+-- @treturn number a in 0-1
+function color.unpack(c)
+	return c[1], c[2], c[3], c[4]
+end
+
+--- Unpack a color into individual components in 0-255.
+-- @tparam color to unpack
+-- @treturn number r in 0-255
+-- @treturn number g in 0-255
+-- @treturn number b in 0-255
+-- @treturn number a in 0-255
+function color.as_255(c)
+	return c[1] * 255, c[2] * 255, c[3] * 255, c[4] * 255
+end
+
+--- Darken a color by a component-wise fixed amount (alpha unchanged)
+-- @tparam color to darken
+-- @tparam number amount to decrease each component by, 0-1 scale
+-- @treturn color out
+function color.darken(c, v)
+	return new(
+		utils.clamp(c[1] - v, 0, 1),
+		utils.clamp(c[2] - v, 0, 1),
+		utils.clamp(c[3] - v, 0, 1),
+		c[4]
+	)
+end
+
+--- Multiply a color's components by a value (alpha unchanged)
+-- @tparam color to multiply
+-- @tparam number to multiply each component by
+-- @treturn color out
+function color.multiply(c, v)
+	local t = color.new()
+	for i = 1, 3 do
+		t[i] = c[i] * v
+	end
+
+	t[4] = c[4]
+	return t
+end
+
+-- directly set alpha channel
+-- @tparam color to alter
+-- @tparam number new alpha 0-1
+-- @treturn color out
+function color.alpha(c, v)
+	local t = color.new()
+	for i = 1, 3 do
+		t[i] = c[i]
+	end
+
+	t[4] = v
+	return t
+end
+
+--- Multiply a color's alpha by a value
+-- @tparam color to multiply
+-- @tparam number to multiply alpha by
+-- @treturn color out
+function color.opacity(c, v)
+	local t = color.new()
+	for i = 1, 3 do
+		t[i] = c[i]
+	end
+
+	t[4] = c[4] * v
+	return t
+end
+
+--- Set a color's hue (saturation, value, alpha unchanged)
+-- @tparam color to alter
+-- @tparam hue to set 0-1
+-- @treturn color out
+function color.hue(col, hue)
+	local c = color_to_hsv(col)
+	c[1] = (hue + 1) % 1
+	return hsv_to_color(c)
+end
+
+--- Set a color's saturation (hue, value, alpha unchanged)
+-- @tparam color to alter
+-- @tparam saturation to set 0-1
+-- @treturn color out
+function color.saturation(col, percent)
+	local c = color_to_hsv(col)
+	c[2] = utils.clamp(percent, 0, 1)
+	return hsv_to_color(c)
+end
+
+--- Set a color's value (saturation, hue, alpha unchanged)
+-- @tparam color to alter
+-- @tparam value to set 0-1
+-- @treturn color out
+function color.value(col, percent)
+	local c = color_to_hsv(col)
+	c[3] = utils.clamp(percent, 0, 1)
+	return hsv_to_color(c)
+end
+
+-- https://en.wikipedia.org/wiki/SRGB#From_sRGB_to_CIE_XYZ
+function color.gamma_to_linear(r, g, b, a)
+	local function convert(c)
+		if c > 1.0 then
+			return 1.0
+		elseif c < 0.0 then
+			return 0.0
+		elseif c <= 0.04045 then
+			return c / 12.92
+		else
+			return math.pow((c + 0.055) / 1.055, 2.4)
+		end
+	end
+
+	if type(r) == "table" then
+		local c = {}
+		for i = 1, 3 do
+			c[i] = convert(r[i])
+		end
+
+		c[4] = r[4]
+		return c
+	else
+		return convert(r), convert(g), convert(b), a or 1
+	end
+end
+
+-- https://en.wikipedia.org/wiki/SRGB#From_CIE_XYZ_to_sRGB
+function color.linear_to_gamma(r, g, b, a)
+	local function convert(c)
+		if c > 1.0 then
+			return 1.0
+		elseif c < 0.0 then
+			return 0.0
+		elseif c < 0.0031308 then
+			return c * 12.92
+		else
+			return 1.055 * math.pow(c, 0.41666) - 0.055
+		end
+	end
+
+	if type(r) == "table" then
+		local c = {}
+		for i = 1, 3 do
+			c[i] = convert(r[i])
+		end
+
+		c[4] = r[4]
+		return c
+	else
+		return convert(r), convert(g), convert(b), a or 1
+	end
+end
+
+--- Check if color is valid
+-- @tparam color to test
+-- @treturn boolean is color
+function color.is_color(a)
+	if type(a) ~= "table" then
+		return false
+	end
+
+	for i = 1, 4 do
+		if type(a[i]) ~= "number" then
+			return false
+		end
+	end
+
+	return true
+end
+
+--- Return a formatted string.
+-- @tparam color a color to be turned into a string
+-- @treturn string formatted
+function color.to_string(a)
+	return string.format("[ %3.0f, %3.0f, %3.0f, %3.0f ]", a[1], a[2], a[3], a[4])
+end
+
+color_mt.__index = color
+color_mt.__tostring = color.to_string
+
+function color_mt.__call(_, r, g, b, a)
+	return color.new(r, g, b, a)
+end
+
+function color_mt.__add(a, b)
+	return new(a[1] + b[1], a[2] + b[2], a[3] + b[3], a[4] + b[4])
+end
+
+function color_mt.__sub(a, b)
+	return new(a[1] - b[1], a[2] - b[2], a[3] - b[3], a[4] - b[4])
+end
+
+function color_mt.__mul(a, b)
+	if type(a) == "number" then
+		return new(a * b[1], a * b[2], a * b[3], a * b[4])
+	elseif type(b) == "number" then
+		return new(b * a[1], b * a[2], b * a[3], b * a[4])
+	else
+		return new(a[1] * b[1], a[2] * b[2], a[3] * b[3], a[4] * b[4])
+	end
+end
+
+return setmetatable({}, color_mt)

libs/cpml/constants.lua

@@ -0,0 +1,20 @@
+--- Various useful constants
+-- @module constants
+
+--- Constants
+-- @table constants
+-- @field FLT_EPSILON Floating point precision breaks down
+-- @field DBL_EPSILON Double-precise floating point precision breaks down
+-- @field DOT_THRESHOLD Close enough to 1 for interpolations to occur
+local constants = {}
+
+-- same as C's FLT_EPSILON
+constants.FLT_EPSILON = 1.19209290e-07
+
+-- same as C's DBL_EPSILON
+constants.DBL_EPSILON = 2.2204460492503131e-16
+
+-- used for quaternion.slerp
+constants.DOT_THRESHOLD = 0.9995
+
+return constants

libs/cpml/intersect.lua

@@ -0,0 +1,709 @@
+--- Various geometric intersections
+-- @module intersect
+
+local modules     = (...):gsub('%.[^%.]+$', '') .. "."
+local constants   = require(modules .. "constants")
+local mat4        = require(modules .. "mat4")
+local vec3        = require(modules .. "vec3")
+local utils       = require(modules .. "utils")
+local DBL_EPSILON = constants.DBL_EPSILON
+local sqrt        = math.sqrt
+local abs         = math.abs
+local min         = math.min
+local max         = math.max
+local intersect   = {}
+
+-- https://blogs.msdn.microsoft.com/rezanour/2011/08/07/barycentric-coordinates-and-point-in-triangle-tests/
+-- point       is a vec3
+-- triangle[1] is a vec3
+-- triangle[2] is a vec3
+-- triangle[3] is a vec3
+function intersect.point_triangle(point, triangle)
+	local u = triangle[2] - triangle[1]
+	local v = triangle[3] - triangle[1]
+	local w = point       - triangle[1]
+
+	local vw = v:cross(w)
+	local vu = v:cross(u)
+
+	if vw:dot(vu) < 0 then
+		return false
+	end
+
+	local uw = u:cross(w)
+	local uv = u:cross(v)
+
+	if uw:dot(uv) < 0 then
+		return false
+	end
+
+	local d = uv:len()
+	local r = vw:len() / d
+	local t = uw:len() / d
+
+	return r + t <= 1
+end
+
+-- point    is a vec3
+-- aabb.min is a vec3
+-- aabb.max is a vec3
+function intersect.point_aabb(point, aabb)
+	return
+		aabb.min.x <= point.x and
+		aabb.max.x >= point.x and
+		aabb.min.y <= point.y and
+		aabb.max.y >= point.y and
+		aabb.min.z <= point.z and
+		aabb.max.z >= point.z
+end
+
+-- point          is a vec3
+-- frustum.left   is a plane { a, b, c, d }
+-- frustum.right  is a plane { a, b, c, d }
+-- frustum.bottom is a plane { a, b, c, d }
+-- frustum.top    is a plane { a, b, c, d }
+-- frustum.near   is a plane { a, b, c, d }
+-- frustum.far    is a plane { a, b, c, d }
+function intersect.point_frustum(point, frustum)
+	local x, y, z = point:unpack()
+	local planes  = {
+		frustum.left,
+		frustum.right,
+		frustum.bottom,
+		frustum.top,
+		frustum.near,
+		frustum.far or false
+	}
+
+	-- Skip the last test for infinite projections, it'll never fail.
+	if not planes[6] then
+		table.remove(planes)
+	end
+
+	local dot
+	for i = 1, #planes do
+		dot = planes[i].a * x + planes[i].b * y + planes[i].c * z + planes[i].d
+		if dot <= 0 then
+			return false
+		end
+	end
+
+	return true
+end
+
+-- http://www.lighthouse3d.com/tutorials/maths/ray-triangle-intersection/
+-- ray.position  is a vec3
+-- ray.direction is a vec3
+-- triangle[1]   is a vec3
+-- triangle[2]   is a vec3
+-- triangle[3]   is a vec3
+-- backface_cull is a boolean (optional)
+function intersect.ray_triangle(ray, triangle, backface_cull)
+	local e1 = triangle[2] - triangle[1]
+	local e2 = triangle[3] - triangle[1]
+	local h  = ray.direction:cross(e2)
+	local a  = h:dot(e1)
+
+	-- if a is negative, ray hits the backface
+	if backface_cull and a < 0 then
+		return false
+	end
+
+	-- if a is too close to 0, ray does not intersect triangle
+	if abs(a) <= DBL_EPSILON then
+		return false
+	end
+
+	local f = 1 / a
+	local s = ray.position - triangle[1]
+	local u = s:dot(h) * f
+
+	-- ray does not intersect triangle
+	if u < 0 or u > 1 then
+		return false
+	end
+
+	local q = s:cross(e1)
+	local v = ray.direction:dot(q) * f
+
+	-- ray does not intersect triangle
+	if v < 0 or u + v > 1 then
+		return false
+	end
+
+	-- at this stage we can compute t to find out where
+	-- the intersection point is on the line
+	local t = q:dot(e2) * f
+
+	-- return position of intersection and distance from ray origin
+	if t >= DBL_EPSILON then
+		return ray.position + ray.direction * t, t
+	end
+
+	-- ray does not intersect triangle
+	return false
+end
+
+-- https://gamedev.stackexchange.com/questions/96459/fast-ray-sphere-collision-code
+-- ray.position    is a vec3
+-- ray.direction   is a vec3
+-- sphere.position is a vec3
+-- sphere.radius   is a number
+function intersect.ray_sphere(ray, sphere)
+	local offset = ray.position - sphere.position
+	local b = offset:dot(ray.direction)
+	local c = offset:dot(offset) - sphere.radius * sphere.radius
+
+	-- ray's position outside sphere (c > 0)
+	-- ray's direction pointing away from sphere (b > 0)
+	if c > 0 and b > 0 then
+		return false
+	end
+
+	local discr = b * b - c
+
+	-- negative discriminant
+	if discr < 0 then
+		return false
+	end
+
+	-- Clamp t to 0
+	local t = -b - sqrt(discr)
+	t = t < 0 and 0 or t
+
+	-- Return collision point and distance from ray origin
+	return ray.position + ray.direction * t, t
+end
+
+-- http://gamedev.stackexchange.com/a/18459
+-- ray.position  is a vec3
+-- ray.direction is a vec3
+-- aabb.min      is a vec3
+-- aabb.max      is a vec3
+function intersect.ray_aabb(ray, aabb)
+	local dir     = ray.direction:normalize()
+	local dirfrac = vec3(
+		1 / dir.x,
+		1 / dir.y,
+		1 / dir.z
+	)
+
+	local t1 = (aabb.min.x - ray.position.x) * dirfrac.x
+	local t2 = (aabb.max.x - ray.position.x) * dirfrac.x
+	local t3 = (aabb.min.y - ray.position.y) * dirfrac.y
+	local t4 = (aabb.max.y - ray.position.y) * dirfrac.y
+	local t5 = (aabb.min.z - ray.position.z) * dirfrac.z
+	local t6 = (aabb.max.z - ray.position.z) * dirfrac.z
+
+	local tmin = max(max(min(t1, t2), min(t3, t4)), min(t5, t6))
+	local tmax = min(min(max(t1, t2), max(t3, t4)), max(t5, t6))
+
+	-- ray is intersecting AABB, but whole AABB is behind us
+	if tmax < 0 then
+		return false
+	end
+
+	-- ray does not intersect AABB
+	if tmin > tmax then
+		return false
+	end
+
+	-- Return collision point and distance from ray origin
+	return ray.position + ray.direction * tmin, tmin
+end
+
+-- http://stackoverflow.com/a/23976134/1190664
+-- ray.position   is a vec3
+-- ray.direction  is a vec3
+-- plane.position is a vec3
+-- plane.normal   is a vec3
+function intersect.ray_plane(ray, plane)
+	local denom = plane.normal:dot(ray.direction)
+
+	-- ray does not intersect plane
+	if abs(denom) < DBL_EPSILON then
+		return false
+	end
+
+	-- distance of direction
+	local d = plane.position - ray.position
+	local t = d:dot(plane.normal) / denom
+
+	if t < DBL_EPSILON then
+		return false
+	end
+
+	-- Return collision point and distance from ray origin
+	return ray.position + ray.direction * t, t
+end
+
+function intersect.ray_capsule(ray, capsule)
+	local dist2, p1, p2 = intersect.closest_point_segment_segment(
+		ray.position,
+		ray.position + ray.direction * 1e10,
+		capsule.a,
+		capsule.b
+	)
+	if dist2 <= capsule.radius^2 then
+		return p1
+	end
+
+	return false
+end
+
+-- https://web.archive.org/web/20120414063459/http://local.wasp.uwa.edu.au/~pbourke//geometry/lineline3d/
+-- a[1] is a vec3
+-- a[2] is a vec3
+-- b[1] is a vec3
+-- b[2] is a vec3
+-- e    is a number
+function intersect.line_line(a, b, e)
+	-- new points
+	local p13 = a[1] - b[1]
+	local p43 = b[2] - b[1]
+	local p21 = a[2] - a[1]
+
+	-- if lengths are negative or too close to 0, lines do not intersect
+	if p43:len2() < DBL_EPSILON or p21:len2() < DBL_EPSILON then
+		return false
+	end
+
+	-- dot products
+	local d1343 = p13:dot(p43)
+	local d4321 = p43:dot(p21)
+	local d1321 = p13:dot(p21)
+	local d4343 = p43:dot(p43)
+	local d2121 = p21:dot(p21)
+	local denom = d2121 * d4343 - d4321 * d4321
+
+	-- if denom is too close to 0, lines do not intersect
+	if abs(denom) < DBL_EPSILON then
+		return false
+	end
+
+	local numer = d1343 * d4321 - d1321 * d4343
+	local mua   = numer / denom
+	local mub   = (d1343 + d4321 * mua) / d4343
+
+	-- return positions of intersection on each line
+	local out1 = a[1] + p21 * mua
+	local out2 = b[1] + p43 * mub
+	local dist = out1:dist(out2)
+
+	-- if distance of the shortest segment between lines is less than threshold
+	if e and dist > e then
+		return false
+	end
+
+	return { out1, out2 }, dist
+end
+
+-- a[1] is a vec3
+-- a[2] is a vec3
+-- b[1] is a vec3
+-- b[2] is a vec3
+-- e    is a number
+function intersect.segment_segment(a, b, e)
+	local c, d = intersect.line_line(a, b, e)
+
+	if c and ((
+		a[1].x <= c[1].x and
+		a[1].y <= c[1].y and
+		a[1].z <= c[1].z and
+		c[1].x <= a[2].x and
+		c[1].y <= a[2].y and
+		c[1].z <= a[2].z
+	) or (
+		a[1].x >= c[1].x and
+		a[1].y >= c[1].y and
+		a[1].z >= c[1].z and
+		c[1].x >= a[2].x and
+		c[1].y >= a[2].y and
+		c[1].z >= a[2].z
+	)) and ((
+		b[1].x <= c[2].x and
+		b[1].y <= c[2].y and
+		b[1].z <= c[2].z and
+		c[2].x <= b[2].x and
+		c[2].y <= b[2].y and
+		c[2].z <= b[2].z
+	) or (
+		b[1].x >= c[2].x and
+		b[1].y >= c[2].y and
+		b[1].z >= c[2].z and
+		c[2].x >= b[2].x and
+		c[2].y >= b[2].y and
+		c[2].z >= b[2].z
+	)) then
+		return c, d
+	end
+
+	-- segments do not intersect
+	return false
+end
+
+-- a.min is a vec3
+-- a.max is a vec3
+-- b.min is a vec3
+-- b.max is a vec3
+function intersect.aabb_aabb(a, b)
+	return
+		a.min.x <= b.max.x and
+		a.max.x >= b.min.x and
+		a.min.y <= b.max.y and
+		a.max.y >= b.min.y and
+		a.min.z <= b.max.z and
+		a.max.z >= b.min.z
+end
+
+-- aabb.position is a vec3
+-- aabb.extent   is a vec3 (half-size)
+-- obb.position  is a vec3
+-- obb.extent    is a vec3 (half-size)
+-- obb.rotation  is a mat4
+function intersect.aabb_obb(aabb, obb)
+	local a   = aabb.extent
+	local b   = obb.extent
+	local T   = obb.position - aabb.position
+	local rot = mat4():transpose(obb.rotation)
+	local B   = {}
+	local t
+
+	for i = 1, 3 do
+		B[i] = {}
+		for j = 1, 3 do
+			assert((i - 1) * 4 + j < 16 and (i - 1) * 4 + j > 0)
+			B[i][j] = abs(rot[(i - 1) * 4 + j]) + 1e-6
+		end
+	end
+
+	t = abs(T.x)
+	if not (t <= (b.x + a.x * B[1][1] + b.y * B[1][2] + b.z * B[1][3])) then return false end
+	t = abs(T.x * B[1][1] + T.y * B[2][1] + T.z * B[3][1])
+	if not (t <= (b.x + a.x * B[1][1] + a.y * B[2][1] + a.z * B[3][1])) then return false end
+	t = abs(T.y)
+	if not (t <= (a.y + b.x * B[2][1] + b.y * B[2][2] + b.z * B[2][3])) then return false end
+	t = abs(T.z)
+	if not (t <= (a.z + b.x * B[3][1] + b.y * B[3][2] + b.z * B[3][3])) then return false end
+	t = abs(T.x * B[1][2] + T.y * B[2][2] + T.z * B[3][2])
+	if not (t <= (b.y + a.x * B[1][2] + a.y * B[2][2] + a.z * B[3][2])) then return false end
+	t = abs(T.x * B[1][3] + T.y * B[2][3] + T.z * B[3][3])
+	if not (t <= (b.z + a.x * B[1][3] + a.y * B[2][3] + a.z * B[3][3])) then return false end
+	t = abs(T.z * B[2][1] - T.y * B[3][1])
+	if not (t <= (a.y * B[3][1] + a.z * B[2][1] + b.y * B[1][3] + b.z * B[1][2])) then return false end
+	t = abs(T.z * B[2][2] - T.y * B[3][2])
+	if not (t <= (a.y * B[3][2] + a.z * B[2][2] + b.x * B[1][3] + b.z * B[1][1])) then return false end
+	t = abs(T.z * B[2][3] - T.y * B[3][3])
+	if not (t <= (a.y * B[3][3] + a.z * B[2][3] + b.x * B[1][2] + b.y * B[1][1])) then return false end
+	t = abs(T.x * B[3][1] - T.z * B[1][1])
+	if not (t <= (a.x * B[3][1] + a.z * B[1][1] + b.y * B[2][3] + b.z * B[2][2])) then return false end
+	t = abs(T.x * B[3][2] - T.z * B[1][2])
+	if not (t <= (a.x * B[3][2] + a.z * B[1][2] + b.x * B[2][3] + b.z * B[2][1])) then return false end
+	t = abs(T.x * B[3][3] - T.z * B[1][3])
+	if not (t <= (a.x * B[3][3] + a.z * B[1][3] + b.x * B[2][2] + b.y * B[2][1])) then return false end
+	t = abs(T.y * B[1][1] - T.x * B[2][1])
+	if not (t <= (a.x * B[2][1] + a.y * B[1][1] + b.y * B[3][3] + b.z * B[3][2])) then return false end
+	t = abs(T.y * B[1][2] - T.x * B[2][2])
+	if not (t <= (a.x * B[2][2] + a.y * B[1][2] + b.x * B[3][3] + b.z * B[3][1])) then return false end
+	t = abs(T.y * B[1][3] - T.x * B[2][3])
+	if not (t <= (a.x * B[2][3] + a.y * B[1][3] + b.x * B[3][2] + b.y * B[3][1])) then return false end
+
+	-- https://gamedev.stackexchange.com/questions/24078/which-side-was-hit
+	-- Minkowski Sum
+	local wy = (aabb.extent * 2 + obb.extent * 2) * (aabb.position.y - obb.position.y)
+	local hx = (aabb.extent * 2 + obb.extent * 2) * (aabb.position.x - obb.position.x)
+
+	if wy.x > hx.x and wy.y > hx.y and wy.z > hx.z then
+		if wy.x > -hx.x and wy.y > -hx.y and wy.z > -hx.z then
+			return vec3(obb.rotation * {  0, -1, 0, 1 })
+		else
+			return vec3(obb.rotation * { -1,  0, 0, 1 })
+		end
+	else
+		if wy.x > -hx.x and wy.y > -hx.y and wy.z > -hx.z then
+			return vec3(obb.rotation * { 1, 0, 0, 1 })
+		else
+			return vec3(obb.rotation * { 0, 1, 0, 1 })
+		end
+	end
+end
+
+-- http://stackoverflow.com/a/4579069/1190664
+-- aabb.min        is a vec3
+-- aabb.max        is a vec3
+-- sphere.position is a vec3
+-- sphere.radius   is a number
+local axes = { "x", "y", "z" }
+function intersect.aabb_sphere(aabb, sphere)
+	local dist2 = sphere.radius ^ 2
+
+	for _, axis in ipairs(axes) do
+		local pos  = sphere.position[axis]
+		local amin = aabb.min[axis]
+		local amax = aabb.max[axis]
+
+		if pos < amin then
+			dist2 = dist2 - (pos - amin) ^ 2
+		elseif pos > amax then
+			dist2 = dist2 - (pos - amax) ^ 2
+		end
+	end
+
+	return dist2 > 0
+end
+
+-- aabb.min       is a vec3
+-- aabb.max       is a vec3
+-- frustum.left   is a plane { a, b, c, d }
+-- frustum.right  is a plane { a, b, c, d }
+-- frustum.bottom is a plane { a, b, c, d }
+-- frustum.top    is a plane { a, b, c, d }
+-- frustum.near   is a plane { a, b, c, d }
+-- frustum.far    is a plane { a, b, c, d }
+function intersect.aabb_frustum(aabb, frustum)
+	-- Indexed for the 'index trick' later
+	local box = {
+		aabb.min,
+		aabb.max
+	}
+
+	-- We have 6 planes defining the frustum, 5 if infinite.
+	local planes = {
+		frustum.left,
+		frustum.right,
+		frustum.bottom,
+		frustum.top,
+		frustum.near,
+		frustum.far or false
+	}
+
+	-- Skip the last test for infinite projections, it'll never fail.
+	if not planes[6] then
+		table.remove(planes)
+	end
+
+	for i = 1, #planes do
+		-- This is the current plane
+		local p = planes[i]
+
+		-- p-vertex selection (with the index trick)
+		-- According to the plane normal we can know the
+		-- indices of the positive vertex
+		local px = p.a > 0.0 and 2 or 1
+		local py = p.b > 0.0 and 2 or 1
+		local pz = p.c > 0.0 and 2 or 1
+
+		-- project p-vertex on plane normal
+		-- (How far is p-vertex from the origin)
+		local dot = (p.a * box[px].x) + (p.b * box[py].y) + (p.c * box[pz].z)
+
+		-- Doesn't intersect if it is behind the plane
+		if dot < -p.d then
+			return false
+		end
+	end
+
+	return true
+end
+
+-- outer.min is a vec3
+-- outer.max is a vec3
+-- inner.min is a vec3
+-- inner.max is a vec3
+function intersect.encapsulate_aabb(outer, inner)
+	return
+		outer.min.x <= inner.min.x and
+		outer.max.x >= inner.max.x and
+		outer.min.y <= inner.min.y and
+		outer.max.y >= inner.max.y and
+		outer.min.z <= inner.min.z and
+		outer.max.z >= inner.max.z
+end
+
+-- a.position is a vec3
+-- a.radius   is a number
+-- b.position is a vec3
+-- b.radius   is a number
+function intersect.circle_circle(a, b)
+	return a.position:dist(b.position) <= a.radius + b.radius
+end
+
+-- a.position is a vec3
+-- a.radius   is a number
+-- b.position is a vec3
+-- b.radius   is a number
+function intersect.sphere_sphere(a, b)
+	return intersect.circle_circle(a, b)
+end
+
+-- http://realtimecollisiondetection.net/blog/?p=103
+-- sphere.position is a vec3
+-- sphere.radius   is a number
+-- triangle[1]     is a vec3
+-- triangle[2]     is a vec3
+-- triangle[3]     is a vec3
+function intersect.sphere_triangle(sphere, triangle)
+	-- Sphere is centered at origin
+	local A  = triangle[1] - sphere.position
+	local B  = triangle[2] - sphere.position
+	local C  = triangle[3] - sphere.position
+
+	-- Compute normal of triangle plane
+	local V  = (B - A):cross(C - A)
+
+	-- Test if sphere lies outside triangle plane
+	local rr = sphere.radius * sphere.radius
+	local d  = A:dot(V)
+	local e  = V:dot(V)
+	local s1 = d * d > rr * e
+
+	-- Test if sphere lies outside triangle vertices
+	local aa = A:dot(A)
+	local ab = A:dot(B)
+	local ac = A:dot(C)
+	local bb = B:dot(B)
+	local bc = B:dot(C)
+	local cc = C:dot(C)
+
+	local s2 = (aa > rr) and (ab > aa) and (ac > aa)
+	local s3 = (bb > rr) and (ab > bb) and (bc > bb)
+	local s4 = (cc > rr) and (ac > cc) and (bc > cc)
+
+	-- Test is sphere lies outside triangle edges
+	local AB = B - A
+	local BC = C - B
+	local CA = A - C
+
+	local d1 = ab - aa
+	local d2 = bc - bb
+	local d3 = ac - cc
+
+	local e1 = AB:dot(AB)
+	local e2 = BC:dot(BC)
+	local e3 = CA:dot(CA)
+
+	local Q1 = A * e1 - AB * d1
+	local Q2 = B * e2 - BC * d2
+	local Q3 = C * e3 - CA * d3
+
+	local QC = C * e1 - Q1
+	local QA = A * e2 - Q2
+	local QB = B * e3 - Q3
+
+	local s5 = (Q1:dot(Q1) > rr * e1 * e1) and (Q1:dot(QC) > 0)
+	local s6 = (Q2:dot(Q2) > rr * e2 * e2) and (Q2:dot(QA) > 0)
+	local s7 = (Q3:dot(Q3) > rr * e3 * e3) and (Q3:dot(QB) > 0)
+
+	-- Return whether or not any of the tests passed
+	return s1 or s2 or s3 or s4 or s5 or s6 or s7
+end
+
+-- sphere.position is a vec3
+-- sphere.radius   is a number
+-- frustum.left    is a plane { a, b, c, d }
+-- frustum.right   is a plane { a, b, c, d }
+-- frustum.bottom  is a plane { a, b, c, d }
+-- frustum.top     is a plane { a, b, c, d }
+-- frustum.near    is a plane { a, b, c, d }
+-- frustum.far     is a plane { a, b, c, d }
+function intersect.sphere_frustum(sphere, frustum)
+	local x, y, z = sphere.position:unpack()
+	local planes  = {
+		frustum.left,
+		frustum.right,
+		frustum.bottom,
+		frustum.top,
+		frustum.near
+	}
+
+	if frustum.far then
+		table.insert(planes, frustum.far, 5)
+	end
+
+	local dot
+	for i = 1, #planes do
+		dot = planes[i].a * x + planes[i].b * y + planes[i].c * z + planes[i].d
+
+		if dot <= -sphere.radius then
+			return false
+		end
+	end
+
+	-- dot + radius is the distance of the object from the near plane.
+	-- make sure that the near plane is the last test!
+	return dot + sphere.radius
+end
+
+function intersect.capsule_capsule(c1, c2)
+	local dist2, p1, p2 = intersect.closest_point_segment_segment(c1.a, c1.b, c2.a, c2.b)
+	local radius = c1.radius + c2.radius
+
+	if dist2 <= radius * radius then
+		return p1, p2
+	end
+
+	return false
+end
+
+function intersect.closest_point_segment_segment(p1, p2, p3, p4)
+	local s  -- Distance of intersection along segment 1
+	local t  -- Distance of intersection along segment 2
+	local c1 -- Collision point on segment 1
+	local c2 -- Collision point on segment 2
+
+	local d1 = p2 - p1 -- Direction of segment 1
+	local d2 = p4 - p3 -- Direction of segment 2
+	local r  = p1 - p3
+	local a  = d1:dot(d1)
+	local e  = d2:dot(d2)
+	local f  = d2:dot(r)
+
+	-- Check if both segments degenerate into points
+	if a <= DBL_EPSILON and e <= DBL_EPSILON then
+		s  = 0
+		t  = 0
+		c1 = p1
+		c2 = p3
+		return (c1 - c2):dot(c1 - c2), s, t, c1, c2
+	end
+
+	-- Check if segment 1 degenerates into a point
+	if a <= DBL_EPSILON then
+		s = 0
+		t = utils.clamp(f / e, 0, 1)
+	else
+		local c = d1:dot(r)
+
+		-- Check is segment 2 degenerates into a point
+		if e <= DBL_EPSILON then
+			t = 0
+			s = utils.clamp(-c / a, 0, 1)
+		else
+			local b     = d1:dot(d2)
+			local denom = a * e - b * b
+
+			if abs(denom) > 0 then
+				s = utils.clamp((b * f - c * e) / denom, 0, 1)
+			else
+				s = 0
+			end
+
+			t = (b * s + f) / e
+
+			if t < 0 then
+				t = 0
+				s = utils.clamp(-c / a, 0, 1)
+			elseif t > 1 then
+				t = 1
+				s = utils.clamp((b - c) / a, 0, 1)
+			end
+		end
+	end
+
+	c1 = p1 + d1 * s
+	c2 = p3 + d2 * t
+
+	return (c1 - c2):dot(c1 - c2), c1, c2, s, t
+end
+
+return intersect

libs/cpml/mat4.lua

@@ -0,0 +1,943 @@
+--- double 4x4, 1-based, column major matrices
+-- @module mat4
+local modules   = (...):gsub('%.[^%.]+$', '') .. "."
+local constants = require(modules .. "constants")
+local vec2      = require(modules .. "vec2")
+local vec3      = require(modules .. "vec3")
+local quat      = require(modules .. "quat")
+local utils     = require(modules .. "utils")
+local precond   = require(modules .. "_private_precond")
+local private   = require(modules .. "_private_utils")
+local sqrt      = math.sqrt
+local cos       = math.cos
+local sin       = math.sin
+local tan       = math.tan
+local rad       = math.rad
+local mat4      = {}
+local mat4_mt   = {}
+
+-- Private constructor.
+local function new(m)
+	m = m or {
+		0, 0, 0, 0,
+		0, 0, 0, 0,
+		0, 0, 0, 0,
+		0, 0, 0, 0
+	}
+	m._m = m
+	return setmetatable(m, mat4_mt)
+end
+
+-- Convert matrix into identity
+local function identity(m)
+	m[1],  m[2],  m[3],  m[4]  = 1, 0, 0, 0
+	m[5],  m[6],  m[7],  m[8]  = 0, 1, 0, 0
+	m[9],  m[10], m[11], m[12] = 0, 0, 1, 0
+	m[13], m[14], m[15], m[16] = 0, 0, 0, 1
+	return m
+end
+
+-- Do the check to see if JIT is enabled. If so use the optimized FFI structs.
+local status, ffi
+if type(jit) == "table" and jit.status() then
+	--  status, ffi = pcall(require, "ffi")
+	if status then
+		ffi.cdef "typedef struct { double _m[16]; } cpml_mat4;"
+		new = ffi.typeof("cpml_mat4")
+	end
+end
+
+-- Statically allocate a temporary variable used in some of our functions.
+local tmp = new()
+local tm4 = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
+local tv4 = { 0, 0, 0, 0 }
+
+--- The public constructor.
+-- @param a Can be of four types: </br>
+-- table Length 16 (4x4 matrix)
+-- table Length 9 (3x3 matrix)
+-- table Length 4 (4 vec4s)
+-- nil
+-- @treturn mat4 out
+function mat4.new(a)
+	local out = new()
+
+	-- 4x4 matrix
+	if type(a) == "table" and #a == 16 then
+		for i = 1, 16 do
+			out[i] = tonumber(a[i])
+		end
+
+	-- 3x3 matrix
+	elseif type(a) == "table" and #a == 9 then
+		out[1], out[2],  out[3]  = a[1], a[2], a[3]
+		out[5], out[6],  out[7]  = a[4], a[5], a[6]
+		out[9], out[10], out[11] = a[7], a[8], a[9]
+		out[16] = 1
+
+	-- 4 vec4s
+	elseif type(a) == "table" and type(a[1]) == "table" then
+		local idx = 1
+		for i = 1, 4 do
+			for j = 1, 4 do
+				out[idx] = a[i][j]
+				idx = idx + 1
+			end
+		end
+
+	-- nil
+	else
+		out[1]  = 1
+		out[6]  = 1
+		out[11] = 1
+		out[16] = 1
+	end
+
+	return out
+end
+
+--- Create an identity matrix.
+-- @tparam mat4 a Matrix to overwrite
+-- @treturn mat4 out
+function mat4.identity(a)
+	return identity(a or new())
+end
+
+--- Create a matrix from an angle/axis pair.
+-- @tparam number angle Angle of rotation
+-- @tparam vec3 axis Axis of rotation
+-- @treturn mat4 out
+function mat4.from_angle_axis(angle, axis)
+	local l = axis:len()
+	if l == 0 then
+		return new()
+	end
+
+	local x, y, z = axis.x / l, axis.y / l, axis.z / l
+	local c = cos(angle)
+	local s = sin(angle)
+
+	return new {
+		x*x*(1-c)+c,   y*x*(1-c)+z*s, x*z*(1-c)-y*s, 0,
+		x*y*(1-c)-z*s, y*y*(1-c)+c,   y*z*(1-c)+x*s, 0,
+		x*z*(1-c)+y*s, y*z*(1-c)-x*s, z*z*(1-c)+c,   0,
+		0, 0, 0, 1
+	}
+end
+
+--- Create a matrix from a quaternion.
+-- @tparam quat q Rotation quaternion
+-- @treturn mat4 out
+function mat4.from_quaternion(q)
+	return mat4.from_angle_axis(q:to_angle_axis())
+end
+
+--- Create a matrix from a direction/up pair.
+-- @tparam vec3 direction Vector direction
+-- @tparam vec3 up Up direction
+-- @treturn mat4 out
+function mat4.from_direction(direction, up)
+	local forward = vec3.normalize(direction)
+	local side = vec3.cross(forward, up):normalize()
+	local new_up = vec3.cross(side, forward):normalize()
+
+	local out = new()
+	out[1]    = side.x
+	out[5]    = side.y
+	out[9]    = side.z
+	out[2]    = new_up.x
+	out[6]    = new_up.y
+	out[10]   = new_up.z
+	out[3]    = forward.x
+	out[7]    = forward.y
+	out[11]   = forward.z
+	out[16]   = 1
+
+	return out
+end
+
+--- Create a matrix from a transform.
+-- @tparam vec3 trans Translation vector
+-- @tparam quat rot Rotation quaternion
+-- @tparam vec3 scale Scale vector
+-- @treturn mat4 out
+function mat4.from_transform(trans, rot, scale)
+	local rx, ry, rz, rw = rot.x, rot.y, rot.z, rot.w
+
+	local sm = new {
+		scale.x, 0,       0,       0,
+		0,       scale.y, 0,       0,
+		0,       0,       scale.z, 0,
+		0,       0,       0,       1,
+	}
+
+	local rm = new {
+		1-2*(ry*ry+rz*rz), 2*(rx*ry-rz*rw), 2*(rx*rz+ry*rw), 0,
+		2*(rx*ry+rz*rw), 1-2*(rx*rx+rz*rz), 2*(ry*rz-rx*rw), 0,
+		2*(rx*rz-ry*rw), 2*(ry*rz+rx*rw), 1-2*(rx*rx+ry*ry), 0,
+		0, 0, 0, 1
+	}
+
+	local rsm = rm * sm
+
+	rsm[13] = trans.x
+	rsm[14] = trans.y
+	rsm[15] = trans.z
+
+	return rsm
+end
+
+--- Create matrix from orthogonal.
+-- @tparam number left
+-- @tparam number right
+-- @tparam number top
+-- @tparam number bottom
+-- @tparam number near
+-- @tparam number far
+-- @treturn mat4 out
+function mat4.from_ortho(left, right, top, bottom, near, far)
+	local out = new()
+	out[1]    =  2 / (right - left)
+	out[6]    =  2 / (top - bottom)
+	out[11]   = -2 / (far - near)
+	out[13]   = -((right + left) / (right - left))
+	out[14]   = -((top + bottom) / (top - bottom))
+	out[15]   = -((far + near) / (far - near))
+	out[16]   =  1
+
+	return out
+end
+
+--- Create matrix from perspective.
+-- @tparam number fovy Field of view
+-- @tparam number aspect Aspect ratio
+-- @tparam number near Near plane
+-- @tparam number far Far plane
+-- @treturn mat4 out
+function mat4.from_perspective(fovy, aspect, near, far)
+	assert(aspect ~= 0)
+	assert(near   ~= far)
+
+	local t   = tan(rad(fovy) / 2)
+	local out = new()
+	out[1]    =  1 / (t * aspect)
+	out[6]    =  1 / t
+	out[11]   = -(far + near) / (far - near)
+	out[12]   = -1
+	out[15]   = -(2 * far * near) / (far - near)
+	out[16]   =  0
+
+	return out
+end
+
+-- Adapted from the Oculus SDK.
+--- Create matrix from HMD perspective.
+-- @tparam number tanHalfFov Tangent of half of the field of view
+-- @tparam number zNear Near plane
+-- @tparam number zFar Far plane
+-- @tparam boolean flipZ Z axis is flipped or not
+-- @tparam boolean farAtInfinity Far plane is infinite or not
+-- @treturn mat4 out
+function mat4.from_hmd_perspective(tanHalfFov, zNear, zFar, flipZ, farAtInfinity)
+	-- CPML is right-handed and intended for GL, so these don't need to be arguments.
+	local rightHanded = true
+	local isOpenGL    = true
+
+	local function CreateNDCScaleAndOffsetFromFov(tanHalfFov)
+		local x_scale  = 2 / (tanHalfFov.LeftTan + tanHalfFov.RightTan)
+		local x_offset =     (tanHalfFov.LeftTan - tanHalfFov.RightTan) * x_scale * 0.5
+		local y_scale  = 2 / (tanHalfFov.UpTan   + tanHalfFov.DownTan )
+		local y_offset =     (tanHalfFov.UpTan   - tanHalfFov.DownTan ) * y_scale * 0.5
+
+		local result = {
+			Scale  = vec2(x_scale, y_scale),
+			Offset = vec2(x_offset, y_offset)
+		}
+
+		-- Hey - why is that Y.Offset negated?
+		-- It's because a projection matrix transforms from world coords with Y=up,
+		-- whereas this is from NDC which is Y=down.
+		 return result
+	end
+
+	if not flipZ and farAtInfinity then
+		print("Error: Cannot push Far Clip to Infinity when Z-order is not flipped")
+		farAtInfinity = false
+	end
+
+	 -- A projection matrix is very like a scaling from NDC, so we can start with that.
+	local scaleAndOffset  = CreateNDCScaleAndOffsetFromFov(tanHalfFov)
+	local handednessScale = rightHanded and -1.0 or 1.0
+	local projection      = new()
+
+	-- Produces X result, mapping clip edges to [-w,+w]
+	projection[1] = scaleAndOffset.Scale.x
+	projection[2] = 0
+	projection[3] = handednessScale * scaleAndOffset.Offset.x
+	projection[4] = 0
+
+	-- Produces Y result, mapping clip edges to [-w,+w]
+	-- Hey - why is that YOffset negated?
+	-- It's because a projection matrix transforms from world coords with Y=up,
+	-- whereas this is derived from an NDC scaling, which is Y=down.
+	projection[5] = 0
+	projection[6] = scaleAndOffset.Scale.y
+	projection[7] = handednessScale * -scaleAndOffset.Offset.y
+	projection[8] = 0
+
+	-- Produces Z-buffer result - app needs to fill this in with whatever Z range it wants.
+	-- We'll just use some defaults for now.
+	projection[9]  = 0
+	projection[10] = 0
+
+	if farAtInfinity then
+		if isOpenGL then
+			-- It's not clear this makes sense for OpenGL - you don't get the same precision benefits you do in D3D.
+			projection[11] = -handednessScale
+			projection[12] = 2.0 * zNear
+		else
+			projection[11] = 0
+			projection[12] = zNear
+		end
+	else
+		if isOpenGL then
+			-- Clip range is [-w,+w], so 0 is at the middle of the range.
+			projection[11] = -handednessScale * (flipZ and -1.0 or 1.0) * (zNear + zFar) / (zNear - zFar)
+			projection[12] = 2.0 * ((flipZ and -zFar or zFar) * zNear) / (zNear - zFar)
+		else
+			-- Clip range is [0,+w], so 0 is at the start of the range.
+			projection[11] = -handednessScale * (flipZ and -zNear or zFar) / (zNear - zFar)
+			projection[12] = ((flipZ and -zFar or zFar) * zNear) / (zNear - zFar)
+		end
+	end
+
+	-- Produces W result (= Z in)
+	projection[13] = 0
+	projection[14] = 0
+	projection[15] = handednessScale
+	projection[16] = 0
+
+	return projection:transpose(projection)
+end
+
+--- Clone a matrix.
+-- @tparam mat4 a Matrix to clone
+-- @treturn mat4 out
+function mat4.clone(a)
+	return new(a)
+end
+
+function mul_internal(out, a, b)
+	tm4[1]  = b[1]  * a[1] + b[2]  * a[5] + b[3]  * a[9]  + b[4]  * a[13]
+	tm4[2]  = b[1]  * a[2] + b[2]  * a[6] + b[3]  * a[10] + b[4]  * a[14]
+	tm4[3]  = b[1]  * a[3] + b[2]  * a[7] + b[3]  * a[11] + b[4]  * a[15]
+	tm4[4]  = b[1]  * a[4] + b[2]  * a[8] + b[3]  * a[12] + b[4]  * a[16]
+	tm4[5]  = b[5]  * a[1] + b[6]  * a[5] + b[7]  * a[9]  + b[8]  * a[13]
+	tm4[6]  = b[5]  * a[2] + b[6]  * a[6] + b[7]  * a[10] + b[8]  * a[14]
+	tm4[7]  = b[5]  * a[3] + b[6]  * a[7] + b[7]  * a[11] + b[8]  * a[15]
+	tm4[8]  = b[5]  * a[4] + b[6]  * a[8] + b[7]  * a[12] + b[8]  * a[16]
+	tm4[9]  = b[9]  * a[1] + b[10] * a[5] + b[11] * a[9]  + b[12] * a[13]
+	tm4[10] = b[9]  * a[2] + b[10] * a[6] + b[11] * a[10] + b[12] * a[14]
+	tm4[11] = b[9]  * a[3] + b[10] * a[7] + b[11] * a[11] + b[12] * a[15]
+	tm4[12] = b[9]  * a[4] + b[10] * a[8] + b[11] * a[12] + b[12] * a[16]
+	tm4[13] = b[13] * a[1] + b[14] * a[5] + b[15] * a[9]  + b[16] * a[13]
+	tm4[14] = b[13] * a[2] + b[14] * a[6] + b[15] * a[10] + b[16] * a[14]
+	tm4[15] = b[13] * a[3] + b[14] * a[7] + b[15] * a[11] + b[16] * a[15]
+	tm4[16] = b[13] * a[4] + b[14] * a[8] + b[15] * a[12] + b[16] * a[16]
+
+	for i = 1, 16 do
+		out[i] = tm4[i]
+	end
+end
+
+--- Multiply N matrices.
+-- @tparam mat4 out Matrix to store the result
+-- @tparam mat4 or {mat4, ...} left hand operand(s)
+-- @tparam mat4 right hand operand if a is not table
+-- @treturn mat4 out multiplied matrix result
+function mat4.mul(out, a, b)
+	if mat4.is_mat4(a) then
+		mul_internal(out, a, b)
+		return out
+	end
+	if #a == 0 then
+		identity(out)
+	elseif #a == 1 then
+		-- only one matrix, just copy
+		for i = 1, 16 do
+			out[i] = a[1][i]
+		end
+	else
+		local ma = a[1]
+		local mb = a[2]
+		for i = 2, #a do
+			mul_internal(out, ma, mb)
+			ma = out
+		end
+	end
+	return out
+end
+
+--- Multiply a matrix and a vec3, with perspective division.
+-- This function uses an implicit 1 for the fourth component.
+-- @tparam vec3 out vec3 to store the result
+-- @tparam mat4 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn vec3 out
+function mat4.mul_vec3_perspective(out, a, b)
+	local v4x = b.x * a[1] + b.y * a[5] + b.z * a[9]  + a[13]
+	local v4y = b.x * a[2] + b.y * a[6] + b.z * a[10] + a[14]
+	local v4z = b.x * a[3] + b.y * a[7] + b.z * a[11] + a[15]
+	local v4w = b.x * a[4] + b.y * a[8] + b.z * a[12] + a[16]
+	local inv_w = 0
+	if v4w ~= 0 then
+		inv_w = utils.sign(v4w) / v4w
+	end
+	out.x = v4x * inv_w
+	out.y = v4y * inv_w
+	out.z = v4z * inv_w
+	return out
+end
+
+--- Multiply a matrix and a vec4.
+-- @tparam table out table to store the result
+-- @tparam mat4 a Left hand operand
+-- @tparam table b Right hand operand
+-- @treturn vec4 out
+function mat4.mul_vec4(out, a, b)
+	tv4[1] = b[1] * a[1] + b[2] * a[5] + b [3] * a[9]  + b[4] * a[13]
+	tv4[2] = b[1] * a[2] + b[2] * a[6] + b [3] * a[10] + b[4] * a[14]
+	tv4[3] = b[1] * a[3] + b[2] * a[7] + b [3] * a[11] + b[4] * a[15]
+	tv4[4] = b[1] * a[4] + b[2] * a[8] + b [3] * a[12] + b[4] * a[16]
+
+	for i = 1, 4 do
+		out[i] = tv4[i]
+	end
+
+	return out
+end
+
+--- Invert a matrix.
+-- @tparam mat4 out Matrix to store the result
+-- @tparam mat4 a Matrix to invert
+-- @treturn mat4 out
+function mat4.invert(out, a)
+	tm4[1]  =  a[6] * a[11] * a[16] - a[6] * a[12] * a[15] - a[10] * a[7] * a[16] + a[10] * a[8] * a[15] + a[14] * a[7] * a[12] - a[14] * a[8] * a[11]
+	tm4[2]  = -a[2] * a[11] * a[16] + a[2] * a[12] * a[15] + a[10] * a[3] * a[16] - a[10] * a[4] * a[15] - a[14] * a[3] * a[12] + a[14] * a[4] * a[11]
+	tm4[3]  =  a[2] * a[7]  * a[16] - a[2] * a[8]  * a[15] - a[6]  * a[3] * a[16] + a[6]  * a[4] * a[15] + a[14] * a[3] * a[8]  - a[14] * a[4] * a[7]
+	tm4[4]  = -a[2] * a[7]  * a[12] + a[2] * a[8]  * a[11] + a[6]  * a[3] * a[12] - a[6]  * a[4] * a[11] - a[10] * a[3] * a[8]  + a[10] * a[4] * a[7]
+	tm4[5]  = -a[5] * a[11] * a[16] + a[5] * a[12] * a[15] + a[9]  * a[7] * a[16] - a[9]  * a[8] * a[15] - a[13] * a[7] * a[12] + a[13] * a[8] * a[11]
+	tm4[6]  =  a[1] * a[11] * a[16] - a[1] * a[12] * a[15] - a[9]  * a[3] * a[16] + a[9]  * a[4] * a[15] + a[13] * a[3] * a[12] - a[13] * a[4] * a[11]
+	tm4[7]  = -a[1] * a[7]  * a[16] + a[1] * a[8]  * a[15] + a[5]  * a[3] * a[16] - a[5]  * a[4] * a[15] - a[13] * a[3] * a[8]  + a[13] * a[4] * a[7]
+	tm4[8]  =  a[1] * a[7]  * a[12] - a[1] * a[8]  * a[11] - a[5]  * a[3] * a[12] + a[5]  * a[4] * a[11] + a[9]  * a[3] * a[8]  - a[9]  * a[4] * a[7]
+	tm4[9]  =  a[5] * a[10] * a[16] - a[5] * a[12] * a[14] - a[9]  * a[6] * a[16] + a[9]  * a[8] * a[14] + a[13] * a[6] * a[12] - a[13] * a[8] * a[10]
+	tm4[10] = -a[1] * a[10] * a[16] + a[1] * a[12] * a[14] + a[9]  * a[2] * a[16] - a[9]  * a[4] * a[14] - a[13] * a[2] * a[12] + a[13] * a[4] * a[10]
+	tm4[11] =  a[1] * a[6]  * a[16] - a[1] * a[8]  * a[14] - a[5]  * a[2] * a[16] + a[5]  * a[4] * a[14] + a[13] * a[2] * a[8]  - a[13] * a[4] * a[6]
+	tm4[12] = -a[1] * a[6]  * a[12] + a[1] * a[8]  * a[10] + a[5]  * a[2] * a[12] - a[5]  * a[4] * a[10] - a[9]  * a[2] * a[8]  + a[9]  * a[4] * a[6]
+	tm4[13] = -a[5] * a[10] * a[15] + a[5] * a[11] * a[14] + a[9]  * a[6] * a[15] - a[9]  * a[7] * a[14] - a[13] * a[6] * a[11] + a[13] * a[7] * a[10]
+	tm4[14] =  a[1] * a[10] * a[15] - a[1] * a[11] * a[14] - a[9]  * a[2] * a[15] + a[9]  * a[3] * a[14] + a[13] * a[2] * a[11] - a[13] * a[3] * a[10]
+	tm4[15] = -a[1] * a[6]  * a[15] + a[1] * a[7]  * a[14] + a[5]  * a[2] * a[15] - a[5]  * a[3] * a[14] - a[13] * a[2] * a[7]  + a[13] * a[3] * a[6]
+	tm4[16] =  a[1] * a[6]  * a[11] - a[1] * a[7]  * a[10] - a[5]  * a[2] * a[11] + a[5]  * a[3] * a[10] + a[9]  * a[2] * a[7]  - a[9]  * a[3] * a[6]
+
+	local det = a[1] * tm4[1] + a[2] * tm4[5] + a[3] * tm4[9] + a[4] * tm4[13]
+
+	if det == 0 then return a end
+
+	det = 1 / det
+
+	for i = 1, 16 do
+		out[i] = tm4[i] * det
+	end
+
+	return out
+end
+
+--- Scale a matrix.
+-- @tparam mat4 out Matrix to store the result
+-- @tparam mat4 a Matrix to scale
+-- @tparam vec3 s Scalar
+-- @treturn mat4 out
+function mat4.scale(out, a, s)
+	identity(tmp)
+	tmp[1]  = s.x
+	tmp[6]  = s.y
+	tmp[11] = s.z
+
+	return out:mul(tmp, a)
+end
+
+--- Rotate a matrix.
+-- @tparam mat4 out Matrix to store the result
+-- @tparam mat4 a Matrix to rotate
+-- @tparam number angle Angle to rotate by (in radians)
+-- @tparam vec3 axis Axis to rotate on
+-- @treturn mat4 out
+function mat4.rotate(out, a, angle, axis)
+	if type(angle) == "table" or type(angle) == "cdata" then
+		angle, axis = angle:to_angle_axis()
+	end
+
+	local l = axis:len()
+
+	if l == 0 then
+		return a
+	end
+
+	local x, y, z = axis.x / l, axis.y / l, axis.z / l
+	local c = cos(angle)
+	local s = sin(angle)
+
+	identity(tmp)
+	tmp[1]  = x * x * (1 - c) + c
+	tmp[2]  = y * x * (1 - c) + z * s
+	tmp[3]  = x * z * (1 - c) - y * s
+	tmp[5]  = x * y * (1 - c) - z * s
+	tmp[6]  = y * y * (1 - c) + c
+	tmp[7]  = y * z * (1 - c) + x * s
+	tmp[9]  = x * z * (1 - c) + y * s
+	tmp[10] = y * z * (1 - c) - x * s
+	tmp[11] = z * z * (1 - c) + c
+
+	return out:mul(tmp, a)
+end
+
+--- Translate a matrix.
+-- @tparam mat4 out Matrix to store the result
+-- @tparam mat4 a Matrix to translate
+-- @tparam vec3 t Translation vector
+-- @treturn mat4 out
+function mat4.translate(out, a, t)
+	identity(tmp)
+	tmp[13] = t.x
+	tmp[14] = t.y
+	tmp[15] = t.z
+
+	return out:mul(tmp, a)
+end
+
+--- Shear a matrix.
+-- @tparam mat4 out Matrix to store the result
+-- @tparam mat4 a Matrix to translate
+-- @tparam number yx
+-- @tparam number zx
+-- @tparam number xy
+-- @tparam number zy
+-- @tparam number xz
+-- @tparam number yz
+-- @treturn mat4 out
+function mat4.shear(out, a, yx, zx, xy, zy, xz, yz)
+	identity(tmp)
+	tmp[2]  = yx or 0
+	tmp[3]  = zx or 0
+	tmp[5]  = xy or 0
+	tmp[7]  = zy or 0
+	tmp[9]  = xz or 0
+	tmp[10] = yz or 0
+
+	return out:mul(tmp, a)
+end
+
+--- Reflect a matrix across a plane.
+-- @tparam mat4 Matrix to store the result
+-- @tparam a Matrix to reflect
+-- @tparam vec3 position A point on the plane
+-- @tparam vec3 normal The (normalized!) normal vector of the plane
+function mat4.reflect(out, a, position, normal)
+	local nx, ny, nz = normal:unpack()
+	local d = -position:dot(normal)
+	tmp[1] = 1 - 2 * nx ^ 2
+	tmp[2] = 2 * nx * ny
+	tmp[3] = -2 * nx * nz
+	tmp[4] = 0
+	tmp[5] = -2 * nx * ny
+	tmp[6] = 1 - 2 * ny ^ 2
+	tmp[7] = -2 * ny * nz
+	tmp[8] = 0
+	tmp[9] = -2 * nx * nz
+	tmp[10] = -2 * ny * nz
+	tmp[11] = 1 - 2 * nz ^ 2
+	tmp[12] = 0
+	tmp[13] = -2 * nx * d
+	tmp[14] = -2 * ny * d
+	tmp[15] = -2 * nz * d
+	tmp[16] = 1
+
+	return out:mul(tmp, a)
+end
+
+--- Transform matrix to look at a point.
+-- @tparam mat4 out Matrix to store result
+-- @tparam vec3 eye Location of viewer's view plane
+-- @tparam vec3 center Location of object to view
+-- @tparam vec3 up Up direction
+-- @treturn mat4 out
+function mat4.look_at(out, eye, look_at, up)
+	local z_axis = (eye - look_at):normalize()
+	local x_axis = up:cross(z_axis):normalize()
+	local y_axis = z_axis:cross(x_axis)
+	out[1] = x_axis.x
+	out[2] = y_axis.x
+	out[3] = z_axis.x
+	out[4] = 0
+	out[5] = x_axis.y
+	out[6] = y_axis.y
+	out[7] = z_axis.y
+	out[8] = 0
+	out[9] = x_axis.z
+	out[10] = y_axis.z
+	out[11] = z_axis.z
+	out[12] = 0
+	out[13] = -out[  1]*eye.x - out[4+1]*eye.y - out[8+1]*eye.z
+	out[14] = -out[  2]*eye.x - out[4+2]*eye.y - out[8+2]*eye.z
+	out[15] = -out[  3]*eye.x - out[4+3]*eye.y - out[8+3]*eye.z
+	out[16] = -out[  4]*eye.x - out[4+4]*eye.y - out[8+4]*eye.z + 1
+	return out
+end
+
+--- Transform matrix to target a point.
+-- @tparam mat4 out Matrix to store result
+-- @tparam vec3 eye Location of viewer's view plane
+-- @tparam vec3 center Location of object to view
+-- @tparam vec3 up Up direction
+-- @treturn mat4 out
+function mat4.target(out, from, to, up)
+	local z_axis = (from - to):normalize()
+	local x_axis = up:cross(z_axis):normalize()
+	local y_axis = z_axis:cross(x_axis)
+	out[1] = x_axis.x
+	out[2] = x_axis.y
+	out[3] = x_axis.z
+	out[4] = 0
+	out[5] = y_axis.x
+	out[6] = y_axis.y
+	out[7] = y_axis.z
+	out[8] = 0
+	out[9] = z_axis.x
+	out[10] = z_axis.y
+	out[11] = z_axis.z
+	out[12] = 0
+	out[13] = from.x
+	out[14] = from.y
+	out[15] = from.z
+	out[16] = 1
+	return out
+end
+
+--- Transpose a matrix.
+-- @tparam mat4 out Matrix to store the result
+-- @tparam mat4 a Matrix to transpose
+-- @treturn mat4 out
+function mat4.transpose(out, a)
+	tm4[1]  = a[1]
+	tm4[2]  = a[5]
+	tm4[3]  = a[9]
+	tm4[4]  = a[13]
+	tm4[5]  = a[2]
+	tm4[6]  = a[6]
+	tm4[7]  = a[10]
+	tm4[8]  = a[14]
+	tm4[9]  = a[3]
+	tm4[10] = a[7]
+	tm4[11] = a[11]
+	tm4[12] = a[15]
+	tm4[13] = a[4]
+	tm4[14] = a[8]
+	tm4[15] = a[12]
+	tm4[16] = a[16]
+
+	for i = 1, 16 do
+		out[i] = tm4[i]
+	end
+
+	return out
+end
+
+--- Project a point into screen space
+-- @tparam vec3 obj Object position in world space
+-- @tparam mat4 mvp Projection matrix
+-- @tparam table viewport XYWH of viewport
+-- @treturn vec3 win
+function mat4.project(obj, mvp, viewport)
+	local point = mat4.mul_vec3_perspective(vec3(), mvp, obj)
+	point.x = point.x * 0.5 + 0.5
+	point.y = point.y * 0.5 + 0.5
+	point.z = point.z * 0.5 + 0.5
+	point.x = point.x * viewport[3] + viewport[1]
+	point.y = point.y * viewport[4] + viewport[2]
+	return point
+end
+
+--- Unproject a point from screen space to world space.
+-- @tparam vec3 win Object position in screen space
+-- @tparam mat4 mvp Projection matrix
+-- @tparam table viewport XYWH of viewport
+-- @treturn vec3 obj
+function mat4.unproject(win, mvp, viewport)
+	local point = vec3.clone(win)
+
+	-- 0..n -> 0..1
+	point.x = (point.x - viewport[1]) / viewport[3]
+	point.y = (point.y - viewport[2]) / viewport[4]
+
+	-- 0..1 -> -1..1
+	point.x = point.x * 2 - 1
+	point.y = point.y * 2 - 1
+	point.z = point.z * 2 - 1
+
+	return mat4.mul_vec3_perspective(point, tmp:invert(mvp), point)
+end
+
+--- Return a boolean showing if a table is or is not a mat4.
+-- @tparam mat4 a Matrix to be tested
+-- @treturn boolean is_mat4
+function mat4.is_mat4(a)
+	if type(a) == "cdata" then
+		return ffi.istype("cpml_mat4", a)
+	end
+
+	if type(a) ~= "table" then
+		return false
+	end
+
+	for i = 1, 16 do
+		if type(a[i]) ~= "number" then
+			return false
+		end
+	end
+
+	return true
+end
+
+--- Return whether any component is NaN
+-- @tparam mat4 a Matrix to be tested
+-- @treturn boolean if any component is NaN
+function vec2.has_nan(a)
+	for i = 1, 16 do
+		if private.is_nan(a[i]) then
+			return true
+		end
+	end
+	return false
+end
+
+--- Return a formatted string.
+-- @tparam mat4 a Matrix to be turned into a string
+-- @treturn string formatted
+function mat4.to_string(a)
+	local str = "[ "
+	for i = 1, 16 do
+		str = str .. string.format("%+0.3f", a[i])
+		if i < 16 then
+			str = str .. ", "
+		end
+	end
+	str = str .. " ]"
+	return str
+end
+
+--- Convert a matrix to row vec4s.
+-- @tparam mat4 a Matrix to be converted
+-- @treturn table vec4s
+function mat4.to_vec4s(a)
+	return {
+		{ a[1],  a[2],  a[3],  a[4]  },
+		{ a[5],  a[6],  a[7],  a[8]  },
+		{ a[9],  a[10], a[11], a[12] },
+		{ a[13], a[14], a[15], a[16] }
+	}
+end
+
+--- Convert a matrix to col vec4s.
+-- @tparam mat4 a Matrix to be converted
+-- @treturn table vec4s
+function mat4.to_vec4s_cols(a)
+	return {
+		{ a[1], a[5], a[9],  a[13] },
+		{ a[2], a[6], a[10], a[14] },
+		{ a[3], a[7], a[11], a[15] },
+		{ a[4], a[8], a[12], a[16] }
+	}
+end
+
+-- http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/
+--- Convert a matrix to a quaternion.
+-- @tparam mat4 a Matrix to be converted
+-- @treturn quat out
+function mat4.to_quat(a)
+	identity(tmp):transpose(a)
+
+	local w     = sqrt(1 + tmp[1] + tmp[6] + tmp[11]) / 2
+	local scale = w * 4
+	local q     = quat.new(
+		tmp[10] - tmp[7] / scale,
+		tmp[3]  - tmp[9] / scale,
+		tmp[5]  - tmp[2] / scale,
+		w
+	)
+
+	return q:normalize(q)
+end
+
+-- http://www.crownandcutlass.com/features/technicaldetails/frustum.html
+--- Convert a matrix to a frustum.
+-- @tparam mat4 a Matrix to be converted (projection * view)
+-- @tparam boolean infinite Infinite removes the far plane
+-- @treturn frustum out
+function mat4.to_frustum(a, infinite)
+	local t
+	local frustum = {}
+
+	-- Extract the LEFT plane
+	frustum.left   = {}
+	frustum.left.a = a[4]  + a[1]
+	frustum.left.b = a[8]  + a[5]
+	frustum.left.c = a[12] + a[9]
+	frustum.left.d = a[16] + a[13]
+
+	-- Normalize the result
+	t = sqrt(frustum.left.a * frustum.left.a + frustum.left.b * frustum.left.b + frustum.left.c * frustum.left.c)
+	frustum.left.a = frustum.left.a / t
+	frustum.left.b = frustum.left.b / t
+	frustum.left.c = frustum.left.c / t
+	frustum.left.d = frustum.left.d / t
+
+	-- Extract the RIGHT plane
+	frustum.right   = {}
+	frustum.right.a = a[4]  - a[1]
+	frustum.right.b = a[8]  - a[5]
+	frustum.right.c = a[12] - a[9]
+	frustum.right.d = a[16] - a[13]
+
+	-- Normalize the result
+	t = sqrt(frustum.right.a * frustum.right.a + frustum.right.b * frustum.right.b + frustum.right.c * frustum.right.c)
+	frustum.right.a = frustum.right.a / t
+	frustum.right.b = frustum.right.b / t
+	frustum.right.c = frustum.right.c / t
+	frustum.right.d = frustum.right.d / t
+
+	-- Extract the BOTTOM plane
+	frustum.bottom   = {}
+	frustum.bottom.a = a[4]  + a[2]
+	frustum.bottom.b = a[8]  + a[6]
+	frustum.bottom.c = a[12] + a[10]
+	frustum.bottom.d = a[16] + a[14]
+
+	-- Normalize the result
+	t = sqrt(frustum.bottom.a * frustum.bottom.a + frustum.bottom.b * frustum.bottom.b + frustum.bottom.c * frustum.bottom.c)
+	frustum.bottom.a = frustum.bottom.a / t
+	frustum.bottom.b = frustum.bottom.b / t
+	frustum.bottom.c = frustum.bottom.c / t
+	frustum.bottom.d = frustum.bottom.d / t
+
+	-- Extract the TOP plane
+	frustum.top   = {}
+	frustum.top.a = a[4]  - a[2]
+	frustum.top.b = a[8]  - a[6]
+	frustum.top.c = a[12] - a[10]
+	frustum.top.d = a[16] - a[14]
+
+	-- Normalize the result
+	t = sqrt(frustum.top.a * frustum.top.a + frustum.top.b * frustum.top.b + frustum.top.c * frustum.top.c)
+	frustum.top.a = frustum.top.a / t
+	frustum.top.b = frustum.top.b / t
+	frustum.top.c = frustum.top.c / t
+	frustum.top.d = frustum.top.d / t
+
+	-- Extract the NEAR plane
+	frustum.near   = {}
+	frustum.near.a = a[4]  + a[3]
+	frustum.near.b = a[8]  + a[7]
+	frustum.near.c = a[12] + a[11]
+	frustum.near.d = a[16] + a[15]
+
+	-- Normalize the result
+	t = sqrt(frustum.near.a * frustum.near.a + frustum.near.b * frustum.near.b + frustum.near.c * frustum.near.c)
+	frustum.near.a = frustum.near.a / t
+	frustum.near.b = frustum.near.b / t
+	frustum.near.c = frustum.near.c / t
+	frustum.near.d = frustum.near.d / t
+
+	if not infinite then
+		-- Extract the FAR plane
+		frustum.far   = {}
+		frustum.far.a = a[4]  - a[3]
+		frustum.far.b = a[8]  - a[7]
+		frustum.far.c = a[12] - a[11]
+		frustum.far.d = a[16] - a[15]
+
+		-- Normalize the result
+		t = sqrt(frustum.far.a * frustum.far.a + frustum.far.b * frustum.far.b + frustum.far.c * frustum.far.c)
+		frustum.far.a = frustum.far.a / t
+		frustum.far.b = frustum.far.b / t
+		frustum.far.c = frustum.far.c / t
+		frustum.far.d = frustum.far.d / t
+	end
+
+	return frustum
+end
+
+function mat4_mt.__index(t, k)
+	if type(t) == "cdata" then
+		if type(k) == "number" then
+			return t._m[k-1]
+		end
+	end
+
+	return rawget(mat4, k)
+end
+
+function mat4_mt.__newindex(t, k, v)
+	if type(t) == "cdata" then
+		if type(k) == "number" then
+			t._m[k-1] = v
+		end
+	end
+end
+
+mat4_mt.__tostring = mat4.to_string
+
+function mat4_mt.__call(_, a)
+	return mat4.new(a)
+end
+
+function mat4_mt.__unm(a)
+	return new():invert(a)
+end
+
+function mat4_mt.__eq(a, b)
+	if not mat4.is_mat4(a) or not mat4.is_mat4(b) then
+		return false
+	end
+
+	for i = 1, 16 do
+		if not utils.tolerance(b[i]-a[i], constants.FLT_EPSILON) then
+			return false
+		end
+	end
+
+	return true
+end
+
+function mat4_mt.__mul(a, b)
+	precond.assert(mat4.is_mat4(a), "__mul: Wrong argument type '%s' for left hand operand. (<cpml.mat4> expected)", type(a))
+
+	if vec3.is_vec3(b) then
+		return mat4.mul_vec3_perspective(vec3(), a, b)
+	end
+
+	assert(mat4.is_mat4(b) or #b == 4, "__mul: Wrong argument type for right hand operand. (<cpml.mat4> or table #4 expected)")
+
+	if mat4.is_mat4(b) then
+		return new():mul(a, b)
+	end
+
+	return mat4.mul_vec4({}, a, b)
+end
+
+if status then
+	xpcall(function() -- Allow this to silently fail; assume failure means someone messed with package.loaded
+		ffi.metatype(new, mat4_mt)
+	end, function() end)
+end
+
+return setmetatable({}, mat4_mt)

libs/cpml/mesh.lua

@@ -0,0 +1,51 @@
+--- Mesh utilities
+-- @module mesh
+
+local modules = (...):gsub('%.[^%.]+$', '') .. "."
+local vec3    = require(modules .. "vec3")
+local mesh    = {}
+
+-- vertices is an arbitrary list of vec3s
+function mesh.average(vertices)
+	local out = vec3()
+	for _, v in ipairs(vertices) do
+		out = out + v
+	end
+	return out / #vertices
+end
+
+-- triangle[1] is a vec3
+-- triangle[2] is a vec3
+-- triangle[3] is a vec3
+function mesh.normal(triangle)
+	local ba = triangle[2] - triangle[1]
+	local ca = triangle[3] - triangle[1]
+	return ba:cross(ca):normalize()
+end
+
+-- triangle[1] is a vec3
+-- triangle[2] is a vec3
+-- triangle[3] is a vec3
+function mesh.plane_from_triangle(triangle)
+	return {
+		origin = triangle[1],
+		normal = mesh.normal(triangle)
+	}
+end
+
+-- plane.origin is a vec3
+-- plane.normal is a vec3
+-- direction    is a vec3
+function mesh.is_front_facing(plane, direction)
+	return plane.normal:dot(direction) >= 0
+end
+
+-- point        is a vec3
+-- plane.origin is a vec3
+-- plane.normal is a vec3
+-- plane.dot    is a number
+function mesh.signed_distance(point, plane)
+	return point:dot(plane.normal) - plane.normal:dot(plane.origin)
+end
+
+return mesh

libs/cpml/octree.lua

@@ -0,0 +1,634 @@
+-- https://github.com/Nition/UnityOctree
+-- https://github.com/Nition/UnityOctree/blob/master/LICENCE
+-- https://github.com/Nition/UnityOctree/blob/master/Scripts/BoundsOctree.cs
+-- https://github.com/Nition/UnityOctree/blob/master/Scripts/BoundsOctreeNode.cs
+
+--- Octree
+-- @module octree
+
+local modules = (...):gsub('%.[^%.]+$', '') .. "."
+local intersect  = require(modules .. "intersect")
+local mat4       = require(modules .. "mat4")
+local utils      = require(modules .. "utils")
+local vec3       = require(modules .. "vec3")
+local Octree     = {}
+local OctreeNode = {}
+local Node
+
+Octree.__index     = Octree
+OctreeNode.__index = OctreeNode
+
+--== Octree ==--
+
+--- Constructor for the bounds octree.
+-- @param initialWorldSize Size of the sides of the initial node, in metres. The octree will never shrink smaller than this
+-- @param initialWorldPos Position of the centre of the initial node
+-- @param minNodeSize Nodes will stop splitting if the new nodes would be smaller than this (metres)
+-- @param looseness Clamped between 1 and 2. Values > 1 let nodes overlap
+local function new(initialWorldSize, initialWorldPos, minNodeSize, looseness)
+	local tree = setmetatable({}, Octree)
+
+	if minNodeSize > initialWorldSize then
+		print("Minimum node size must be at least as big as the initial world size. Was: " .. minNodeSize .. " Adjusted to: " .. initialWorldSize)
+		minNodeSize = initialWorldSize
+	end
+
+	-- The total amount of objects currently in the tree
+	tree.count = 0
+
+	-- Size that the octree was on creation
+	tree.initialSize = initialWorldSize
+
+	-- Minimum side length that a node can be - essentially an alternative to having a max depth
+	tree.minSize = minNodeSize
+
+	-- Should be a value between 1 and 2. A multiplier for the base size of a node.
+	-- 1.0 is a "normal" octree, while values > 1 have overlap
+	tree.looseness = utils.clamp(looseness, 1, 2)
+
+	-- Root node of the octree
+	tree.rootNode = Node(tree.initialSize, tree.minSize, tree.looseness, initialWorldPos)
+
+	return tree
+end
+
+--- Used when growing the octree. Works out where the old root node would fit inside a new, larger root node.
+-- @param xDir X direction of growth. 1 or -1
+-- @param yDir Y direction of growth. 1 or -1
+-- @param zDir Z direction of growth. 1 or -1
+-- @return Octant where the root node should be
+local function get_root_pos_index(xDir, yDir, zDir)
+	local result = xDir > 0 and 1 or 0
+	if yDir < 0 then return result + 4 end
+	if zDir > 0 then return result + 2 end
+end
+
+--- Add an object.
+-- @param obj Object to add
+-- @param objBounds 3D bounding box around the object
+function Octree:add(obj, objBounds)
+	-- Add object or expand the octree until it can be added
+	local count = 0 -- Safety check against infinite/excessive growth
+
+	while not self.rootNode:add(obj, objBounds) do
+		count = count + 1
+		self:grow(objBounds.center - self.rootNode.center)
+
+		if count > 20 then
+			print("Aborted Add operation as it seemed to be going on forever (" .. count - 1 .. ") attempts at growing the octree.")
+			return
+		end
+
+		self.count = self.count + 1
+	end
+end
+
+--- Remove an object. Makes the assumption that the object only exists once in the tree.
+-- @param obj Object to remove
+-- @return bool True if the object was removed successfully
+function Octree:remove(obj)
+	local removed = self.rootNode:remove(obj)
+
+	-- See if we can shrink the octree down now that we've removed the item
+	if removed then
+		self.count = self.count - 1
+		self:shrink()
+	end
+
+	return removed
+end
+
+--- Check if the specified bounds intersect with anything in the tree. See also: get_colliding.
+-- @param checkBounds bounds to check
+-- @return bool True if there was a collision
+function Octree:is_colliding(checkBounds)
+	return self.rootNode:is_colliding(checkBounds)
+end
+
+--- Returns an array of objects that intersect with the specified bounds, if any. Otherwise returns an empty array. See also: is_colliding.
+-- @param checkBounds bounds to check
+-- @return table Objects that intersect with the specified bounds
+function Octree:get_colliding(checkBounds)
+	return self.rootNode:get_colliding(checkBounds)
+end
+
+--- Cast a ray through the node and its children
+-- @param ray Ray with a position and a direction
+-- @param func Function to execute on any objects within child nodes
+-- @param out Table to store results of func in
+-- @return boolean True if an intersect detected
+function Octree:cast_ray(ray, func, out)
+	assert(func)
+	return self.rootNode:cast_ray(ray, func, out)
+end
+
+--- Draws node boundaries visually for debugging.
+function Octree:draw_bounds(cube)
+	self.rootNode:draw_bounds(cube)
+end
+
+--- Draws the bounds of all objects in the tree visually for debugging.
+function Octree:draw_objects(cube, filter)
+	self.rootNode:draw_objects(cube, filter)
+end
+
+--- Grow the octree to fit in all objects.
+-- @param direction Direction to grow
+function Octree:grow(direction)
+	local xDirection = direction.x >= 0 and 1 or -1
+	local yDirection = direction.y >= 0 and 1 or -1
+	local zDirection = direction.z >= 0 and 1 or -1
+
+	local oldRoot   = self.rootNode
+	local half      = self.rootNode.baseLength / 2
+	local newLength = self.rootNode.baseLength * 2
+	local newCenter = self.rootNode.center + vec3(xDirection * half, yDirection * half, zDirection * half)
+
+	-- Create a new, bigger octree root node
+	self.rootNode = Node(newLength, self.minSize, self.looseness, newCenter)
+
+	-- Create 7 new octree children to go with the old root as children of the new root
+	local rootPos  = get_root_pos_index(xDirection, yDirection, zDirection)
+	local children = {}
+
+	for i = 0, 7 do
+		if i == rootPos then
+			children[i+1] = oldRoot
+		else
+			xDirection  = i % 2 == 0 and -1 or 1
+			yDirection  = i > 3 and -1 or 1
+			zDirection  = (i < 2 or (i > 3 and i < 6)) and -1 or 1
+			children[i+1] = Node(self.rootNode.baseLength, self.minSize, self.looseness, newCenter + vec3(xDirection * half, yDirection * half, zDirection * half))
+		end
+	end
+
+	-- Attach the new children to the new root node
+	self.rootNode:set_children(children)
+end
+
+--- Shrink the octree if possible, else leave it the same.
+function Octree:shrink()
+	self.rootNode = self.rootNode:shrink_if_possible(self.initialSize)
+end
+
+--== Octree Node ==--
+
+--- Constructor.
+-- @param baseLength Length of this node, not taking looseness into account
+-- @param minSize Minimum size of nodes in this octree
+-- @param looseness Multiplier for baseLengthVal to get the actual size
+-- @param center Centre position of this node
+local function new_node(baseLength, minSize, looseness, center)
+	local node = setmetatable({}, OctreeNode)
+
+	-- Objects in this node
+	node.objects = {}
+
+	-- Child nodes
+	node.children = {}
+
+	-- If there are already numObjectsAllowed in a node, we split it into children
+	-- A generally good number seems to be something around 8-15
+	node.numObjectsAllowed = 8
+
+	node:set_values(baseLength, minSize, looseness, center)
+
+	return node
+end
+
+local function new_bound(center, size)
+	return {
+		center = center,
+		size   = size,
+		min    = center - (size / 2),
+		max    = center + (size / 2)
+	}
+end
+
+--- Add an object.
+-- @param obj Object to add
+-- @param objBounds 3D bounding box around the object
+-- @return boolean True if the object fits entirely within this node
+function OctreeNode:add(obj, objBounds)
+	if not intersect.encapsulate_aabb(self.bounds, objBounds) then
+		return false
+	end
+
+	-- We know it fits at this level if we've got this far
+	-- Just add if few objects are here, or children would be below min size
+	if #self.objects < self.numObjectsAllowed
+	or self.baseLength / 2 < self.minSize then
+		table.insert(self.objects, {
+			data   = obj,
+			bounds = objBounds
+		})
+	else
+		-- Fits at this level, but we can go deeper. Would it fit there?
+
+		local best_fit_child
+
+		-- Create the 8 children
+		if #self.children == 0 then
+			self:split()
+
+			if #self.children == 0 then
+				print("Child creation failed for an unknown reason. Early exit.")
+				return false
+			end
+
+			-- Now that we have the new children, see if this node's existing objects would fit there
+			for i = #self.objects, 1, -1 do
+				local object = self.objects[i]
+				-- Find which child the object is closest to based on where the
+				-- object's center is located in relation to the octree's center.
+				best_fit_child = self:best_fit_child(object.bounds)
+
+				-- Does it fit?
+				if intersect.encapsulate_aabb(self.children[best_fit_child].bounds, object.bounds) then
+					self.children[best_fit_child]:add(object.data, object.bounds) -- Go a level deeper
+					table.remove(self.objects, i) -- Remove from here
+				end
+			end
+		end
+
+		-- Now handle the new object we're adding now
+		best_fit_child = self:best_fit_child(objBounds)
+
+		if intersect.encapsulate_aabb(self.children[best_fit_child].bounds, objBounds) then
+			self.children[best_fit_child]:add(obj, objBounds)
+		else
+			table.insert(self.objects, {
+				data   = obj,
+				bounds = objBounds
+			})
+		end
+	end
+
+	return true
+end
+
+--- Remove an object. Makes the assumption that the object only exists once in the tree.
+-- @param obj Object to remove
+-- @return boolean True if the object was removed successfully
+function OctreeNode:remove(obj)
+	local removed = false
+
+	for i, object in ipairs(self.objects) do
+		if object == obj then
+			removed = table.remove(self.objects, i) and true or false
+			break
+		end
+	end
+
+	if not removed then
+		for _, child in ipairs(self.children) do
+			removed = child:remove(obj)
+			if removed then break end
+		end
+	end
+
+	if removed then
+		-- Check if we should merge nodes now that we've removed an item
+		if self:should_merge() then
+			self:merge()
+		end
+	end
+
+	return removed
+end
+
+--- Check if the specified bounds intersect with anything in the tree. See also: get_colliding.
+-- @param checkBounds Bounds to check
+-- @return boolean True if there was a collision
+function OctreeNode:is_colliding(checkBounds)
+	-- Are the input bounds at least partially in this node?
+	if not intersect.aabb_aabb(self.bounds, checkBounds) then
+		return false
+	end
+
+	-- Check against any objects in this node
+	for _, object in ipairs(self.objects) do
+		if intersect.aabb_aabb(object.bounds, checkBounds) then
+			return true
+		end
+	end
+
+	-- Check children
+	for _, child in ipairs(self.children) do
+		if child:is_colliding(checkBounds) then
+			return true
+		end
+	end
+
+	return false
+end
+
+--- Returns an array of objects that intersect with the specified bounds, if any. Otherwise returns an empty array. See also: is_colliding.
+-- @param checkBounds Bounds to check. Passing by ref as it improve performance with structs
+-- @param results List results
+-- @return table Objects that intersect with the specified bounds
+function OctreeNode:get_colliding(checkBounds, results)
+	results = results or {}
+
+	-- Are the input bounds at least partially in this node?
+	if not intersect.aabb_aabb(self.bounds, checkBounds) then
+		return results
+	end
+
+	-- Check against any objects in this node
+	for _, object in ipairs(self.objects) do
+		if intersect.aabb_aabb(object.bounds, checkBounds) then
+			table.insert(results, object.data)
+		end
+	end
+
+	-- Check children
+	for _, child in ipairs(self.children) do
+		results = child:get_colliding(checkBounds, results)
+	end
+
+	return results
+end
+
+--- Cast a ray through the node and its children
+-- @param ray Ray with a position and a direction
+-- @param func Function to execute on any objects within child nodes
+-- @param out Table to store results of func in
+-- @param depth (used internally)
+-- @return boolean True if an intersect is detected
+function OctreeNode:cast_ray(ray, func, out, depth)
+	depth = depth or 1
+
+	if intersect.ray_aabb(ray, self.bounds) then
+		if #self.objects > 0 then
+			local hit = func(ray, self.objects, out)
+
+			if hit then
+				return hit
+			end
+		end
+
+		for _, child in ipairs(self.children) do
+			local hit = child:cast_ray(ray, func, out, depth + 1)
+
+			if hit then
+				return hit
+			end
+		end
+	end
+
+	return false
+end
+
+--- Set the 8 children of this octree.
+-- @param childOctrees The 8 new child nodes
+function OctreeNode:set_children(childOctrees)
+	if #childOctrees ~= 8 then
+		print("Child octree array must be length 8. Was length: " .. #childOctrees)
+		return
+	end
+
+	self.children = childOctrees
+end
+
+--- We can shrink the octree if:
+--- - This node is >= double minLength in length
+--- - All objects in the root node are within one octant
+--- - This node doesn't have children, or does but 7/8 children are empty
+--- We can also shrink it if there are no objects left at all!
+-- @param minLength Minimum dimensions of a node in this octree
+-- @return table The new root, or the existing one if we didn't shrink
+function OctreeNode:shrink_if_possible(minLength)
+	if self.baseLength < 2 * minLength then
+		return self
+	end
+
+	if #self.objects == 0 and #self.children == 0 then
+		return self
+	end
+
+	-- Check objects in root
+	local bestFit = 0
+
+	for i, object in ipairs(self.objects) do
+		local newBestFit = self:best_fit_child(object.bounds)
+
+		if i == 1 or newBestFit == bestFit then
+			-- In same octant as the other(s). Does it fit completely inside that octant?
+			if intersect.encapsulate_aabb(self.childBounds[newBestFit], object.bounds) then
+				if bestFit < 1 then
+					bestFit = newBestFit
+				end
+			else
+				-- Nope, so we can't reduce. Otherwise we continue
+				return self
+			end
+		else
+			return self -- Can't reduce - objects fit in different octants
+		end
+	end
+
+	-- Check objects in children if there are any
+	if #self.children > 0 then
+		local childHadContent = false
+
+		for i, child in ipairs(self.children) do
+			if child:has_any_objects() then
+				if childHadContent then
+					return self -- Can't shrink - another child had content already
+				end
+
+				if bestFit > 0 and bestFit ~= i then
+					return self -- Can't reduce - objects in root are in a different octant to objects in child
+				end
+
+				childHadContent = true
+				bestFit = i
+			end
+		end
+	end
+
+	-- Can reduce
+	if #self.children == 0 then
+		-- We don't have any children, so just shrink this node to the new size
+		-- We already know that everything will still fit in it
+		self:set_values(self.baseLength / 2, self.minSize, self.looseness, self.childBounds[bestFit].center)
+		return self
+	end
+
+	-- We have children. Use the appropriate child as the new root node
+	return self.children[bestFit]
+end
+
+--- Set values for this node.
+-- @param baseLength Length of this node, not taking looseness into account
+-- @param minSize Minimum size of nodes in this octree
+-- @param looseness Multiplier for baseLengthVal to get the actual size
+-- @param center Centre position of this node
+function OctreeNode:set_values(baseLength, minSize, looseness, center)
+	-- Length of this node if it has a looseness of 1.0
+	self.baseLength = baseLength
+
+	-- Minimum size for a node in this octree
+	self.minSize = minSize
+
+	-- Looseness value for this node
+	self.looseness = looseness
+
+	-- Centre of this node
+	self.center = center
+
+	-- Actual length of sides, taking the looseness value into account
+	self.adjLength = self.looseness * self.baseLength
+
+	-- Create the bounding box.
+	self.size = vec3(self.adjLength, self.adjLength, self.adjLength)
+
+	-- Bounding box that represents this node
+	self.bounds = new_bound(self.center, self.size)
+
+	self.quarter           = self.baseLength / 4
+	self.childActualLength = (self.baseLength / 2) * self.looseness
+	self.childActualSize   = vec3(self.childActualLength, self.childActualLength, self.childActualLength)
+
+	-- Bounds of potential children to this node. These are actual size (with looseness taken into account), not base size
+	self.childBounds =  {
+		new_bound(self.center + vec3(-self.quarter,  self.quarter, -self.quarter), self.childActualSize),
+		new_bound(self.center + vec3( self.quarter,  self.quarter, -self.quarter), self.childActualSize),
+		new_bound(self.center + vec3(-self.quarter,  self.quarter,  self.quarter), self.childActualSize),
+		new_bound(self.center + vec3( self.quarter,  self.quarter,  self.quarter), self.childActualSize),
+		new_bound(self.center + vec3(-self.quarter, -self.quarter, -self.quarter), self.childActualSize),
+		new_bound(self.center + vec3( self.quarter, -self.quarter, -self.quarter), self.childActualSize),
+		new_bound(self.center + vec3(-self.quarter, -self.quarter,  self.quarter), self.childActualSize),
+		new_bound(self.center + vec3( self.quarter, -self.quarter,  self.quarter), self.childActualSize)
+	}
+end
+
+--- Splits the octree into eight children.
+function OctreeNode:split()
+	if #self.children > 0 then return end
+
+	local quarter   = self.baseLength / 4
+	local newLength = self.baseLength / 2
+
+	table.insert(self.children, Node(newLength, self.minSize, self.looseness, self.center + vec3(-quarter,  quarter, -quarter)))
+	table.insert(self.children, Node(newLength, self.minSize, self.looseness, self.center + vec3( quarter,  quarter, -quarter)))
+	table.insert(self.children, Node(newLength, self.minSize, self.looseness, self.center + vec3(-quarter,  quarter,  quarter)))
+	table.insert(self.children, Node(newLength, self.minSize, self.looseness, self.center + vec3( quarter,  quarter,  quarter)))
+	table.insert(self.children, Node(newLength, self.minSize, self.looseness, self.center + vec3(-quarter, -quarter, -quarter)))
+	table.insert(self.children, Node(newLength, self.minSize, self.looseness, self.center + vec3( quarter, -quarter, -quarter)))
+	table.insert(self.children, Node(newLength, self.minSize, self.looseness, self.center + vec3(-quarter, -quarter,  quarter)))
+	table.insert(self.children, Node(newLength, self.minSize, self.looseness, self.center + vec3( quarter, -quarter,  quarter)))
+end
+
+--- Merge all children into this node - the opposite of Split.
+--- Note: We only have to check one level down since a merge will never happen if the children already have children,
+--- since THAT won't happen unless there are already too many objects to merge.
+function OctreeNode:merge()
+	for _, child in ipairs(self.children) do
+		for _, object in ipairs(child.objects) do
+			table.insert(self.objects, object)
+		end
+	end
+
+	-- Remove the child nodes (and the objects in them - they've been added elsewhere now)
+	self.children = {}
+end
+
+--- Find which child node this object would be most likely to fit in.
+-- @param objBounds The object's bounds
+-- @return number One of the eight child octants
+function OctreeNode:best_fit_child(objBounds)
+	return (objBounds.center.x <= self.center.x and 0 or 1) + (objBounds.center.y >= self.center.y and 0 or 4) + (objBounds.center.z <= self.center.z and 0 or 2) + 1
+end
+
+--- Checks if there are few enough objects in this node and its children that the children should all be merged into this.
+-- @return boolean True there are less or the same abount of objects in this and its children than numObjectsAllowed
+function OctreeNode:should_merge()
+	local totalObjects = #self.objects
+
+	for _, child in ipairs(self.children) do
+		if #child.children > 0 then
+			-- If any of the *children* have children, there are definitely too many to merge,
+			-- or the child would have been merged already
+			return false
+		end
+
+		totalObjects = totalObjects + #child.objects
+	end
+
+	return totalObjects <= self.numObjectsAllowed
+end
+
+--- Checks if this node or anything below it has something in it.
+-- @return boolean True if this node or any of its children, grandchildren etc have something in the
+function OctreeNode:has_any_objects()
+	if #self.objects > 0 then return true end
+
+	for _, child in ipairs(self.children) do
+		if child:has_any_objects() then return true end
+	end
+
+	return false
+end
+
+--- Draws node boundaries visually for debugging.
+-- @param cube Cube model to draw
+-- @param depth Used for recurcive calls to this method
+function OctreeNode:draw_bounds(cube, depth)
+	depth = depth or 0
+	local tint = depth / 7 -- Will eventually get values > 1. Color rounds to 1 automatically
+
+	love.graphics.setColor(tint * 255, 0, (1 - tint) * 255)
+	local m = mat4()
+		:translate(self.center)
+		:scale(vec3(self.adjLength, self.adjLength, self.adjLength))
+
+	love.graphics.updateMatrix("transform", m)
+	love.graphics.setWireframe(true)
+	love.graphics.draw(cube)
+	love.graphics.setWireframe(false)
+
+	for _, child in ipairs(self.children) do
+		child:draw_bounds(cube, depth + 1)
+	end
+
+	love.graphics.setColor(255, 255, 255)
+end
+
+--- Draws the bounds of all objects in the tree visually for debugging.
+-- @param cube Cube model to draw
+-- @param filter a function returning true or false to determine visibility.
+function OctreeNode:draw_objects(cube, filter)
+	local tint = self.baseLength / 20
+	love.graphics.setColor(0, (1 - tint) * 255, tint * 255, 63)
+
+	for _, object in ipairs(self.objects) do
+		if filter and filter(object.data) or not filter then
+			local m = mat4()
+				:translate(object.bounds.center)
+				:scale(object.bounds.size)
+
+			love.graphics.updateMatrix("transform", m)
+			love.graphics.draw(cube)
+		end
+	end
+
+	for _, child in ipairs(self.children) do
+		child:draw_objects(cube, filter)
+	end
+
+	love.graphics.setColor(255, 255, 255)
+end
+
+Node = setmetatable({
+	new = new_node
+}, {
+	__call = function(_, ...) return new_node(...) end
+})
+
+return setmetatable({
+	new = new
+}, {
+	__call = function(_, ...) return new(...) end
+})

libs/cpml/quat.lua

@@ -0,0 +1,498 @@
+--- A quaternion and associated utilities.
+-- @module quat
+
+local modules       = (...):gsub('%.[^%.]+$', '') .. "."
+local constants     = require(modules .. "constants")
+local vec3          = require(modules .. "vec3")
+local precond       = require(modules .. "_private_precond")
+local private       = require(modules .. "_private_utils")
+local DOT_THRESHOLD = constants.DOT_THRESHOLD
+local DBL_EPSILON   = constants.DBL_EPSILON
+local acos          = math.acos
+local cos           = math.cos
+local sin           = math.sin
+local min           = math.min
+local max           = math.max
+local sqrt          = math.sqrt
+local quat          = {}
+local quat_mt       = {}
+
+-- Private constructor.
+local function new(x, y, z, w)
+	return setmetatable({
+		x = x or 0,
+		y = y or 0,
+		z = z or 0,
+		w = w or 1
+	}, quat_mt)
+end
+
+-- Do the check to see if JIT is enabled. If so use the optimized FFI structs.
+local status, ffi
+if type(jit) == "table" and jit.status() then
+	status, ffi = pcall(require, "ffi")
+	if status then
+		ffi.cdef "typedef struct { double x, y, z, w;} cpml_quat;"
+		new = ffi.typeof("cpml_quat")
+	end
+end
+
+-- Statically allocate a temporary variable used in some of our functions.
+local tmp = new()
+local qv, uv, uuv = vec3(), vec3(), vec3()
+
+--- Constants
+-- @table quat
+-- @field unit Unit quaternion
+-- @field zero Empty quaternion
+quat.unit = new(0, 0, 0, 1)
+quat.zero = new(0, 0, 0, 0)
+
+--- The public constructor.
+-- @param x Can be of two types: </br>
+-- number x X component
+-- table {x, y, z, w} or {x=x, y=y, z=z, w=w}
+-- @tparam number y Y component
+-- @tparam number z Z component
+-- @tparam number w W component
+-- @treturn quat out
+function quat.new(x, y, z, w)
+	-- number, number, number, number
+	if x and y and z and w then
+		precond.typeof(x, "number", "new: Wrong argument type for x")
+		precond.typeof(y, "number", "new: Wrong argument type for y")
+		precond.typeof(z, "number", "new: Wrong argument type for z")
+		precond.typeof(w, "number", "new: Wrong argument type for w")
+
+		return new(x, y, z, w)
+
+	-- {x, y, z, w} or {x=x, y=y, z=z, w=w}
+	elseif type(x) == "table" then
+		local xx, yy, zz, ww = x.x or x[1], x.y or x[2], x.z or x[3], x.w or x[4]
+		precond.typeof(xx, "number", "new: Wrong argument type for x")
+		precond.typeof(yy, "number", "new: Wrong argument type for y")
+		precond.typeof(zz, "number", "new: Wrong argument type for z")
+		precond.typeof(ww, "number", "new: Wrong argument type for w")
+
+		return new(xx, yy, zz, ww)
+	end
+
+	return new(0, 0, 0, 1)
+end
+
+--- Create a quaternion from an angle/axis pair.
+-- @tparam number angle Angle (in radians)
+-- @param axis/x -- Can be of two types, a vec3 axis, or the x component of that axis
+-- @param y axis -- y component of axis (optional, only if x component param used)
+-- @param z axis -- z component of axis (optional, only if x component param used)
+-- @treturn quat out
+function quat.from_angle_axis(angle, axis, a3, a4)
+	if axis and a3 and a4 then
+		local x, y, z = axis, a3, a4
+		local s = sin(angle * 0.5)
+		local c = cos(angle * 0.5)
+		return new(x * s, y * s, z * s, c)
+	else
+		return quat.from_angle_axis(angle, axis.x, axis.y, axis.z)
+	end
+end
+
+--- Create a quaternion from a normal/up vector pair.
+-- @tparam vec3 normal
+-- @tparam vec3 up (optional)
+-- @treturn quat out
+function quat.from_direction(normal, up)
+	local u = up or vec3.unit_z
+	local n = normal:normalize()
+	local a = u:cross(n)
+	local d = u:dot(n)
+	return new(a.x, a.y, a.z, d + 1)
+end
+
+--- Clone a quaternion.
+-- @tparam quat a Quaternion to clone
+-- @treturn quat out
+function quat.clone(a)
+	return new(a.x, a.y, a.z, a.w)
+end
+
+--- Add two quaternions.
+-- @tparam quat a Left hand operand
+-- @tparam quat b Right hand operand
+-- @treturn quat out
+function quat.add(a, b)
+	return new(
+		a.x + b.x,
+		a.y + b.y,
+		a.z + b.z,
+		a.w + b.w
+	)
+end
+
+--- Subtract a quaternion from another.
+-- @tparam quat a Left hand operand
+-- @tparam quat b Right hand operand
+-- @treturn quat out
+function quat.sub(a, b)
+	return new(
+		a.x - b.x,
+		a.y - b.y,
+		a.z - b.z,
+		a.w - b.w
+	)
+end
+
+--- Multiply two quaternions.
+-- @tparam quat a Left hand operand
+-- @tparam quat b Right hand operand
+-- @treturn quat quaternion equivalent to "apply b, then a"
+function quat.mul(a, b)
+	return new(
+		a.x * b.w + a.w * b.x + a.y * b.z - a.z * b.y,
+		a.y * b.w + a.w * b.y + a.z * b.x - a.x * b.z,
+		a.z * b.w + a.w * b.z + a.x * b.y - a.y * b.x,
+		a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z
+	)
+end
+
+--- Multiply a quaternion and a vec3.
+-- @tparam quat a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn vec3 out
+function quat.mul_vec3(a, b)
+	qv.x = a.x
+	qv.y = a.y
+	qv.z = a.z
+	uv   = qv:cross(b)
+	uuv  = qv:cross(uv)
+	return b + ((uv * a.w) + uuv) * 2
+end
+
+--- Raise a normalized quaternion to a scalar power.
+-- @tparam quat a Left hand operand (should be a unit quaternion)
+-- @tparam number s Right hand operand
+-- @treturn quat out
+function quat.pow(a, s)
+	-- Do it as a slerp between identity and a (code borrowed from slerp)
+	if a.w < 0 then
+		a   = -a
+	end
+	local dot = a.w
+
+	dot = min(max(dot, -1), 1)
+
+	local theta = acos(dot) * s
+	local c = new(a.x, a.y, a.z, 0):normalize() * sin(theta)
+	c.w = cos(theta)
+	return c
+end
+
+--- Normalize a quaternion.
+-- @tparam quat a Quaternion to normalize
+-- @treturn quat out
+function quat.normalize(a)
+	if a:is_zero() then
+		return new(0, 0, 0, 0)
+	end
+	return a:scale(1 / a:len())
+end
+
+--- Get the dot product of two quaternions.
+-- @tparam quat a Left hand operand
+-- @tparam quat b Right hand operand
+-- @treturn number dot
+function quat.dot(a, b)
+	return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w
+end
+
+--- Return the length of a quaternion.
+-- @tparam quat a Quaternion to get length of
+-- @treturn number len
+function quat.len(a)
+	return sqrt(a.x * a.x + a.y * a.y + a.z * a.z + a.w * a.w)
+end
+
+--- Return the squared length of a quaternion.
+-- @tparam quat a Quaternion to get length of
+-- @treturn number len
+function quat.len2(a)
+	return a.x * a.x + a.y * a.y + a.z * a.z + a.w * a.w
+end
+
+--- Multiply a quaternion by a scalar.
+-- @tparam quat a Left hand operand
+-- @tparam number s Right hand operand
+-- @treturn quat out
+function quat.scale(a, s)
+	return new(
+		a.x * s,
+		a.y * s,
+		a.z * s,
+		a.w * s
+	)
+end
+
+--- Alias of from_angle_axis.
+-- @tparam number angle Angle (in radians)
+-- @param axis/x -- Can be of two types, a vec3 axis, or the x component of that axis
+-- @param y axis -- y component of axis (optional, only if x component param used)
+-- @param z axis -- z component of axis (optional, only if x component param used)
+-- @treturn quat out
+function quat.rotate(angle, axis, a3, a4)
+	return quat.from_angle_axis(angle, axis, a3, a4)
+end
+
+--- Return the conjugate of a quaternion.
+-- @tparam quat a Quaternion to conjugate
+-- @treturn quat out
+function quat.conjugate(a)
+	return new(-a.x, -a.y, -a.z, a.w)
+end
+
+--- Return the inverse of a quaternion.
+-- @tparam quat a Quaternion to invert
+-- @treturn quat out
+function quat.inverse(a)
+	tmp.x = -a.x
+	tmp.y = -a.y
+	tmp.z = -a.z
+	tmp.w =  a.w
+	return tmp:normalize()
+end
+
+--- Return the reciprocal of a quaternion.
+-- @tparam quat a Quaternion to reciprocate
+-- @treturn quat out
+function quat.reciprocal(a)
+	if a:is_zero() then
+		error("Cannot reciprocate a zero quaternion")
+		return false
+	end
+
+	tmp.x = -a.x
+	tmp.y = -a.y
+	tmp.z = -a.z
+	tmp.w =  a.w
+
+	return tmp:scale(1 / a:len2())
+end
+
+--- Lerp between two quaternions.
+-- @tparam quat a Left hand operand
+-- @tparam quat b Right hand operand
+-- @tparam number s Step value
+-- @treturn quat out
+function quat.lerp(a, b, s)
+	return (a + (b - a) * s):normalize()
+end
+
+--- Slerp between two quaternions.
+-- @tparam quat a Left hand operand
+-- @tparam quat b Right hand operand
+-- @tparam number s Step value
+-- @treturn quat out
+function quat.slerp(a, b, s)
+	local dot = a:dot(b)
+
+	if dot < 0 then
+		a   = -a
+		dot = -dot
+	end
+
+	if dot > DOT_THRESHOLD then
+		return a:lerp(b, s)
+	end
+
+	dot = min(max(dot, -1), 1)
+
+	local theta = acos(dot) * s
+	local c = (b - a * dot):normalize()
+	return a * cos(theta) + c * sin(theta)
+end
+
+--- Unpack a quaternion into individual components.
+-- @tparam quat a Quaternion to unpack
+-- @treturn number x
+-- @treturn number y
+-- @treturn number z
+-- @treturn number w
+function quat.unpack(a)
+	return a.x, a.y, a.z, a.w
+end
+
+--- Return a boolean showing if a table is or is not a quat.
+-- @tparam quat a Quaternion to be tested
+-- @treturn boolean is_quat
+function quat.is_quat(a)
+	if type(a) == "cdata" then
+		return ffi.istype("cpml_quat", a)
+	end
+
+	return
+		type(a)   == "table"  and
+		type(a.x) == "number" and
+		type(a.y) == "number" and
+		type(a.z) == "number" and
+		type(a.w) == "number"
+end
+
+--- Return a boolean showing if a table is or is not a zero quat.
+-- @tparam quat a Quaternion to be tested
+-- @treturn boolean is_zero
+function quat.is_zero(a)
+	return
+		a.x == 0 and
+		a.y == 0 and
+		a.z == 0 and
+		a.w == 0
+end
+
+--- Return a boolean showing if a table is or is not a real quat.
+-- @tparam quat a Quaternion to be tested
+-- @treturn boolean is_real
+function quat.is_real(a)
+	return
+		a.x == 0 and
+		a.y == 0 and
+		a.z == 0
+end
+
+--- Return a boolean showing if a table is or is not an imaginary quat.
+-- @tparam quat a Quaternion to be tested
+-- @treturn boolean is_imaginary
+function quat.is_imaginary(a)
+	return a.w == 0
+end
+
+--- Return whether any component is NaN
+-- @tparam quat a Quaternion to be tested
+-- @treturn boolean if x,y,z, or w is NaN
+function quat.has_nan(a)
+	return private.is_nan(a.x) or
+		private.is_nan(a.y) or
+		private.is_nan(a.z) or
+		private.is_nan(a.w)
+end
+
+--- Convert a quaternion into an angle plus axis components.
+-- @tparam quat a Quaternion to convert
+-- @tparam identityAxis vec3 of axis to use on identity/degenerate quaternions (optional, default returns 0,0,0,1)
+-- @treturn number angle
+-- @treturn x axis-x
+-- @treturn y axis-y
+-- @treturn z axis-z
+function quat.to_angle_axis_unpack(a, identityAxis)
+	if a.w > 1 or a.w < -1 then
+		a = a:normalize()
+	end
+
+	-- If length of xyz components is less than DBL_EPSILON, this is zero or close enough (an identity quaternion)
+	-- Normally an identity quat would return a nonsense answer, so we return an arbitrary zero rotation early.
+	-- FIXME: Is it safe to assume there are *no* valid quaternions with nonzero degenerate lengths?
+	if a.x*a.x + a.y*a.y + a.z*a.z < constants.DBL_EPSILON*constants.DBL_EPSILON then
+		if identityAxis then
+			return 0,identityAxis:unpack()
+		else
+			return 0,0,0,1
+		end
+	end
+
+	local x, y, z
+	local angle = 2 * acos(a.w)
+	local s     = sqrt(1 - a.w * a.w)
+
+	if s < DBL_EPSILON then
+		x = a.x
+		y = a.y
+		z = a.z
+	else
+		x = a.x / s
+		y = a.y / s
+		z = a.z / s
+	end
+
+	return angle, x, y, z
+end
+
+--- Convert a quaternion into an angle/axis pair.
+-- @tparam quat a Quaternion to convert
+-- @tparam identityAxis vec3 of axis to use on identity/degenerate quaternions (optional, default returns 0,vec3(0,0,1))
+-- @treturn number angle
+-- @treturn vec3 axis
+function quat.to_angle_axis(a, identityAxis)
+	local angle, x, y, z = a:to_angle_axis_unpack(identityAxis)
+	return angle, vec3(x, y, z)
+end
+
+--- Convert a quaternion into a vec3.
+-- @tparam quat a Quaternion to convert
+-- @treturn vec3 out
+function quat.to_vec3(a)
+	return vec3(a.x, a.y, a.z)
+end
+
+--- Return a formatted string.
+-- @tparam quat a Quaternion to be turned into a string
+-- @treturn string formatted
+function quat.to_string(a)
+	return string.format("(%+0.3f,%+0.3f,%+0.3f,%+0.3f)", a.x, a.y, a.z, a.w)
+end
+
+quat_mt.__index    = quat
+quat_mt.__tostring = quat.to_string
+
+function quat_mt.__call(_, x, y, z, w)
+	return quat.new(x, y, z, w)
+end
+
+function quat_mt.__unm(a)
+	return a:scale(-1)
+end
+
+function quat_mt.__eq(a,b)
+	if not quat.is_quat(a) or not quat.is_quat(b) then
+		return false
+	end
+	return a.x == b.x and a.y == b.y and a.z == b.z and a.w == b.w
+end
+
+function quat_mt.__add(a, b)
+	precond.assert(quat.is_quat(a), "__add: Wrong argument type '%s' for left hand operand. (<cpml.quat> expected)", type(a))
+	precond.assert(quat.is_quat(b), "__add: Wrong argument type '%s' for right hand operand. (<cpml.quat> expected)", type(b))
+	return a:add(b)
+end
+
+function quat_mt.__sub(a, b)
+	precond.assert(quat.is_quat(a), "__sub: Wrong argument type '%s' for left hand operand. (<cpml.quat> expected)", type(a))
+	precond.assert(quat.is_quat(b), "__sub: Wrong argument type '%s' for right hand operand. (<cpml.quat> expected)", type(b))
+	return a:sub(b)
+end
+
+function quat_mt.__mul(a, b)
+	precond.assert(quat.is_quat(a), "__mul: Wrong argument type '%s' for left hand operand. (<cpml.quat> expected)", type(a))
+	assert(quat.is_quat(b) or vec3.is_vec3(b) or type(b) == "number", "__mul: Wrong argument type for right hand operand. (<cpml.quat> or <cpml.vec3> or <number> expected)")
+
+	if quat.is_quat(b) then
+		return a:mul(b)
+	end
+
+	if type(b) == "number" then
+		return a:scale(b)
+	end
+
+	return a:mul_vec3(b)
+end
+
+function quat_mt.__pow(a, n)
+	precond.assert(quat.is_quat(a), "__pow: Wrong argument type '%s' for left hand operand. (<cpml.quat> expected)", type(a))
+	precond.typeof(n, "number", "__pow: Wrong argument type for right hand operand.")
+	return a:pow(n)
+end
+
+if status then
+	xpcall(function() -- Allow this to silently fail; assume failure means someone messed with package.loaded
+		ffi.metatype(new, quat_mt)
+	end, function() end)
+end
+
+return setmetatable({}, quat_mt)

libs/cpml/simplex.lua

@@ -0,0 +1,349 @@
+--- Simplex Noise
+-- @module simplex
+
+--
+-- Based on code in "Simplex noise demystified", by Stefan Gustavson
+-- www.itn.liu.se/~stegu/simplexnoise/simplexnoise.pdf
+--
+-- Thanks to Mike Pall for some cleanup and improvements (and for LuaJIT!)
+--
+-- Permission is hereby granted, free of charge, to any person obtaining
+-- a copy of this software and associated documentation files (the
+-- "Software"), to deal in the Software without restriction, including
+-- without limitation the rights to use, copy, modify, merge, publish,
+-- distribute, sublicense, and/or sell copies of the Software, and to
+-- permit persons to whom the Software is furnished to do so, subject to
+-- the following conditions:
+--
+-- The above copyright notice and this permission notice shall be
+-- included in all copies or substantial portions of the Software.
+--
+-- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+-- EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+-- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+-- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+-- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+-- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+-- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+--
+-- [ MIT license: http://www.opensource.org/licenses/mit-license.php ]
+--
+
+if _G.love and _G.love.math then
+	return love.math.noise
+end
+
+-- Bail out with dummy module if FFI is missing.
+local has_ffi, ffi = pcall(require, "ffi")
+if not has_ffi then
+	return function()
+		return 0
+	end
+end
+
+-- Modules --
+local bit = require("bit")
+
+-- Imports --
+local band   = bit.band
+local bor    = bit.bor
+local floor  = math.floor
+local lshift = bit.lshift
+local max    = math.max
+local rshift = bit.rshift
+
+-- Permutation of 0-255, replicated to allow easy indexing with sums of two bytes --
+local Perms = ffi.new("uint8_t[512]", {
+	151, 160, 137,  91,  90,  15, 131,  13, 201,  95,  96,  53, 194, 233,   7, 225,
+	140,  36, 103,  30,  69, 142,   8,  99,  37, 240,  21,  10,  23, 190,   6, 148,
+	247, 120, 234,  75,   0,  26, 197,  62,  94, 252, 219, 203, 117,  35,  11,  32,
+	 57, 177,  33,  88, 237, 149,  56,  87, 174,  20, 125, 136, 171, 168,  68, 175,
+	 74, 165,  71, 134, 139,  48,  27, 166,  77, 146, 158, 231,  83, 111, 229, 122,
+	 60, 211, 133, 230, 220, 105,  92,  41,  55,  46, 245,  40, 244, 102, 143,  54,
+	 65,  25,  63, 161,   1, 216,  80,  73, 209,  76, 132, 187, 208,  89,  18, 169,
+	200, 196, 135, 130, 116, 188, 159,  86, 164, 100, 109, 198, 173, 186,   3,  64,
+	 52, 217, 226, 250, 124, 123,   5, 202,  38, 147, 118, 126, 255,  82,  85, 212,
+	207, 206,  59, 227,  47,  16,  58,  17, 182, 189,  28,  42, 223, 183, 170, 213,
+	119, 248, 152,   2,  44, 154, 163,  70, 221, 153, 101, 155, 167,  43, 172,   9,
+	129,  22,  39, 253,  19,  98, 108, 110,  79, 113, 224, 232, 178, 185, 112, 104,
+	218, 246,  97, 228, 251,  34, 242, 193, 238, 210, 144,  12, 191, 179, 162, 241,
+	 81,  51, 145, 235, 249,  14, 239, 107,  49, 192, 214,  31, 181, 199, 106, 157,
+	184,  84, 204, 176, 115, 121,  50,  45, 127,   4, 150, 254, 138, 236, 205,  93,
+	222, 114,  67,  29,  24,  72, 243, 141, 128, 195,  78,  66, 215,  61, 156, 180
+})
+
+-- The above, mod 12 for each element --
+local Perms12 = ffi.new("uint8_t[512]")
+
+for i = 0, 255 do
+	local x = Perms[i] % 12
+
+	Perms[i + 256], Perms12[i], Perms12[i + 256] = Perms[i], x, x
+end
+
+-- Gradients for 2D, 3D case --
+local Grads3 = ffi.new("const double[12][3]",
+	{ 1, 1, 0 }, { -1, 1, 0 }, { 1, -1, 0 }, { -1, -1, 0 },
+	{ 1, 0, 1 }, { -1, 0, 1 }, { 1, 0, -1 }, { -1, 0, -1 },
+	{ 0, 1, 1 }, { 0, -1, 1 }, { 0, 1, -1 }, { 0, -1, -1 }
+)
+
+-- 2D weight contribution
+local function GetN2(bx, by, x, y)
+	local t = .5 - x * x - y * y
+	local index = Perms12[bx + Perms[by]]
+
+	return max(0, (t * t) * (t * t)) * (Grads3[index][0] * x + Grads3[index][1] * y)
+end
+
+local function simplex_2d(x, y)
+	--[[
+		2D skew factors:
+		F = (math.sqrt(3) - 1) / 2
+		G = (3 - math.sqrt(3)) / 6
+		G2 = 2 * G - 1
+	]]
+
+	-- Skew the input space to determine which simplex cell we are in.
+	local s = (x + y) * 0.366025403 -- F
+	local ix, iy = floor(x + s), floor(y + s)
+
+	-- Unskew the cell origin back to (x, y) space.
+	local t = (ix + iy) * 0.211324865 -- G
+	local x0 = x + t - ix
+	local y0 = y + t - iy
+
+	-- Calculate the contribution from the two fixed corners.
+	-- A step of (1,0) in (i,j) means a step of (1-G,-G) in (x,y), and
+	-- A step of (0,1) in (i,j) means a step of (-G,1-G) in (x,y).
+	ix, iy = band(ix, 255), band(iy, 255)
+
+	local n0 = GetN2(ix, iy, x0, y0)
+	local n2 = GetN2(ix + 1, iy + 1, x0 - 0.577350270, y0 - 0.577350270) -- G2
+
+	--[[
+		Determine other corner based on simplex (equilateral triangle) we are in:
+		if x0 > y0 then
+			ix, x1 = ix + 1, x1 - 1
+		else
+			iy, y1 = iy + 1, y1 - 1
+		end
+	]]
+	local xi = rshift(floor(y0 - x0), 31) -- y0 < x0
+	local n1 = GetN2(ix + xi, iy + (1 - xi), x0 + 0.211324865 - xi, y0 - 0.788675135 + xi) -- x0 + G - xi, y0 + G - (1 - xi)
+
+	-- Add contributions from each corner to get the final noise value.
+	-- The result is scaled to return values in the interval [-1,1].
+	return 70.1480580019 * (n0 + n1 + n2)
+end
+
+-- 3D weight contribution
+local function GetN3(ix, iy, iz, x, y, z)
+	local t = .6 - x * x - y * y - z * z
+	local index = Perms12[ix + Perms[iy + Perms[iz]]]
+
+	return max(0, (t * t) * (t * t)) * (Grads3[index][0] * x + Grads3[index][1] * y + Grads3[index][2] * z)
+end
+
+local function simplex_3d(x, y, z)
+	--[[
+		3D skew factors:
+		F = 1 / 3
+		G = 1 / 6
+		G2 = 2 * G
+		G3 = 3 * G - 1
+	]]
+
+	-- Skew the input space to determine which simplex cell we are in.
+	local s = (x + y + z) * 0.333333333 -- F
+	local ix, iy, iz = floor(x + s), floor(y + s), floor(z + s)
+
+	-- Unskew the cell origin back to (x, y, z) space.
+	local t = (ix + iy + iz) * 0.166666667 -- G
+	local x0 = x + t - ix
+	local y0 = y + t - iy
+	local z0 = z + t - iz
+
+	-- Calculate the contribution from the two fixed corners.
+	-- A step of (1,0,0) in (i,j,k) means a step of (1-G,-G,-G) in (x,y,z);
+	-- a step of (0,1,0) in (i,j,k) means a step of (-G,1-G,-G) in (x,y,z);
+	-- a step of (0,0,1) in (i,j,k) means a step of (-G,-G,1-G) in (x,y,z).
+	ix, iy, iz = band(ix, 255), band(iy, 255), band(iz, 255)
+
+	local n0 = GetN3(ix, iy, iz, x0, y0, z0)
+	local n3 = GetN3(ix + 1, iy + 1, iz + 1, x0 - 0.5, y0 - 0.5, z0 - 0.5) -- G3
+
+	--[[
+		Determine other corners based on simplex (skewed tetrahedron) we are in:
+
+		if x0 >= y0 then -- ~A
+			if y0 >= z0 then -- ~A and ~B
+				i1, j1, k1, i2, j2, k2 = 1, 0, 0, 1, 1, 0
+			elseif x0 >= z0 then -- ~A and B and ~C
+				i1, j1, k1, i2, j2, k2 = 1, 0, 0, 1, 0, 1
+			else -- ~A and B and C
+				i1, j1, k1, i2, j2, k2 = 0, 0, 1, 1, 0, 1
+			end
+		else -- A
+			if y0 < z0 then -- A and B
+				i1, j1, k1, i2, j2, k2 = 0, 0, 1, 0, 1, 1
+			elseif x0 < z0 then -- A and ~B and C
+				i1, j1, k1, i2, j2, k2 = 0, 1, 0, 0, 1, 1
+			else -- A and ~B and ~C
+				i1, j1, k1, i2, j2, k2 = 0, 1, 0, 1, 1, 0
+			end
+		end
+	]]
+
+	local xLy = rshift(floor(x0 - y0), 31) -- x0 < y0
+	local yLz = rshift(floor(y0 - z0), 31) -- y0 < z0
+	local xLz = rshift(floor(x0 - z0), 31) -- x0 < z0
+
+	local i1 = band(1 - xLy, bor(1 - yLz, 1 - xLz)) -- x0 >= y0 and (y0 >= z0 or x0 >= z0)
+	local j1 = band(xLy, 1 - yLz) -- x0 < y0 and y0 >= z0
+	local k1 = band(yLz, bor(xLy, xLz)) -- y0 < z0 and (x0 < y0 or x0 < z0)
+
+	local i2 = bor(1 - xLy, band(1 - yLz, 1 - xLz)) -- x0 >= y0 or (y0 >= z0 and x0 >= z0)
+	local j2 = bor(xLy, 1 - yLz) -- x0 < y0 or y0 >= z0
+	local k2 = bor(band(1 - xLy, yLz), band(xLy, bor(yLz, xLz))) -- (x0 >= y0 and y0 < z0) or (x0 < y0 and (y0 < z0 or x0 < z0))
+
+	local n1 = GetN3(ix + i1, iy + j1, iz + k1, x0 + 0.166666667 - i1, y0 + 0.166666667 - j1, z0 + 0.166666667 - k1) -- G
+	local n2 = GetN3(ix + i2, iy + j2, iz + k2, x0 + 0.333333333 - i2, y0 + 0.333333333 - j2, z0 + 0.333333333 - k2) -- G2
+
+	-- Add contributions from each corner to get the final noise value.
+	-- The result is scaled to stay just inside [-1,1]
+	return 28.452842 * (n0 + n1 + n2 + n3)
+end
+
+-- Gradients for 4D case --
+local Grads4 = ffi.new("const double[32][4]",
+	{ 0, 1, 1, 1 }, { 0, 1, 1, -1 }, { 0, 1, -1, 1 }, { 0, 1, -1, -1 },
+	{ 0, -1, 1, 1 }, { 0, -1, 1, -1 }, { 0, -1, -1, 1 }, { 0, -1, -1, -1 },
+	{ 1, 0, 1, 1 }, { 1, 0, 1, -1 }, { 1, 0, -1, 1 }, { 1, 0, -1, -1 },
+	{ -1, 0, 1, 1 }, { -1, 0, 1, -1 }, { -1, 0, -1, 1 }, { -1, 0, -1, -1 },
+	{ 1, 1, 0, 1 }, { 1, 1, 0, -1 }, { 1, -1, 0, 1 }, { 1, -1, 0, -1 },
+	{ -1, 1, 0, 1 }, { -1, 1, 0, -1 }, { -1, -1, 0, 1 }, { -1, -1, 0, -1 },
+	{ 1, 1, 1, 0 }, { 1, 1, -1, 0 }, { 1, -1, 1, 0 }, { 1, -1, -1, 0 },
+	{ -1, 1, 1, 0 }, { -1, 1, -1, 0 }, { -1, -1, 1, 0 }, { -1, -1, -1, 0 }
+)
+
+-- 4D weight contribution
+local function GetN4(ix, iy, iz, iw, x, y, z, w)
+	local t = .6 - x * x - y * y - z * z - w * w
+	local index = band(Perms[ix + Perms[iy + Perms[iz + Perms[iw]]]], 0x1F)
+
+	return max(0, (t * t) * (t * t)) * (Grads4[index][0] * x + Grads4[index][1] * y + Grads4[index][2] * z + Grads4[index][3] * w)
+end
+
+-- A lookup table to traverse the simplex around a given point in 4D.
+-- Details can be found where this table is used, in the 4D noise method.
+local Simplex = ffi.new("uint8_t[64][4]",
+	{ 0, 1, 2, 3 }, { 0, 1, 3, 2 }, {}, { 0, 2, 3, 1 }, {}, {}, {}, { 1, 2, 3 },
+	{ 0, 2, 1, 3 }, {}, { 0, 3, 1, 2 }, { 0, 3, 2, 1 }, {}, {}, {}, { 1, 3, 2 },
+	{}, {}, {}, {}, {}, {}, {}, {},
+	{ 1, 2, 0, 3 }, {}, { 1, 3, 0, 2 }, {}, {}, {}, { 2, 3, 0, 1 }, { 2, 3, 1 },
+	{ 1, 0, 2, 3 }, { 1, 0, 3, 2 }, {}, {}, {}, { 2, 0, 3, 1 }, {}, { 2, 1, 3 },
+	{}, {}, {}, {}, {}, {}, {}, {},
+	{ 2, 0, 1, 3 }, {}, {}, {}, { 3, 0, 1, 2 }, { 3, 0, 2, 1 }, {}, { 3, 1, 2 },
+	{ 2, 1, 0, 3 }, {}, {}, {}, { 3, 1, 0, 2 }, {}, { 3, 2, 0, 1 }, { 3, 2, 1 }
+)
+
+-- Convert the above indices to masks that can be shifted / anded into offsets --
+for i = 0, 63 do
+	Simplex[i][0] = lshift(1, Simplex[i][0]) - 1
+	Simplex[i][1] = lshift(1, Simplex[i][1]) - 1
+	Simplex[i][2] = lshift(1, Simplex[i][2]) - 1
+	Simplex[i][3] = lshift(1, Simplex[i][3]) - 1
+end
+
+local function simplex_4d(x, y, z, w)
+	--[[
+		4D skew factors:
+		F = (math.sqrt(5) - 1) / 4
+		G = (5 - math.sqrt(5)) / 20
+		G2 = 2 * G
+		G3 = 3 * G
+		G4 = 4 * G - 1
+	]]
+
+	-- Skew the input space to determine which simplex cell we are in.
+	local s = (x + y + z + w) * 0.309016994 -- F
+	local ix, iy, iz, iw = floor(x + s), floor(y + s), floor(z + s), floor(w + s)
+
+	-- Unskew the cell origin back to (x, y, z) space.
+	local t = (ix + iy + iz + iw) * 0.138196601 -- G
+	local x0 = x + t - ix
+	local y0 = y + t - iy
+	local z0 = z + t - iz
+	local w0 = w + t - iw
+
+	-- For the 4D case, the simplex is a 4D shape I won't even try to describe.
+	-- To find out which of the 24 possible simplices we're in, we need to
+	-- determine the magnitude ordering of x0, y0, z0 and w0.
+	-- The method below is a good way of finding the ordering of x,y,z,w and
+	-- then find the correct traversal order for the simplex we<77>re in.
+	-- First, six pair-wise comparisons are performed between each possible pair
+	-- of the four coordinates, and the results are used to add up binary bits
+	-- for an integer index.
+	local c1 = band(rshift(floor(y0 - x0), 26), 32)
+	local c2 = band(rshift(floor(z0 - x0), 27), 16)
+	local c3 = band(rshift(floor(z0 - y0), 28), 8)
+	local c4 = band(rshift(floor(w0 - x0), 29), 4)
+	local c5 = band(rshift(floor(w0 - y0), 30), 2)
+	local c6 = rshift(floor(w0 - z0), 31)
+
+	-- Simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
+	-- Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w
+	-- impossible. Only the 24 indices which have non-zero entries make any sense.
+	-- We use a thresholding to set the coordinates in turn from the largest magnitude.
+	local c = c1 + c2 + c3 + c4 + c5 + c6
+
+	-- The number 3 (i.e. bit 2) in the "simplex" array is at the position of the largest coordinate.
+	local i1 = rshift(Simplex[c][0], 2)
+	local j1 = rshift(Simplex[c][1], 2)
+	local k1 = rshift(Simplex[c][2], 2)
+	local l1 = rshift(Simplex[c][3], 2)
+
+	-- The number 2 (i.e. bit 1) in the "simplex" array is at the second largest coordinate.
+	local i2 = band(rshift(Simplex[c][0], 1), 1)
+	local j2 = band(rshift(Simplex[c][1], 1), 1)
+	local k2 = band(rshift(Simplex[c][2], 1), 1)
+	local l2 = band(rshift(Simplex[c][3], 1), 1)
+
+	-- The number 1 (i.e. bit 0) in the "simplex" array is at the second smallest coordinate.
+	local i3 = band(Simplex[c][0], 1)
+	local j3 = band(Simplex[c][1], 1)
+	local k3 = band(Simplex[c][2], 1)
+	local l3 = band(Simplex[c][3], 1)
+
+	-- Work out the hashed gradient indices of the five simplex corners
+	-- Sum up and scale the result to cover the range [-1,1]
+	ix, iy, iz, iw = band(ix, 255), band(iy, 255), band(iz, 255), band(iw, 255)
+
+	local n0 = GetN4(ix, iy, iz, iw, x0, y0, z0, w0)
+	local n1 = GetN4(ix + i1, iy + j1, iz + k1, iw + l1, x0 + 0.138196601 - i1, y0 + 0.138196601 - j1, z0 + 0.138196601 - k1, w0 + 0.138196601 - l1) -- G
+	local n2 = GetN4(ix + i2, iy + j2, iz + k2, iw + l2, x0 + 0.276393202 - i2, y0 + 0.276393202 - j2, z0 + 0.276393202 - k2, w0 + 0.276393202 - l2) -- G2
+	local n3 = GetN4(ix + i3, iy + j3, iz + k3, iw + l3, x0 + 0.414589803 - i3, y0 + 0.414589803 - j3, z0 + 0.414589803 - k3, w0 + 0.414589803 - l3) -- G3
+	local n4 = GetN4(ix + 1, iy + 1, iz + 1, iw + 1, x0 - 0.447213595, y0 - 0.447213595, z0 - 0.447213595, w0 - 0.447213595) -- G4
+
+	return 2.210600293 * (n0 + n1 + n2 + n3 + n4)
+end
+
+--- Simplex Noise
+-- @param x
+-- @param y
+-- @param z optional
+-- @param w optional
+-- @return Noise value in the range [-1, +1]
+return function(x, y, z, w)
+	if w then
+		return simplex_4d(x, y, z, w)
+	end
+	if z then
+		return simplex_3d(x, y, z)
+	end
+	if y then
+		return simplex_2d(x, y)
+	end
+	error "Simplex requires at least two arguments"
+end

libs/cpml/utils.lua

@@ -0,0 +1,228 @@
+--- Various utility functions
+-- @module utils
+
+local modules = (...): gsub('%.[^%.]+$', '') .. "."
+local vec2    = require(modules .. "vec2")
+local vec3    = require(modules .. "vec3")
+local private = require(modules .. "_private_utils")
+local abs     = math.abs
+local ceil    = math.ceil
+local floor   = math.floor
+local log     = math.log
+local utils   = {}
+
+-- reimplementation of math.frexp, due to its removal from Lua 5.3 :(
+-- courtesy of airstruck
+local log2 = log(2)
+
+local frexp = math.frexp or function(x)
+	if x == 0 then return 0, 0 end
+	local e = floor(log(abs(x)) / log2 + 1)
+	return x / 2 ^ e, e
+end
+
+--- Clamps a value within the specified range.
+-- @param value Input value
+-- @param min Minimum output value
+-- @param max Maximum output value
+-- @return number
+function utils.clamp(value, min, max)
+	return math.max(math.min(value, max), min)
+end
+
+--- Returns `value` if it is equal or greater than |`size`|, or 0.
+-- @param value
+-- @param size
+-- @return number
+function utils.deadzone(value, size)
+	return abs(value) >= size and value or 0
+end
+
+--- Check if value is equal or greater than threshold.
+-- @param value
+-- @param threshold
+-- @return boolean
+function utils.threshold(value, threshold)
+	-- I know, it barely saves any typing at all.
+	return abs(value) >= threshold
+end
+
+--- Check if value is equal or less than threshold.
+-- @param value
+-- @param threshold
+-- @return boolean
+function utils.tolerance(value, threshold)
+	-- I know, it barely saves any typing at all.
+	return abs(value) <= threshold
+end
+
+--- Scales a value from one range to another.
+-- @param value Input value
+-- @param min_in Minimum input value
+-- @param max_in Maximum input value
+-- @param min_out Minimum output value
+-- @param max_out Maximum output value
+-- @return number
+function utils.map(value, min_in, max_in, min_out, max_out)
+	return ((value) - (min_in)) * ((max_out) - (min_out)) / ((max_in) - (min_in)) + (min_out)
+end
+
+--- Linear interpolation.
+-- Performs linear interpolation between 0 and 1 when `low` < `progress` < `high`.
+-- @param low value to return when `progress` is 0
+-- @param high value to return when `progress` is 1
+-- @param progress (0-1)
+-- @return number
+function utils.lerp(low, high, progress)
+	return low * (1 - progress) + high * progress
+end
+
+--- Exponential decay
+-- @param low initial value
+-- @param high target value
+-- @param rate portion of the original value remaining per second
+-- @param dt time delta
+-- @return number
+function utils.decay(low, high, rate, dt)
+	return utils.lerp(low, high, 1.0 - math.exp(-rate * dt))
+end
+
+--- Hermite interpolation.
+-- Performs smooth Hermite interpolation between 0 and 1 when `low` < `progress` < `high`.
+-- @param progress (0-1)
+-- @param low value to return when `progress` is 0
+-- @param high value to return when `progress` is 1
+-- @return number
+function utils.smoothstep(progress, low, high)
+	local t = utils.clamp((progress - low) / (high - low), 0.0, 1.0)
+	return t * t * (3.0 - 2.0 * t)
+end
+
+--- Round number at a given precision.
+-- Truncates `value` at `precision` points after the decimal (whole number if
+-- left unspecified).
+-- @param value
+-- @param precision
+-- @return number
+utils.round = private.round
+
+--- Wrap `value` around if it exceeds `limit`.
+-- @param value
+-- @param limit
+-- @return number
+function utils.wrap(value, limit)
+	if value < 0 then
+		value = value + utils.round(((-value/limit)+1))*limit
+	end
+	return value % limit
+end
+
+--- Check if a value is a power-of-two.
+-- Returns true if a number is a valid power-of-two, otherwise false.
+-- @author undef
+-- @param value
+-- @return boolean
+function utils.is_pot(value)
+	-- found here: https://love2d.org/forums/viewtopic.php?p=182219#p182219
+	-- check if a number is a power-of-two
+	return (frexp(value)) == 0.5
+end
+
+--- Check if a value is NaN
+-- Returns true if a number is not a valid number
+-- @param value
+-- @return boolean
+utils.is_nan = private.is_nan
+
+-- Originally from vec3
+function utils.project_on(a, b)
+	local s =
+		(a.x * b.x + a.y * b.y + a.z or 0 * b.z or 0) /
+		(b.x * b.x + b.y * b.y + b.z or 0 * b.z or 0)
+
+	if a.z and b.z then
+		return vec3(
+			b.x * s,
+			b.y * s,
+			b.z * s
+		)
+	end
+
+	return vec2(
+		b.x * s,
+		b.y * s
+	)
+end
+
+-- Originally from vec3
+function utils.project_from(a, b)
+	local s =
+		(b.x * b.x + b.y * b.y + b.z or 0 * b.z or 0) /
+		(a.x * b.x + a.y * b.y + a.z or 0 * b.z or 0)
+
+	if a.z and b.z then
+		return vec3(
+			b.x * s,
+			b.y * s,
+			b.z * s
+		)
+	end
+
+	return vec2(
+		b.x * s,
+		b.y * s
+	)
+end
+
+-- Originally from vec3
+function utils.mirror_on(a, b)
+	local s =
+		(a.x * b.x + a.y * b.y + a.z or 0 * b.z or 0) /
+		(b.x * b.x + b.y * b.y + b.z or 0 * b.z or 0) * 2
+
+	if a.z and b.z then
+		return vec3(
+			b.x * s - a.x,
+			b.y * s - a.y,
+			b.z * s - a.z
+		)
+	end
+
+	return vec2(
+		b.x * s - a.x,
+		b.y * s - a.y
+	)
+end
+
+-- Originally from vec3
+function utils.reflect(i, n)
+	return i - (n * (2 * n:dot(i)))
+end
+
+-- Originally from vec3
+function utils.refract(i, n, ior)
+	local d = n:dot(i)
+	local k = 1 - ior * ior * (1 - d * d)
+
+	if k >= 0 then
+		return (i * ior) - (n * (ior * d + k ^ 0.5))
+	end
+
+	return vec3()
+end
+
+--- Get the sign of a number
+-- returns 1 for positive values, -1 for negative and 0 for zero.
+-- @param value
+-- @return number
+function utils.sign(n)
+	if n > 0 then
+		return 1
+	elseif n < 0 then
+		return -1
+	else
+		return 0
+	end
+end
+
+return utils

libs/cpml/vec2.lua

@@ -0,0 +1,444 @@
+--- A 2 component vector.
+-- @module vec2
+
+local modules = (...):gsub('%.[^%.]+$', '') .. "."
+local vec3    = require(modules .. "vec3")
+local precond = require(modules .. "_private_precond")
+local private = require(modules .. "_private_utils")
+local acos    = math.acos
+local atan2   = math.atan2 or math.atan
+local sqrt    = math.sqrt
+local cos     = math.cos
+local sin     = math.sin
+local vec2    = {}
+local vec2_mt = {}
+
+-- Private constructor.
+local function new(x, y)
+	return setmetatable({
+		x = x or 0,
+		y = y or 0
+	}, vec2_mt)
+end
+
+-- Do the check to see if JIT is enabled. If so use the optimized FFI structs.
+local status, ffi
+if type(jit) == "table" and jit.status() then
+	status, ffi = pcall(require, "ffi")
+	if status then
+		ffi.cdef "typedef struct { double x, y;} cpml_vec2;"
+		new = ffi.typeof("cpml_vec2")
+	end
+end
+
+--- Constants
+-- @table vec2
+-- @field unit_x X axis of rotation
+-- @field unit_y Y axis of rotation
+-- @field zero Empty vector
+vec2.unit_x = new(1, 0)
+vec2.unit_y = new(0, 1)
+vec2.zero   = new(0, 0)
+
+--- The public constructor.
+-- @param x Can be of three types: </br>
+-- number X component
+-- table {x, y} or {x = x, y = y}
+-- scalar to fill the vector eg. {x, x}
+-- @tparam number y Y component
+-- @treturn vec2 out
+function vec2.new(x, y)
+	-- number, number
+	if x and y then
+		precond.typeof(x, "number", "new: Wrong argument type for x")
+		precond.typeof(y, "number", "new: Wrong argument type for y")
+
+		return new(x, y)
+
+	-- {x, y} or {x=x, y=y}
+	elseif type(x) == "table" or type(x) == "cdata" then -- table in vanilla lua, cdata in luajit
+		local xx, yy = x.x or x[1], x.y or x[2]
+		precond.typeof(xx, "number", "new: Wrong argument type for x")
+		precond.typeof(yy, "number", "new: Wrong argument type for y")
+
+		return new(xx, yy)
+
+	-- number
+	elseif type(x) == "number" then
+		return new(x, x)
+	else
+		return new()
+	end
+end
+
+--- Convert point from polar to cartesian.
+-- @tparam number radius Radius of the point
+-- @tparam number theta Angle of the point (in radians)
+-- @treturn vec2 out
+function vec2.from_cartesian(radius, theta)
+	return new(radius * cos(theta), radius * sin(theta))
+end
+
+--- Clone a vector.
+-- @tparam vec2 a Vector to be cloned
+-- @treturn vec2 out
+function vec2.clone(a)
+	return new(a.x, a.y)
+end
+
+--- Add two vectors.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn vec2 out
+function vec2.add(a, b)
+	return new(
+		a.x + b.x,
+		a.y + b.y
+	)
+end
+
+--- Subtract one vector from another.
+-- Order: If a and b are positions, computes the direction and distance from b
+-- to a.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn vec2 out
+function vec2.sub(a, b)
+	return new(
+		a.x - b.x,
+		a.y - b.y
+	)
+end
+
+--- Multiply a vector by another vector.
+-- Component-size multiplication not matrix multiplication.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn vec2 out
+function vec2.mul(a, b)
+	return new(
+		a.x * b.x,
+		a.y * b.y
+	)
+end
+
+--- Divide a vector by another vector.
+-- Component-size inv multiplication. Like a non-uniform scale().
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn vec2 out
+function vec2.div(a, b)
+	return new(
+		a.x / b.x,
+		a.y / b.y
+	)
+end
+
+--- Get the normal of a vector.
+-- @tparam vec2 a Vector to normalize
+-- @treturn vec2 out
+function vec2.normalize(a)
+	if a:is_zero() then
+		return new()
+	end
+	return a:scale(1 / a:len())
+end
+
+--- Trim a vector to a given length.
+-- @tparam vec2 a Vector to be trimmed
+-- @tparam number len Length to trim the vector to
+-- @treturn vec2 out
+function vec2.trim(a, len)
+	return a:normalize():scale(math.min(a:len(), len))
+end
+
+--- Get the cross product of two vectors.
+-- Order: Positive if a is clockwise from b. Magnitude is the area spanned by
+-- the parallelograms that a and b span.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn number magnitude
+function vec2.cross(a, b)
+	return a.x * b.y - a.y * b.x
+end
+
+--- Get the dot product of two vectors.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn number dot
+function vec2.dot(a, b)
+	return a.x * b.x + a.y * b.y
+end
+
+--- Get the length of a vector.
+-- @tparam vec2 a Vector to get the length of
+-- @treturn number len
+function vec2.len(a)
+	return sqrt(a.x * a.x + a.y * a.y)
+end
+
+--- Get the squared length of a vector.
+-- @tparam vec2 a Vector to get the squared length of
+-- @treturn number len
+function vec2.len2(a)
+	return a.x * a.x + a.y * a.y
+end
+
+--- Get the distance between two vectors.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn number dist
+function vec2.dist(a, b)
+	local dx = a.x - b.x
+	local dy = a.y - b.y
+	return sqrt(dx * dx + dy * dy)
+end
+
+--- Get the squared distance between two vectors.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn number dist
+function vec2.dist2(a, b)
+	local dx = a.x - b.x
+	local dy = a.y - b.y
+	return dx * dx + dy * dy
+end
+
+--- Scale a vector by a scalar.
+-- @tparam vec2 a Left hand operand
+-- @tparam number b Right hand operand
+-- @treturn vec2 out
+function vec2.scale(a, b)
+	return new(
+		a.x * b,
+		a.y * b
+	)
+end
+
+--- Rotate a vector.
+-- @tparam vec2 a Vector to rotate
+-- @tparam number phi Angle to rotate vector by (in radians)
+-- @treturn vec2 out
+function vec2.rotate(a, phi)
+	local c = cos(phi)
+	local s = sin(phi)
+	return new(
+		c * a.x - s * a.y,
+		s * a.x + c * a.y
+	)
+end
+
+--- Get the perpendicular vector of a vector.
+-- @tparam vec2 a Vector to get perpendicular axes from
+-- @treturn vec2 out
+function vec2.perpendicular(a)
+	return new(-a.y, a.x)
+end
+
+--- Signed angle from one vector to another.
+-- Rotations from +x to +y are positive.
+-- @tparam vec2 a Vector
+-- @tparam vec2 b Vector
+-- @treturn number angle in (-pi, pi]
+function vec2.angle_to(a, b)
+	if b then
+		local angle = atan2(b.y, b.x) - atan2(a.y, a.x)
+		-- convert to (-pi, pi]
+		if angle > math.pi       then
+			angle = angle - 2 * math.pi
+		elseif angle <= -math.pi then
+			angle = angle + 2 * math.pi
+		end
+		return angle
+	end
+
+	return atan2(a.y, a.x)
+end
+
+--- Unsigned angle between two vectors.
+-- Directionless and thus commutative.
+-- @tparam vec2 a Vector
+-- @tparam vec2 b Vector
+-- @treturn number angle in [0, pi]
+function vec2.angle_between(a, b)
+	if b then
+		if vec2.is_vec2(a) then
+			return acos(a:dot(b) / (a:len() * b:len()))
+		end
+
+		return acos(vec3.dot(a, b) / (vec3.len(a) * vec3.len(b)))
+	end
+
+	return 0
+end
+
+--- Lerp between two vectors.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @tparam number s Step value
+-- @treturn vec2 out
+function vec2.lerp(a, b, s)
+	return a + (b - a) * s
+end
+
+--- Unpack a vector into individual components.
+-- @tparam vec2 a Vector to unpack
+-- @treturn number x
+-- @treturn number y
+function vec2.unpack(a)
+	return a.x, a.y
+end
+
+--- Return the component-wise minimum of two vectors.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn vec2 A vector where each component is the lesser value for that component between the two given vectors.
+function vec2.component_min(a, b)
+	return new(math.min(a.x, b.x), math.min(a.y, b.y))
+end
+
+--- Return the component-wise maximum of two vectors.
+-- @tparam vec2 a Left hand operand
+-- @tparam vec2 b Right hand operand
+-- @treturn vec2 A vector where each component is the lesser value for that component between the two given vectors.
+function vec2.component_max(a, b)
+	return new(math.max(a.x, b.x), math.max(a.y, b.y))
+end
+
+
+--- Return a boolean showing if a table is or is not a vec2.
+-- @tparam vec2 a Vector to be tested
+-- @treturn boolean is_vec2
+function vec2.is_vec2(a)
+	if type(a) == "cdata" then
+		return ffi.istype("cpml_vec2", a)
+	end
+
+	return
+		type(a)   == "table"  and
+		type(a.x) == "number" and
+		type(a.y) == "number"
+end
+
+--- Return a boolean showing if a table is or is not a zero vec2.
+-- @tparam vec2 a Vector to be tested
+-- @treturn boolean is_zero
+function vec2.is_zero(a)
+	return a.x == 0 and a.y == 0
+end
+
+--- Return whether either value is NaN
+-- @tparam vec2 a Vector to be tested
+-- @treturn boolean if x or y is nan
+function vec2.has_nan(a)
+	return private.is_nan(a.x) or
+		private.is_nan(a.y)
+end
+
+--- Convert point from cartesian to polar.
+-- @tparam vec2 a Vector to convert
+-- @treturn number radius
+-- @treturn number theta
+function vec2.to_polar(a)
+	local radius = sqrt(a.x^2 + a.y^2)
+	local theta  = atan2(a.y, a.x)
+	theta = theta > 0 and theta or theta + 2 * math.pi
+	return radius, theta
+end
+
+-- Round all components to nearest int (or other precision).
+-- @tparam vec2 a Vector to round.
+-- @tparam precision Digits after the decimal (integer if unspecified)
+-- @treturn vec2 Rounded vector
+function vec2.round(a, precision)
+	return vec2.new(private.round(a.x, precision), private.round(a.y, precision))
+end
+
+-- Negate x axis only of vector.
+-- @tparam vec2 a Vector to x-flip.
+-- @treturn vec2 x-flipped vector
+function vec2.flip_x(a)
+	return vec2.new(-a.x, a.y)
+end
+
+-- Negate y axis only of vector.
+-- @tparam vec2 a Vector to y-flip.
+-- @treturn vec2 y-flipped vector
+function vec2.flip_y(a)
+	return vec2.new(a.x, -a.y)
+end
+
+-- Convert vec2 to vec3.
+-- @tparam vec2 a Vector to convert.
+-- @tparam number the new z component, or nil for 0
+-- @treturn vec3 Converted vector
+function vec2.to_vec3(a, z)
+	return vec3(a.x, a.y, z or 0)
+end
+
+--- Return a formatted string.
+-- @tparam vec2 a Vector to be turned into a string
+-- @treturn string formatted
+function vec2.to_string(a)
+	return string.format("(%+0.3f,%+0.3f)", a.x, a.y)
+end
+
+vec2_mt.__index    = vec2
+vec2_mt.__tostring = vec2.to_string
+
+function vec2_mt.__call(_, x, y)
+	return vec2.new(x, y)
+end
+
+function vec2_mt.__unm(a)
+	return new(-a.x, -a.y)
+end
+
+function vec2_mt.__eq(a, b)
+	if not vec2.is_vec2(a) or not vec2.is_vec2(b) then
+		return false
+	end
+	return a.x == b.x and a.y == b.y
+end
+
+function vec2_mt.__add(a, b)
+	precond.assert(vec2.is_vec2(a), "__add: Wrong argument type '%s' for left hand operand. (<cpml.vec2> expected)", type(a))
+	precond.assert(vec2.is_vec2(b), "__add: Wrong argument type '%s' for right hand operand. (<cpml.vec2> expected)", type(b))
+	return a:add(b)
+end
+
+function vec2_mt.__sub(a, b)
+	precond.assert(vec2.is_vec2(a), "__add: Wrong argument type '%s' for left hand operand. (<cpml.vec2> expected)", type(a))
+	precond.assert(vec2.is_vec2(b), "__add: Wrong argument type '%s' for right hand operand. (<cpml.vec2> expected)", type(b))
+	return a:sub(b)
+end
+
+function vec2_mt.__mul(a, b)
+	precond.assert(vec2.is_vec2(a), "__mul: Wrong argument type '%s' for left hand operand. (<cpml.vec2> expected)", type(a))
+	assert(vec2.is_vec2(b) or type(b) == "number", "__mul: Wrong argument type for right hand operand. (<cpml.vec2> or <number> expected)")
+
+	if vec2.is_vec2(b) then
+		return a:mul(b)
+	end
+
+	return a:scale(b)
+end
+
+function vec2_mt.__div(a, b)
+	precond.assert(vec2.is_vec2(a), "__div: Wrong argument type '%s' for left hand operand. (<cpml.vec2> expected)", type(a))
+	assert(vec2.is_vec2(b) or type(b) == "number", "__div: Wrong argument type for right hand operand. (<cpml.vec2> or <number> expected)")
+
+	if vec2.is_vec2(b) then
+		return a:div(b)
+	end
+
+	return a:scale(1 / b)
+end
+
+if status then
+	xpcall(function() -- Allow this to silently fail; assume failure means someone messed with package.loaded
+		ffi.metatype(new, vec2_mt)
+	end, function() end)
+end
+
+return setmetatable({}, vec2_mt)

libs/cpml/vec3.lua

@@ -0,0 +1,434 @@
+--- A 3 component vector.
+-- @module vec3
+
+local modules = (...):gsub('%.[^%.]+$', '') .. "."
+local precond = require(modules .. "_private_precond")
+local private = require(modules .. "_private_utils")
+local sqrt    = math.sqrt
+local cos     = math.cos
+local sin     = math.sin
+local vec3    = {}
+local vec3_mt = {}
+
+-- Private constructor.
+local function new(x, y, z)
+	return setmetatable({
+		x = x or 0,
+		y = y or 0,
+		z = z or 0
+	}, vec3_mt)
+end
+
+-- Do the check to see if JIT is enabled. If so use the optimized FFI structs.
+local status, ffi
+if type(jit) == "table" and jit.status() then
+	status, ffi = pcall(require, "ffi")
+	if status then
+		ffi.cdef "typedef struct { double x, y, z;} cpml_vec3;"
+		new = ffi.typeof("cpml_vec3")
+	end
+end
+
+--- Constants
+-- @table vec3
+-- @field unit_x X axis of rotation
+-- @field unit_y Y axis of rotation
+-- @field unit_z Z axis of rotation
+-- @field zero Empty vector
+vec3.unit_x = new(1, 0, 0)
+vec3.unit_y = new(0, 1, 0)
+vec3.unit_z = new(0, 0, 1)
+vec3.zero   = new(0, 0, 0)
+
+--- The public constructor.
+-- @param x Can be of three types: </br>
+-- number X component
+-- table {x, y, z} or {x=x, y=y, z=z}
+-- scalar To fill the vector eg. {x, x, x}
+-- @tparam number y Y component
+-- @tparam number z Z component
+-- @treturn vec3 out
+function vec3.new(x, y, z)
+	-- number, number, number
+	if x and y and z then
+		precond.typeof(x, "number", "new: Wrong argument type for x")
+		precond.typeof(y, "number", "new: Wrong argument type for y")
+		precond.typeof(z, "number", "new: Wrong argument type for z")
+
+		return new(x, y, z)
+
+	-- {x, y, z} or {x=x, y=y, z=z}
+	elseif type(x) == "table" or type(x) == "cdata" then -- table in vanilla lua, cdata in luajit
+		local xx, yy, zz = x.x or x[1], x.y or x[2], x.z or x[3]
+		precond.typeof(xx, "number", "new: Wrong argument type for x")
+		precond.typeof(yy, "number", "new: Wrong argument type for y")
+		precond.typeof(zz, "number", "new: Wrong argument type for z")
+
+		return new(xx, yy, zz)
+
+	-- number
+	elseif type(x) == "number" then
+		return new(x, x, x)
+	else
+		return new()
+	end
+end
+
+--- Clone a vector.
+-- @tparam vec3 a Vector to be cloned
+-- @treturn vec3 out
+function vec3.clone(a)
+	return new(a.x, a.y, a.z)
+end
+
+--- Add two vectors.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn vec3 out
+function vec3.add(a, b)
+	return new(
+		a.x + b.x,
+		a.y + b.y,
+		a.z + b.z
+	)
+end
+
+--- Subtract one vector from another.
+-- Order: If a and b are positions, computes the direction and distance from b
+-- to a.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn vec3 out
+function vec3.sub(a, b)
+	return new(
+		a.x - b.x,
+		a.y - b.y,
+		a.z - b.z
+	)
+end
+
+--- Multiply a vector by another vector.
+-- Component-wise multiplication not matrix multiplication.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn vec3 out
+function vec3.mul(a, b)
+	return new(
+		a.x * b.x,
+		a.y * b.y,
+		a.z * b.z
+	)
+end
+
+--- Divide a vector by another.
+-- Component-wise inv multiplication. Like a non-uniform scale().
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn vec3 out
+function vec3.div(a, b)
+	return new(
+		a.x / b.x,
+		a.y / b.y,
+		a.z / b.z
+	)
+end
+
+--- Scale a vector to unit length (1).
+-- @tparam vec3 a vector to normalize
+-- @treturn vec3 out
+function vec3.normalize(a)
+	if a:is_zero() then
+		return new()
+	end
+	return a:scale(1 / a:len())
+end
+
+--- Scale a vector to unit length (1), and return the input length.
+-- @tparam vec3 a vector to normalize
+-- @treturn vec3 out
+-- @treturn number input vector length
+function vec3.normalize_len(a)
+	if a:is_zero() then
+		return new(), 0
+	end
+	local len = a:len()
+	return a:scale(1 / len), len
+end
+
+--- Trim a vector to a given length
+-- @tparam vec3 a vector to be trimmed
+-- @tparam number len Length to trim the vector to
+-- @treturn vec3 out
+function vec3.trim(a, len)
+	return a:normalize():scale(math.min(a:len(), len))
+end
+
+--- Get the cross product of two vectors.
+-- Resulting direction is right-hand rule normal of plane defined by a and b.
+-- Magnitude is the area spanned by the parallelograms that a and b span.
+-- Order: Direction determined by right-hand rule.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn vec3 out
+function vec3.cross(a, b)
+	return new(
+		a.y * b.z - a.z * b.y,
+		a.z * b.x - a.x * b.z,
+		a.x * b.y - a.y * b.x
+	)
+end
+
+--- Get the dot product of two vectors.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn number dot
+function vec3.dot(a, b)
+	return a.x * b.x + a.y * b.y + a.z * b.z
+end
+
+--- Get the length of a vector.
+-- @tparam vec3 a Vector to get the length of
+-- @treturn number len
+function vec3.len(a)
+	return sqrt(a.x * a.x + a.y * a.y + a.z * a.z)
+end
+
+--- Get the squared length of a vector.
+-- @tparam vec3 a Vector to get the squared length of
+-- @treturn number len
+function vec3.len2(a)
+	return a.x * a.x + a.y * a.y + a.z * a.z
+end
+
+--- Get the distance between two vectors.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn number dist
+function vec3.dist(a, b)
+	local dx = a.x - b.x
+	local dy = a.y - b.y
+	local dz = a.z - b.z
+	return sqrt(dx * dx + dy * dy + dz * dz)
+end
+
+--- Get the squared distance between two vectors.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn number dist
+function vec3.dist2(a, b)
+	local dx = a.x - b.x
+	local dy = a.y - b.y
+	local dz = a.z - b.z
+	return dx * dx + dy * dy + dz * dz
+end
+
+--- Scale a vector by a scalar.
+-- @tparam vec3 a Left hand operand
+-- @tparam number b Right hand operand
+-- @treturn vec3 out
+function vec3.scale(a, b)
+	return new(
+		a.x * b,
+		a.y * b,
+		a.z * b
+	)
+end
+
+--- Rotate vector about an axis.
+-- @tparam vec3 a Vector to rotate
+-- @tparam number phi Angle to rotate vector by (in radians)
+-- @tparam vec3 axis Axis to rotate by
+-- @treturn vec3 out
+function vec3.rotate(a, phi, axis)
+	if not vec3.is_vec3(axis) then
+		return a
+	end
+
+	local u = axis:normalize()
+	local c = cos(phi)
+	local s = sin(phi)
+
+	-- Calculate generalized rotation matrix
+	local m1 = new((c + u.x * u.x * (1 - c)),       (u.x * u.y * (1 - c) - u.z * s), (u.x * u.z * (1 - c) + u.y * s))
+	local m2 = new((u.y * u.x * (1 - c) + u.z * s), (c + u.y * u.y * (1 - c)),       (u.y * u.z * (1 - c) - u.x * s))
+	local m3 = new((u.z * u.x * (1 - c) - u.y * s), (u.z * u.y * (1 - c) + u.x * s), (c + u.z * u.z * (1 - c))      )
+
+	return new(
+		a:dot(m1),
+		a:dot(m2),
+		a:dot(m3)
+	)
+end
+
+--- Get the perpendicular vector of a vector.
+-- @tparam vec3 a Vector to get perpendicular axes from
+-- @treturn vec3 out
+function vec3.perpendicular(a)
+	return new(-a.y, a.x, 0)
+end
+
+--- Lerp between two vectors.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @tparam number s Step value
+-- @treturn vec3 out
+function vec3.lerp(a, b, s)
+	return a + (b - a) * s
+end
+
+-- Round all components to nearest int (or other precision).
+-- @tparam vec3 a Vector to round.
+-- @tparam precision Digits after the decimal (round numebr if unspecified)
+-- @treturn vec3 Rounded vector
+function vec3.round(a, precision)
+	return vec3.new(private.round(a.x, precision), private.round(a.y, precision), private.round(a.z, precision))
+end
+
+--- Unpack a vector into individual components.
+-- @tparam vec3 a Vector to unpack
+-- @treturn number x
+-- @treturn number y
+-- @treturn number z
+function vec3.unpack(a)
+	return a.x, a.y, a.z
+end
+
+--- Return the component-wise minimum of two vectors.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn vec3 A vector where each component is the lesser value for that component between the two given vectors.
+function vec3.component_min(a, b)
+	return new(math.min(a.x, b.x), math.min(a.y, b.y), math.min(a.z, b.z))
+end
+
+--- Return the component-wise maximum of two vectors.
+-- @tparam vec3 a Left hand operand
+-- @tparam vec3 b Right hand operand
+-- @treturn vec3 A vector where each component is the lesser value for that component between the two given vectors.
+function vec3.component_max(a, b)
+	return new(math.max(a.x, b.x), math.max(a.y, b.y), math.max(a.z, b.z))
+end
+
+-- Negate x axis only of vector.
+-- @tparam vec3 a Vector to x-flip.
+-- @treturn vec3 x-flipped vector
+function vec3.flip_x(a)
+	return vec3.new(-a.x, a.y, a.z)
+end
+
+-- Negate y axis only of vector.
+-- @tparam vec3 a Vector to y-flip.
+-- @treturn vec3 y-flipped vector
+function vec3.flip_y(a)
+	return vec3.new(a.x, -a.y, a.z)
+end
+
+-- Negate z axis only of vector.
+-- @tparam vec3 a Vector to z-flip.
+-- @treturn vec3 z-flipped vector
+function vec3.flip_z(a)
+	return vec3.new(a.x, a.y, -a.z)
+end
+
+function vec3.angle_to(a, b)
+	local v = a:normalize():dot(b:normalize())
+	return math.acos(v)
+end
+
+--- Return a boolean showing if a table is or is not a vec3.
+-- @tparam vec3 a Vector to be tested
+-- @treturn boolean is_vec3
+function vec3.is_vec3(a)
+	if type(a) == "cdata" then
+		return ffi.istype("cpml_vec3", a)
+	end
+
+	return
+		type(a)   == "table"  and
+		type(a.x) == "number" and
+		type(a.y) == "number" and
+		type(a.z) == "number"
+end
+
+--- Return a boolean showing if a table is or is not a zero vec3.
+-- @tparam vec3 a Vector to be tested
+-- @treturn boolean is_zero
+function vec3.is_zero(a)
+	return a.x == 0 and a.y == 0 and a.z == 0
+end
+
+--- Return whether any component is NaN
+-- @tparam vec3 a Vector to be tested
+-- @treturn boolean if x,y, or z are nan
+function vec3.has_nan(a)
+	return private.is_nan(a.x) or
+		private.is_nan(a.y) or
+		private.is_nan(a.z)
+end
+
+--- Return a formatted string.
+-- @tparam vec3 a Vector to be turned into a string
+-- @treturn string formatted
+function vec3.to_string(a)
+	return string.format("(%+0.3f,%+0.3f,%+0.3f)", a.x, a.y, a.z)
+end
+
+vec3_mt.__index    = vec3
+vec3_mt.__tostring = vec3.to_string
+
+function vec3_mt.__call(_, x, y, z)
+	return vec3.new(x, y, z)
+end
+
+function vec3_mt.__unm(a)
+	return new(-a.x, -a.y, -a.z)
+end
+
+function vec3_mt.__eq(a, b)
+	if not vec3.is_vec3(a) or not vec3.is_vec3(b) then
+		return false
+	end
+	return a.x == b.x and a.y == b.y and a.z == b.z
+end
+
+function vec3_mt.__add(a, b)
+	precond.assert(vec3.is_vec3(a), "__add: Wrong argument type '%s' for left hand operand. (<cpml.vec3> expected)", type(a))
+	precond.assert(vec3.is_vec3(b), "__add: Wrong argument type '%s' for right hand operand. (<cpml.vec3> expected)", type(b))
+	return a:add(b)
+end
+
+function vec3_mt.__sub(a, b)
+	precond.assert(vec3.is_vec3(a), "__sub: Wrong argument type '%s' for left hand operand. (<cpml.vec3> expected)", type(a))
+	precond.assert(vec3.is_vec3(b), "__sub: Wrong argument type '%s' for right hand operand. (<cpml.vec3> expected)", type(b))
+	return a:sub(b)
+end
+
+function vec3_mt.__mul(a, b)
+	precond.assert(vec3.is_vec3(a), "__mul: Wrong argument type '%s' for left hand operand. (<cpml.vec3> expected)", type(a))
+	precond.assert(vec3.is_vec3(b) or type(b) == "number", "__mul: Wrong argument type '%s' for right hand operand. (<cpml.vec3> or <number> expected)", type(b))
+
+	if vec3.is_vec3(b) then
+		return a:mul(b)
+	end
+
+	return a:scale(b)
+end
+
+function vec3_mt.__div(a, b)
+	precond.assert(vec3.is_vec3(a), "__div: Wrong argument type '%s' for left hand operand. (<cpml.vec3> expected)", type(a))
+	precond.assert(vec3.is_vec3(b) or type(b) == "number", "__div: Wrong argument type '%s' for right hand operand. (<cpml.vec3> or <number> expected)", type(b))
+
+	if vec3.is_vec3(b) then
+		return a:div(b)
+	end
+
+	return a:scale(1 / b)
+end
+
+if status then
+	xpcall(function() -- Allow this to silently fail; assume failure means someone messed with package.loaded
+		ffi.metatype(new, vec3_mt)
+	end, function() end)
+end
+
+return setmetatable({}, vec3_mt)

load/cpml.lua

@@ -0,0 +1,29 @@
+-- adapted from the init.lua CPML comes with
+
+cpml = {
+	_LICENSE = "CPML is distributed under the terms of the MIT license. See libs/cpml/LICENSE.md.",
+	_URL = "https://github.com/excessive/cpml",
+	_VERSION = "1.2.9",
+	_DESCRIPTION = "Cirno's Perfect Math Library: Just about everything you need for 3D games. Hopefully."
+}
+
+local files = {
+	"bvh",
+	"color",
+	"constants",
+	"intersect",
+	"mat4",
+	"mesh",
+	"octree",
+	"quat",
+	"simplex",
+	"utils",
+	"vec2",
+	"vec3",
+	"bound2",
+	"bound3"
+}
+
+for _, file in ipairs(files) do
+	cpml[file] = require('libs.cpml.' .. file)
+end

main.lua

@@ -9,6 +9,7 @@ function love.load()
 	require "load.bgm"
 	require "load.save"
 	require "load.bigint"
+	require 'load.cpml'
 	require "load.version"
 	loadSave()
 	require "funcs"
@@ -26,6 +27,7 @@ function love.load()
 
 	-- init config
 	initConfig()
+	config.depth_3d = 100 -- TODO add a setting for this to the menu
 
 	love.window.setFullscreen(config["fullscreen"])
 	if config.secret then playSE("welcome") end
@@ -60,8 +62,8 @@ function initModules()
 end
 
 function love.draw()
-	love.graphics.setCanvas(GLOBAL_CANVAS)
+	love.graphics.setCanvas{GLOBAL_CANVAS, depth = config.depth_3d > 0}
-	love.graphics.clear()
+	love.graphics.clear(0, 0, 0, 1, false, config.depth_3d > 0)
 
 	love.graphics.push()
 

scene/title.lua

@@ -2,6 +2,7 @@ local TitleScene = Scene:extend()
 
 TitleScene.title = "Title"
 TitleScene.restart_message = false
+local drawBlock = require "tetris.components.draw_block"
 
 local main_menu_screens = {
 	ModeSelectScene,
@@ -81,11 +82,20 @@ function TitleScene:render()
 
 	-- 490, 192
 	for _, b in ipairs(block_offsets) do
+		--[[
 		love.graphics.draw(
 			blocks["2tie"][b.color],
 			490 + b.x, 192 + b.y, 0,
 			2, 2
 		)
+		]]
+		drawBlock(
+			1, 1, 1, 1,
+			blocks["2tie"][b.color],
+			490 + b.x, 192 + b.y,
+			false, false, false, false,
+			true
+		)
 	end
 
 	--[[

tetris/components/draw_block.lua

@@ -0,0 +1,300 @@
+local mat4 = cpml.mat4
+local vec3 = cpml.vec3
+
+local cube = {
+	front = {},
+	top = {},
+	bottom = {},
+	left = {},
+	right = {}
+}
+for scale = 1,2 do
+	cube.front[scale] = love.graphics.newMesh(
+		{
+			{ "VertexPosition", "float", 3 },
+			{ "VertexTexCoord", "float", 2 },
+			{ "VertexColor", "byte", 4}
+		},
+		{
+			{
+				0, 0, 0,
+				0, 0,
+				1, 1, 1
+			},
+
+			{
+				16*scale, 0, 0,
+				1, 0,
+				1, 1, 1,
+			},
+
+			{
+				0, 16*scale, 0,
+				0, 1,
+				1, 1, 1,
+			},
+
+			{
+				16*scale, 16*scale, 0,
+				1, 1,
+				1, 1, 1,
+			},
+		},
+		"strip",
+		"static"
+	)
+	cube.top[scale] = love.graphics.newMesh(
+		{
+			{ "VertexPosition", "float", 3 },
+			{ "VertexTexCoord", "float", 2 },
+			{ "VertexColor", "byte", 4}
+		},
+		{
+			{
+				0, 0, 0,
+				0, 0,
+				1, 1, 1
+			},
+
+			{
+				16*scale, 0, 0,
+				1, 0,
+				1, 1, 1,
+			},
+
+			{
+				0, 0, 16*scale,
+				0, 1,
+				1, 1, 1,
+			},
+
+			{
+				16*scale, 0, 16*scale,
+				1, 1,
+				1, 1, 1,
+			},
+		},
+		"strip",
+		"static"
+	)
+	cube.bottom[scale] = love.graphics.newMesh(
+		{
+			{ "VertexPosition", "float", 3 },
+			{ "VertexTexCoord", "float", 2 },
+			{ "VertexColor", "byte", 4}
+		},
+		{
+			{
+				0, 16*scale, 16*scale,
+				0, 0,
+				1, 1, 1
+			},
+
+			{
+				16*scale, 16*scale, 16*scale,
+				1, 0,
+				1, 1, 1,
+			},
+
+			{
+				0, 16*scale, 0,
+				0, 1,
+				1, 1, 1,
+			},
+
+			{
+				16*scale, 16*scale, 0,
+				1, 1,
+				1, 1, 1,
+			},
+		},
+		"strip",
+		"static"
+	)
+	cube.left[scale] = love.graphics.newMesh(
+		{
+			{ "VertexPosition", "float", 3 },
+			{ "VertexTexCoord", "float", 2 },
+			{ "VertexColor", "byte", 4}
+		},
+		{
+			{
+				0, 0, 0,
+				0, 0,
+				1, 1, 1
+			},
+
+			{
+				0, 16*scale, 0,
+				1, 0,
+				1, 1, 1,
+			},
+
+			{
+				0, 0, 16*scale,
+				0, 1,
+				1, 1, 1,
+			},
+
+			{
+				0, 16*scale, 16*scale,
+				1, 1,
+				1, 1, 1,
+			},
+		},
+		"strip",
+		"static"
+	)
+	cube.right[scale] = love.graphics.newMesh(
+		{
+			{ "VertexPosition", "float", 3 },
+			{ "VertexTexCoord", "float", 2 },
+			{ "VertexColor", "byte", 4}
+		},
+		{
+			{
+				16*scale, 16*scale, 16*scale,
+				0, 0,
+				1, 1, 1
+			},
+
+			{
+				16*scale, 0, 16*scale,
+				1, 0,
+				1, 1, 1,
+			},
+
+			{
+				16*scale, 16*scale, 0,
+				0, 1,
+				1, 1, 1,
+			},
+
+			{
+				16*scale, 0, 0,
+				1, 1,
+				1, 1, 1,
+			},
+		},
+		"strip",
+		"static"
+	)
+end
+
+local shader_no_discard = {
+	shader = love.graphics.newShader(
+		[[
+		vec4 effect(vec4 color, Image tex, vec2 texture_coords, vec2 screen_coords) {
+			vec4 texcolor = Texel(tex, texture_coords);
+			return texcolor * color;
+		}
+		]],
+		[[
+		uniform mat4 transform;
+		vec4 position(mat4 transform_projection, vec4 vertex_position) {
+			return transform * TransformMatrix * (vertex_position - vec4(vec2(640.0, 480.0) / 2.0, 0.0, 0.0));
+		}
+		]]
+	),
+	width = -1,
+	height = -1,
+	depth_3d = -1
+}
+local shader_discard = {
+	shader = love.graphics.newShader(
+		[[
+		uniform vec4 rect;
+		vec4 effect(vec4 color, Image tex, vec2 texture_coords, vec2 screen_coords) {
+			if (screen_coords.x < rect.x || screen_coords.x > rect.y || screen_coords.y < rect.z || screen_coords.y > rect.w) discard;
+			vec4 texcolor = Texel(tex, texture_coords);
+			return texcolor * color;
+		}
+		]],
+		[[
+		uniform mat4 transform;
+		vec4 position(mat4 transform_projection, vec4 vertex_position) {
+			return transform * TransformMatrix * (vertex_position - vec4(vec2(640.0, 480.0) / 2.0, 0.0, 0.0));
+		}
+		]]
+	),
+	width = -1,
+	height = -1,
+	depth_3d = -1
+}
+
+local function set_shader_transform(shader, width, height)
+	local depth = (480 / 2) / math.tan(math.rad(config.depth_3d) / 2)
+	local look_at = mat4():look_at(vec3(0, 0, depth), vec3(0, 0, 0), vec3(0, 1, 0))
+	local perspective = mat4.from_perspective(config.depth_3d, width / height, -0.5, 0.5)
+	local tall = height / 480 > width / 640
+	if tall then
+		local scale_factor = (480 * width) / (640 * height)
+		shader:send("transform", mat4.scale(mat4.new(), perspective * look_at, vec3(scale_factor)))
+	else
+		shader:send("transform", perspective * look_at)
+	end
+end
+
+return function(R, G, B, A, image, x, y, left, right, top, bottom, big)
+	local scale = big and 2 or 1
+	if config.depth_3d > 0 then
+		love.graphics.setCanvas{GLOBAL_CANVAS, depth = true}
+		love.graphics.push()
+		love.graphics.origin()
+		love.graphics.setDepthMode("lequal", true)
+		local shader
+		local current_width = love.graphics.getWidth()
+		local current_height = love.graphics.getHeight()
+		if left and right and top and bottom then
+			love.graphics.setShader(shader_discard.shader)
+			local scale_factor = math.min(current_width / 640, current_height / 480)
+			local rect_left = (current_width - scale_factor * 640) / 2 + left * scale_factor
+			local rect_right = (current_width - scale_factor * 640) / 2 + right * scale_factor
+			local rect_top = (current_height - scale_factor * 480) / 2 + top * scale_factor
+			local rect_bottom = (current_height - scale_factor * 480) / 2 + bottom * scale_factor
+			shader_discard.shader:send("rect", {rect_left, rect_right, rect_top, rect_bottom})
+			shader = shader_discard
+		else
+			love.graphics.setShader(shader_no_discard.shader)
+			shader = shader_no_discard
+		end
+		if current_width ~= shader.width or current_height ~= shader.height or config.depth_3d ~= shader.depth_3d then
+			set_shader_transform(shader.shader, current_width, current_height)
+			shader.width = current_width
+			shader.height = current_height
+			shader.depth_3d = config.depth_3d
+		end
+
+		if A ~= 1 or y > 480 / 2 then
+			love.graphics.setColor(R, G, B, A)
+			cube.top[scale]:setTexture(image)
+			love.graphics.draw(cube.top[scale], x, y)
+		end
+
+		love.graphics.setColor(R / 2, G / 2, B / 2, A)
+		if A ~= 1 or y < 480 / 2 then
+			cube.bottom[scale]:setTexture(image)
+			love.graphics.draw(cube.bottom[scale], x, y)
+		end
+		if A ~= 1 or x > 640 / 2 then
+			cube.left[scale]:setTexture(image)
+			love.graphics.draw(cube.left[scale], x, y)
+		end
+		if A ~= 1 or x < 640 / 2 then
+			cube.right[scale]:setTexture(image)
+			love.graphics.draw(cube.right[scale], x, y)
+		end
+
+		love.graphics.setColor(R, G, B, A)
+		cube.front[scale]:setTexture(image)
+		love.graphics.draw(cube.front[scale], x, y)
+
+		love.graphics.setShader()
+		love.graphics.setDepthMode()
+		love.graphics.pop()
+		love.graphics.setCanvas(GLOBAL_CANVAS)
+	else
+		love.graphics.setColor(R, G, B, A)
+		love.graphics.draw(image, x, y, 0, scale, scale)
+	end
+end

tetris/components/grid.lua

@@ -2,6 +2,7 @@ local Object = require 'libs.classic'
 
 local Grid = Object:extend()
 
+local drawBlock = require 'tetris.components.draw_block'
 local empty = { skin = "", colour = "" }
 local oob = { skin = "", colour = "" }
 local block = { skin = "2tie", colour = "A" }
@@ -405,23 +406,47 @@ function Grid:update()
 end
 
 function Grid:draw()
+	is_3d = is_3d == nil and false or is_3d
 	for y = 5, self.height do
 		for x = 1, self.width do
 			if blocks[self.grid[y][x].skin] and
 			blocks[self.grid[y][x].skin][self.grid[y][x].colour] then
 				if self.grid_age[y][x] < 2 then
-					love.graphics.setColor(1, 1, 1, 1)
+					drawBlock(1, 1, 1, 1, blocks[self.grid[y][x].skin]["F"], 48+x*16, y*16,
-					love.graphics.draw(blocks[self.grid[y][x].skin]["F"], 48+x*16, y*16)
+						48, 48+(self.width+1)*16,
+						5*16, (self.height+1)*16
+					)
 				else
+					local R, G, B, A
 					if self.grid[y][x].colour == "X" then
-						love.graphics.setColor(0, 0, 0, 0)
+						R = 0
+						G = 0
+						B = 0
+						A = 0
 					elseif self.grid[y][x].skin == "bone" then
-						love.graphics.setColor(1, 1, 1, 1)
+						r = 1
+						G = 1
+						B = 1
+						A = 1
 					else
-						love.graphics.setColor(0.5, 0.5, 0.5, 1)
+						R = 0.5
+						G = 0.5
+						B = 0.5
+						A = 1
 					end
-					love.graphics.draw(blocks[self.grid[y][x].skin][self.grid[y][x].colour], 48+x*16, y*16)
+					drawBlock(R, G, B, A, blocks[self.grid[y][x].skin][self.grid[y][x].colour], 48+x*16, y*16,
+						48, 48+(self.width+1)*16,
+						5*16, (self.height+1)*16
+					)
 				end
+			end
+		end
+	end
+	-- everything that needs to be drawn over the blocks *must* be drawn after all blocks have been drawn, due to how 3D works
+	for y = 5, self.height do
+		for x = 1, self.width do
+			if blocks[self.grid[y][x].skin] and
+			blocks[self.grid[y][x].skin][self.grid[y][x].colour] then
 				if self.grid[y][x].skin ~= "bone" and self.grid[y][x].colour ~= "X" then
 					love.graphics.setColor(0.8, 0.8, 0.8, 1)
 					love.graphics.setLineWidth(1)
@@ -468,9 +493,10 @@ function Grid:drawOutline()
 	end
 end
 
-function Grid:drawInvisible(opacity_function, garbage_opacity_function, lock_flash, brightness)
+function Grid:drawInvisible(opacity_function, garbage_opacity_function, lock_flash, brightness, is_3d)
 	lock_flash = lock_flash == nil and true or lock_flash
 	brightness = brightness == nil and 0.5 or brightness
+	is_3d = is_3d == nil and false or is_3d
 	for y = 5, self.height do
 		for x = 1, self.width do
 			if self.grid[y][x] ~= empty then
@@ -481,8 +507,17 @@ function Grid:drawInvisible(opacity_function, garbage_opacity_function, lock_fla
 				else
 					opacity = opacity_function(self.grid_age[y][x])
 				end
-				love.graphics.setColor(brightness, brightness, brightness, opacity)
+				drawBlock(brightness, brightness, brightness, opacity, blocks[self.grid[y][x].skin][self.grid[y][x].colour], 48+x*16, y*16,
-				love.graphics.draw(blocks[self.grid[y][x].skin][self.grid[y][x].colour], 48+x*16, y*16)
+					48, 48+(self.width+1)*16,
+					5*16, (self.height+1)*16
+				)
+			end
+		end
+	end
+	-- everything that needs to be drawn over the blocks *must* be drawn after all blocks have been drawn, due to how 3D works
+	for y = 5, self.height do
+		for x = 1, self.width do
+			if self.grid[y][x] ~= empty then
 				if lock_flash then
 					if opacity > 0 and self.grid[y][x].colour ~= "X" then
 						love.graphics.setColor(0.64, 0.64, 0.64)
@@ -507,29 +542,40 @@ function Grid:drawInvisible(opacity_function, garbage_opacity_function, lock_fla
 	end
 end
 
-function Grid:drawCustom(colour_function, gamestate)
+function Grid:drawCustom(colour_function, gamestate, is_3d)
-    --[[
+	--[[
-        colour_function: (game, block, x, y, age) -> (R, G, B, A, outlineA)
+		colour_function: (game, block, x, y, age) -> (R, G, B, A, outlineA)
-        When called, calls the supplied function on every block passing the block itself as argument
+		When called, calls the supplied function on every block passing the block itself as argument
-        as well as coordinates and the grid_age value of the same cell.
+		as well as coordinates and the grid_age value of the same cell.
-        Should return a RGBA colour for the block, as well as the opacity of the stack outline (0 for no outline).
+		Should return a RGBA colour for the block, as well as the opacity of the stack outline (0 for no outline).
-        
+		
-        gamestate: the gamemode instance itself to pass in colour_function
+		gamestate: the gamemode instance itself to pass in colour_function
-    ]]
+	]]
+	is_3d = is_3d == nil and false or is_3d
 	for y = 5, self.height do
 		for x = 1, self.width do
-            local block = self.grid[y][x]
+			local block = self.grid[y][x]
 			if block ~= empty then
-                local R, G, B, A, outline = colour_function(gamestate, block, x, y, self.grid_age[y][x])
+				local R, G, B, A, outline = colour_function(gamestate, block, x, y, self.grid_age[y][x])
 				if self.grid[y][x].colour == "X" then
 					A = 0
 				end
-				love.graphics.setColor(R, G, B, A)
+				drawBlock(R, G, B, A, blocks[self.grid[y][x].skin][self.grid[y][x].colour], 48+x*16, y*16,
-				love.graphics.draw(blocks[self.grid[y][x].skin][self.grid[y][x].colour], 48+x*16, y*16)
+					48, 48+(self.width+1)*16,
-                if outline > 0 and self.grid[y][x].colour ~= "X" then
+					5*16, (self.height+1)*16
-                    love.graphics.setColor(0.64, 0.64, 0.64, outline)
+				)
-                    love.graphics.setLineWidth(1)
+			end
-                    if y > 5 and self.grid[y-1][x] == empty or self.grid[y-1][x].colour == "X" then
+		end
+	end
+	-- everything that needs to be drawn over the blocks *must* be drawn after all blocks have been drawn, due to how 3D works
+	for y = 5, self.height do
+		for x = 1, self.width do
+			local block = self.grid[y][x]
+			if block ~= empty then
+				if outline > 0 and self.grid[y][x].colour ~= "X" then
+					love.graphics.setColor(0.64, 0.64, 0.64, outline)
+					love.graphics.setLineWidth(1)
+					if y > 5 and self.grid[y-1][x] == empty or self.grid[y-1][x].colour == "X" then
 						love.graphics.line(48.0+x*16, -0.5+y*16, 64.0+x*16, -0.5+y*16)
 					end
 					if y < self.height and self.grid[y+1][x] == empty or
@@ -542,7 +588,7 @@ function Grid:drawCustom(colour_function, gamestate)
 					if x < self.width and self.grid[y][x+1] == empty then
 						love.graphics.line(64.5+x*16, -0.0+y*16, 64.5+x*16, 16.0+y*16)
 					end
-                end
+				end
 			end
 		end
 	end

tetris/components/piece.lua

@@ -2,6 +2,8 @@ local Object = require 'libs.classic'
 
 local Piece = Object:extend()
 
+local drawBlock = require 'tetris.components.draw_block'
+
 function Piece:new(shape, rotation, position, block_offsets, gravity, lock_delay, skin, colour, big)
 	self.shape = shape
 	self.rotation = rotation
@@ -157,7 +159,6 @@ end
 function Piece:draw(opacity, brightness, grid, partial_das)
 	if opacity == nil then opacity = 1 end
 	if brightness == nil then brightness = 1 end
-	love.graphics.setColor(brightness, brightness, brightness, opacity)
 	local offsets = self:getBlockOffsets()
 	local gravity_offset = 0
 	if config.gamesettings.smooth_movement == 1 and 
@@ -168,16 +169,24 @@ function Piece:draw(opacity, brightness, grid, partial_das)
 	for index, offset in pairs(offsets) do
 		local x = self.position.x + offset.x
 		local y = self.position.y + offset.y
-		if self.big then
+		local scale = self.big and 2 or 1
-			love.graphics.draw(
+		if grid ~= nil then
+			drawBlock(
+				brightness, brightness, brightness, opacity,
 				blocks[self.skin][self.colour],
-				64+x*32+partial_das*2, 16+y*32+gravity_offset*2,
+				64 + (x*16+partial_das)*scale, 16 + (y*16+gravity_offset)*scale,
-				0, 2, 2
+				48, 48+(grid.width+1)*16,
+				0, (grid.height+1)*16,
+				self.big
 			)
 		else
-			love.graphics.draw(
+			drawBlock(
+				brightness, brightness, brightness, opacity,
 				blocks[self.skin][self.colour],
-				64+x*16+partial_das, 16+y*16+gravity_offset
+				64 + (x*16+partial_das)*scale, 16 + (y*16+gravity_offset)*scale,
+				false, false,
+				false, false,
+				self.big
 			)
 		end
 	end

tetris/modes/gamemode.lua

@@ -5,6 +5,7 @@ local playedReadySE = false
 local playedGoSE = false
 
 local Grid = require 'tetris.components.grid'
+local drawBlock = require 'tetris.components.draw_block'
 local Randomizer = require 'tetris.randomizers.randomizer'
 local BagRandomizer = require 'tetris.randomizers.bag'
 local binser = require 'libs.binser'
@@ -787,13 +788,12 @@ function GameMode:drawLineClearAnimation()
 	for y, row in pairs(self.cleared_block_table) do
 		for x, block in pairs(row) do
 			local animation_table = self:animation(x, y, block.skin, block.colour)
-			love.graphics.setColor(
+			drawBlock(
-				animation_table[1], animation_table[2],
+				animation_table[1], animation_table[2], animation_table[3], animation_table[4],
-				animation_table[3], animation_table[4]
-			)
-			love.graphics.draw(
 				blocks[animation_table[5]][animation_table[6]],
-				animation_table[7], animation_table[8]
+				animation_table[7], animation_table[8],
+				48, 48+(self.grid.width+1)*16,
+				5*16, (self.grid.height+1)*16
 			)
 		end
 	end
@@ -820,7 +820,7 @@ function GameMode:drawGhostPiece(ruleset)
 	local ghost_piece = self.piece:withOffset({x=0, y=0})
 	ghost_piece.ghost = true
 	ghost_piece:dropToBottom(self.grid)
-	ghost_piece:draw(0.5)
+	ghost_piece:draw(0.5, 1, self.grid)
 end
 
 function GameMode:drawNextQueue(ruleset)
@@ -837,7 +837,11 @@ function GameMode:drawNextQueue(ruleset)
 		for index, offset in pairs(offsets) do
 			local x = offset.x + ruleset:getDrawOffset(piece, rotation).x + ruleset.spawn_positions[piece].x
 			local y = offset.y + ruleset:getDrawOffset(piece, rotation).y + 4.7
-			love.graphics.draw(blocks[skin][colourscheme[piece]], pos_x+x*16, pos_y+y*16)
+			drawBlock(
+				1, 1, 1, 1,
+				blocks[skin][colourscheme[piece]],
+				pos_x+x*16, pos_y+y*16
+			)
 		end
 	end
 	for i = 1, self.next_queue_length do