diff --git a/libs/cpml/LICENSE.md b/libs/cpml/LICENSE.md new file mode 100644 index 0000000..b9af9d9 --- /dev/null +++ b/libs/cpml/LICENSE.md @@ -0,0 +1,60 @@ +# Licenses + +CPML is Copyright (c) 2016 Colby Klein . + +CPML is Copyright (c) 2016 Landon Manning . + +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. diff --git a/libs/cpml/_private_precond.lua b/libs/cpml/_private_precond.lua new file mode 100644 index 0000000..af4106f --- /dev/null +++ b/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 diff --git a/libs/cpml/_private_utils.lua b/libs/cpml/_private_utils.lua new file mode 100644 index 0000000..fa64b02 --- /dev/null +++ b/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 diff --git a/libs/cpml/bound2.lua b/libs/cpml/bound2.lua new file mode 100644 index 0000000..2da58fe --- /dev/null +++ b/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:
+-- 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) diff --git a/libs/cpml/bound3.lua b/libs/cpml/bound3.lua new file mode 100644 index 0000000..00a7b85 --- /dev/null +++ b/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:
+-- 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) diff --git a/libs/cpml/bvh.lua b/libs/cpml/bvh.lua new file mode 100644 index 0000000..acac611 --- /dev/null +++ b/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 +}) diff --git a/libs/cpml/color.lua b/libs/cpml/color.lua new file mode 100644 index 0000000..e78e9e7 --- /dev/null +++ b/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:
+-- 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) diff --git a/libs/cpml/constants.lua b/libs/cpml/constants.lua new file mode 100644 index 0000000..9b9dab0 --- /dev/null +++ b/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 diff --git a/libs/cpml/intersect.lua b/libs/cpml/intersect.lua new file mode 100644 index 0000000..c8d7086 --- /dev/null +++ b/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 diff --git a/libs/cpml/mat4.lua b/libs/cpml/mat4.lua new file mode 100644 index 0000000..ef4f331 --- /dev/null +++ b/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:
+-- 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. ( 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. ( 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) diff --git a/libs/cpml/mesh.lua b/libs/cpml/mesh.lua new file mode 100644 index 0000000..1d25ae7 --- /dev/null +++ b/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 diff --git a/libs/cpml/octree.lua b/libs/cpml/octree.lua new file mode 100644 index 0000000..fbc15f1 --- /dev/null +++ b/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 +}) diff --git a/libs/cpml/quat.lua b/libs/cpml/quat.lua new file mode 100644 index 0000000..999f970 --- /dev/null +++ b/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:
+-- 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. ( expected)", type(a)) + precond.assert(quat.is_quat(b), "__add: Wrong argument type '%s' for right hand operand. ( 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. ( expected)", type(a)) + precond.assert(quat.is_quat(b), "__sub: Wrong argument type '%s' for right hand operand. ( 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. ( 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. ( or or 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. ( 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) diff --git a/libs/cpml/simplex.lua b/libs/cpml/simplex.lua new file mode 100644 index 0000000..74df4d9 --- /dev/null +++ b/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�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= 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 diff --git a/libs/cpml/vec2.lua b/libs/cpml/vec2.lua new file mode 100644 index 0000000..c392444 --- /dev/null +++ b/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:
+-- 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. ( expected)", type(a)) + precond.assert(vec2.is_vec2(b), "__add: Wrong argument type '%s' for right hand operand. ( 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. ( expected)", type(a)) + precond.assert(vec2.is_vec2(b), "__add: Wrong argument type '%s' for right hand operand. ( 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. ( expected)", type(a)) + assert(vec2.is_vec2(b) or type(b) == "number", "__mul: Wrong argument type for right hand operand. ( or 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. ( expected)", type(a)) + assert(vec2.is_vec2(b) or type(b) == "number", "__div: Wrong argument type for right hand operand. ( or 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) diff --git a/libs/cpml/vec3.lua b/libs/cpml/vec3.lua new file mode 100644 index 0000000..b653506 --- /dev/null +++ b/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:
+-- 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. ( expected)", type(a)) + precond.assert(vec3.is_vec3(b), "__add: Wrong argument type '%s' for right hand operand. ( 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. ( expected)", type(a)) + precond.assert(vec3.is_vec3(b), "__sub: Wrong argument type '%s' for right hand operand. ( 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. ( expected)", type(a)) + precond.assert(vec3.is_vec3(b) or type(b) == "number", "__mul: Wrong argument type '%s' for right hand operand. ( or 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. ( expected)", type(a)) + precond.assert(vec3.is_vec3(b) or type(b) == "number", "__div: Wrong argument type '%s' for right hand operand. ( or 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) diff --git a/load/cpml.lua b/load/cpml.lua new file mode 100644 index 0000000..d846fe8 --- /dev/null +++ b/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 diff --git a/main.lua b/main.lua index 52fecf9..f0e17c0 100644 --- a/main.lua +++ b/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.clear() + love.graphics.setCanvas{GLOBAL_CANVAS, depth = config.depth_3d > 0} + love.graphics.clear(0, 0, 0, 1, false, config.depth_3d > 0) love.graphics.push() diff --git a/scene/title.lua b/scene/title.lua index 35bbc82..a7973c5 100644 --- a/scene/title.lua +++ b/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 --[[ diff --git a/tetris/components/draw_block.lua b/tetris/components/draw_block.lua new file mode 100644 index 0000000..397ff91 --- /dev/null +++ b/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 diff --git a/tetris/components/grid.lua b/tetris/components/grid.lua index 2015e2b..c4659c3 100644 --- a/tetris/components/grid.lua +++ b/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) - love.graphics.draw(blocks[self.grid[y][x].skin]["F"], 48+x*16, y*16) + drawBlock(1, 1, 1, 1, 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) - love.graphics.draw(blocks[self.grid[y][x].skin][self.grid[y][x].colour], 48+x*16, y*16) + drawBlock(brightness, brightness, brightness, opacity, 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) - --[[ - 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 - 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). - - gamestate: the gamemode instance itself to pass in colour_function - ]] +function Grid:drawCustom(colour_function, gamestate, is_3d) + --[[ + 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 + 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). + + 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) - love.graphics.draw(blocks[self.grid[y][x].skin][self.grid[y][x].colour], 48+x*16, y*16) - 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 + 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 + -- 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 diff --git a/tetris/components/piece.lua b/tetris/components/piece.lua index 9be8743..059050c 100644 --- a/tetris/components/piece.lua +++ b/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 - love.graphics.draw( + local scale = self.big and 2 or 1 + 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, - 0, 2, 2 + 64 + (x*16+partial_das)*scale, 16 + (y*16+gravity_offset)*scale, + 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 diff --git a/tetris/modes/gamemode.lua b/tetris/modes/gamemode.lua index 2734806..5d38e3c 100644 --- a/tetris/modes/gamemode.lua +++ b/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( - animation_table[1], animation_table[2], - animation_table[3], animation_table[4] - ) - love.graphics.draw( + drawBlock( + animation_table[1], animation_table[2], animation_table[3], animation_table[4], 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