-- lua-tikz3dtools-occlusion.lua
--[[
We use homogeneous vector convention throughout
{{x, y, z, w, 1}}
]]
--- The dot product.
--- @param u table
> a vector
--- @param v table> another vector
--- @return number the dot product
local function dot_product(u, v)
local a, b = u[1], v[1]
return a[1]*b[1] + a[2]*b[2] + a[3]*b[3]
end
--- Scalar multiplication
local function scalar_multiplication(a, v)
local tmp = v[1]
return { {a*tmp[1], a*tmp[2], a*tmp[3], 1} }
end
local function vector_subtraction(u, v)
local U = u[1]
local V = v[1]
return { {U[1]-V[1], U[2]-V[2], U[3]-V[3], 1} }
end
local function vector_addition(u, v)
local U = u[1]
local V = v[1]
return { {U[1]+V[1], U[2]+V[2], U[3]+V[3], 1} }
end
local function distance(P1, P2)
local A = P1[1]
local B = P2[1]
return math.sqrt(
(A[1]-B[1])^2 +
(A[2]-B[2])^2 +
(A[3]-B[3])^2
)
end
--- The orthogonal vector projection.
--- @param u table> the vector being projected onto
--- @param v table> the vector being projected
--- @return table> the resulting projection
local function orthogonal_vector_projection(u, v)
return scalar_multiplication(
dot_product(v, u) / dot_product(u, u),
u
)
end
local function length(v)
local V = v[1]
return math.sqrt((V[1])^2 + (V[2])^2 + (V[3])^2)
end
local function cross_product(u,v)
local x = u[1][2]*v[1][3]-u[1][3]*v[1][2]
local y = u[1][3]*v[1][1]-u[1][1]*v[1][3]
local z = u[1][1]*v[1][2]-u[1][2]*v[1][1]
local result = { {x, y, z, 1} }
return result
end
local function sign(number)
if math.abs(number) < 0.001 then return "zero" end
if number > 0 then return "positive" end
return "negative"
end
--- Occlusion sorts two points.
--- @param P1 table> a point
--- @param P2 table> another point
--- @return boolean or nil whether P1 is in front
local function point_point_occlusion_sort(P1, P2)
--[[
If the points are not the same point, then calculate their
depth directly.
]]
if distance(P1, P2) > 0.001 then
if distance({{P1[1][1], P1[1][2], 0, 1}}, {{P2[1][1], P2[1][2], 0, 1}}) < 0.001 then
return P1[1][3] > P2[1][3]
end
end
return nil
end
-- Gauss-Jordan solver (column-major)
-- AugCols: array of n+1 columns, each column is an array of n numbers.
-- Options (optional):
-- opts.copy (default true) -> if true, solver works on a copy and does not mutate AugCols
-- opts.eps (default 1e-12) -> pivot tolerance
-- Returns: x (array 1..n) on success, or nil, reason on failure.
local function gauss_jordan_cols(AugCols, opts)
opts = opts or {}
local copy_input = (opts.copy == nil) and true or not not opts.copy
local eps = opts.eps or 1e-12
-- basic validation
local m = #AugCols
if m == 0 then return {}, nil end
local n = #AugCols[1]
if m ~= n + 1 then return nil, "bad_matrix_size (need n+1 columns for n rows)" end
for c = 1, m do
if #AugCols[c] ~= n then return nil, "bad_column_length" end
end
-- clone columns if requested
local cols = AugCols
if copy_input then
cols = {}
for c = 1, m do
local col = {}
for r = 1, n do col[r] = AugCols[c][r] end
cols[c] = col
end
end
-- Gauss-Jordan elimination
for i = 1, n do
-- find pivot row (max abs value in column i for rows i..n)
local pivot_row = i
local maxval = math.abs(cols[i][i])
for r = i+1, n do
local v = math.abs(cols[i][r])
if v > maxval then
maxval = v
pivot_row = r
end
end
if maxval < eps then
return nil, "singular"
end
-- swap rows i and pivot_row across all columns if needed
if pivot_row ~= i then
for c = 1, m do
cols[c][i], cols[c][pivot_row] = cols[c][pivot_row], cols[c][i]
end
end
-- normalize pivot row so pivot becomes 1
local pivot = cols[i][i]
for c = 1, m do
cols[c][i] = cols[c][i] / pivot
end
-- eliminate all other rows at column i
for r = 1, n do
if r ~= i then
local factor = cols[i][r]
if factor ~= 0 then
for c = 1, m do
cols[c][r] = cols[c][r] - factor * cols[c][i]
end
-- numerical cleanup
cols[i][r] = 0
end
end
end
end
-- solution is in the last column (RHS), rows 1..n
local x = {}
for i = 1, n do
x[i] = cols[m][i]
end
return x
end
--[[ Example usage:
local aug = {
{2, 1}, -- col 1 (a11, a21)
{1, -1}, -- col 2 (a12, a22)
{3, 0} -- RHS (b1, b2)
}
local sol, err = gauss_jordan_cols(aug)
-- sol -> {1, 1} (if system solvable)
]]
--- point_line_segment_occlusion_sort
--- @param P table> the point
--- @param L table> the line
--- @return true|false|nil depth relationship
local function point_line_segment_occlusion_sort(P, L)
local line_origin = { L[1] }
local line_direction_vector = vector_subtraction( { L[2] }, { L[1] } )
-- The line has affine basis {line_origin, line_direction_vector}
local point_from_the_lines_origin = vector_subtraction( P, { L[1] } )
local projection = orthogonal_vector_projection(
line_direction_vector,
point_from_the_lines_origin
)
local true_projection = vector_addition(
line_origin,
projection
)
-- We orthogonally project the point onto the line segment.
--[[
If the xy-projections of the point and true_projection are
close, then proceed. Otherwise, they do not occlude one another.
]]
local F = {{P[1][1], P[1][2], 0, 1}}
local G = {{true_projection[1][1], true_projection[1][2], 0, 1}}
if distance(F, G) < 0.001 then
local test = length(projection) / length(line_direction_vector)
if (0 0.001 and dist2 > 0.001 and dist3 > 0.001)
if (
sign_T1T2xT1P1 == sign_T2T3xT2P1 and
sign_T2T3xT2P1 == sign_T3T1xT3P1 and tmp
) then
local o = {T[1]}
local u = vector_subtraction({T[2]}, {T[1]})
local v = vector_subtraction({T[3]}, {T[1]})
local augmented_matrix = {
{u[1][1], u[1][2]} -- t
,{v[1][1], v[1][2]} -- s
,{P[1][1]-o[1][1], P[1][2]-o[1][2]} -- t+s
}
local sol = gauss_jordan_cols(augmented_matrix)
local t, s
if sol then
t = sol[1]
s = sol[2]
end
if t == nil or s == nil then return nil end
if 0 0 then
--[=[
L1A + (t)*L1_dir = L2A + (s)*L2_dir -->...
...--> [(t)*L1_dir] - [(s)*L2_dir] = [L2A - L1A]
]=]
local augmented_matrix = {
{L1_dir[1][1], L1_dir[1][2]}, --> t
{-L2_dir[1][1], -L2_dir[1][2]}, --> s
{L2A[1][1]-L1A[1][1], L2A[1][2]-L1A[1][2]} --> t+s
}
local sol = gauss_jordan_cols(augmented_matrix)
local t, s
if sol then
t = sol[1]
s = sol[2]
end
if t == nil or s == nil then return nil end
if (0 < t and t<1 and 0 < s and s < 1) then
local L1_inverse = vector_addition({L1[1]}, scalar_multiplication(t, vector_subtraction({L1[2]}, {L1[1]})))
local L2_inverse = vector_addition({L2[1]}, scalar_multiplication(s, vector_subtraction({L2[2]}, {L2[1]})))
return point_point_occlusion_sort(L1_inverse, L2_inverse)
end
else
--return point_point_occlusion_sort({L1[1]}, {L2[1]})
local M1 = point_line_segment_occlusion_sort({L1[1]}, L2)
local M2 = point_line_segment_occlusion_sort({L1[2]}, L2)
local M3 = point_line_segment_occlusion_sort({L2[1]}, L1)
if M3 ~= nil then M3 = not M3 end
local M4 = point_line_segment_occlusion_sort({L2[2]}, L1)
if M4 ~= nil then M4 = not M4 end
if (M1 == true or M2 == true or M3 == true or M4 == true) then
return true
end
if (M1 == false or M2 == false or M3 == false or M4 == false) then
return false
end
end
return nil
end
local function line_segment_triangle_occlusion_sort(L, T)
local A = point_triangle_occlusion_sort({L[1]}, T)
local B = point_triangle_occlusion_sort({L[2]}, T)
local T1, T2, T3 = {T[1], T[2]}, {T[2], T[3]}, {T[3], T[1]}
local C = line_segment_line_segment_occlusion_sort(L, T1)
local D = line_segment_line_segment_occlusion_sort(L, T2)
local E = line_segment_line_segment_occlusion_sort(L, T3)
if (
A == true or
B == true or
C == true or
D == true or
E == true
) then
return true
end
if (
A == false or
B == false or
C == false or
D == false or
E == false
) then
return false
end
return nil
end
local function triangle_triangle_occlusion_sort(T1, T2)
local T1AB, T1BC, T1CA = {T1[1],T1[2]}, {T1[2],T1[3]}, {T1[3],T1[1]}
local A = line_segment_triangle_occlusion_sort(T1AB, T2)
local B = line_segment_triangle_occlusion_sort(T1BC, T2)
local C = line_segment_triangle_occlusion_sort(T1CA, T2)
if (
A == true or
B == true or
C == true
) then
return true
end
if (
A == false or
B == false or
C == false
) then
return false
end
local D = point_triangle_occlusion_sort({T1[1]}, T2)
local E = point_triangle_occlusion_sort({T1[2]}, T2)
local F = point_triangle_occlusion_sort({T1[3]}, T2)
local G = point_triangle_occlusion_sort({T2[1]}, T1)
local H = point_triangle_occlusion_sort({T2[2]}, T1)
local I = point_triangle_occlusion_sort({T2[3]}, T1)
if (
D == true or
E == true or
F == true or
G == false or
H == false or
I == false
) then
return true
end
if (
D == false or
E == false or
F == false or
G == true or
H == true or
I == true
) then
return false
end
return nil
end
--- Occlusion sorts an arbitrary set of segments.
--- @param S1 table> a segment
--- @param S2 table> another segment
--- @return boolean or nil
local function occlusion_sort_segments(S1, S2)
if (S1.type == "point" and S2.type == "point") then
return point_point_occlusion_sort(S1.segment, S2.segment)
elseif (S1.type == "point" and S2.type == "line segment") then
return point_line_segment_occlusion_sort(S1.segment, S2.segment)
elseif (S1.type == "point" and S2.type == "triangle") then
return point_triangle_occlusion_sort(S1.segment, S2.segment)
elseif (S1.type == "line segment" and S2.type == "point") then
local ans = point_line_segment_occlusion_sort(S2.segment, S1.segment)
if ans == nil then return nil end
return not ans
elseif (S1.type == "line segment" and S2.type == "line segment") then
return line_segment_line_segment_occlusion_sort(S1.segment, S2.segment)
elseif (S1.type == "line segment" and S2.type == "triangle") then
return line_segment_triangle_occlusion_sort(S1.segment, S2.segment)
elseif (S1.type == "triangle" and S2.type == "point") then
local ans = point_triangle_occlusion_sort(S2.segment, S1.segment)
if ans == nil then return nil end
return not ans
elseif (S1.type == "triangle" and S2.type == "line segment") then
local ans = line_segment_triangle_occlusion_sort(S2.segment, S1.segment)
if ans == nil then return nil end
return not ans
elseif (S1.type == "triangle" and S2.type == "triangle") then
return triangle_triangle_occlusion_sort(S1.segment, S2.segment)
end
end
return occlusion_sort_segments