-- lua-tikz3dtools-implementation.lua

--- indiscernable distance
local eps = 0.0000001

--- the dot product
--- @param u table<table<number>> a vector
--- @param v table<table<number>> 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
--- @param a number the scalar coefficient
--- @param v table<table<number>> the vector
--- @return table<table<number>> the scaled vector
local function scalar_multiplication(a, v)
    local tmp = v[1]
    return { {a*tmp[1], a*tmp[2], a*tmp[3], tmp[4]} }
end

--- normalizes a vector
--- @param v table<table<number>> the vector
--- @return table<table<number>> the vectors normalization
local function normalize(v)
    local u = v[1]
    local len = math.sqrt(u[1]^2 + u[2]^2 + u[3]^2)
    return { {u[1]/len, u[2]/len, u[3]/len, u[4]} }
end

--- sign function
--- @param number number a number
--- @return string the sign|zero
local function sign(number)
    if number > 0 then return "positive" end
    if number == 0 then return "zero" end
    return "negative"
end

--- distance between two points
--- @param P1 table<table<number>> a point
--- @param P2 table<table<number>> another point
--- @return number their distance
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

--- vector subtraction
--- @param u table<table<number>> the first vector
--- @param v table<table<number>> the second vector
--- @return table<table<number>> their difference
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], U[4]} }
end

--- test for indistinguishable points
--- @param P1 table<table<number>> a point
--- @param P2 table<table<number>> another point
--- @return bool whether they intersect
local function point_point_intersecting(P1, P2)
    if distance(P1, P2) < eps then 
        return true 
    else 
        return false 
    end
end


--- CAUTION: CHATGPT GENERATED
-- Correct column-major Gauss-Jordan solver (drop-in replacement)
-- M is an array of columns; each column is an array of n rows.
-- Last column is RHS. Returns:
--   {solution = {..}, freevars = {..}}  on success
--   nil on inconsistency
--- @param M table<table<number>> an augmented matrix
--- @return table<table<number>> the solution set and free variables
local function gauss_jordan(M)
    -- basic validation
    local m = #M        -- number of columns (vars + RHS)
    if m == 0 then return {} end
    local n = #M[1]     -- number of rows (equations)
    local vars = m - 1
    if vars < 1 then return nil end
    for c = 1, m do
        if #M[c] ~= n then return nil end
    end

    -- Work on a local copy
    local cols = {}
    for c = 1, m do
        cols[c] = {}
        for r = 1, n do cols[c][r] = M[c][r] end
    end

    local rank = 0
    local row = 1
    local pivot_cols = {}  -- track pivot columns

    -- Gauss–Jordan elimination
    for col = 1, vars do
        -- find pivot row in column col, rows row..n
        local pivot_row = nil
        local maxval = eps
        for r = row, n do
            local v = math.abs(cols[col][r])
            if v > maxval then
                maxval = v
                pivot_row = r
            end
        end

        if pivot_row then
            -- swap pivot_row with current row
            if pivot_row ~= row then
                for c = 1, m do
                    cols[c][row], cols[c][pivot_row] = cols[c][pivot_row], cols[c][row]
                end
            end

            -- normalize pivot row
            local pivot = cols[col][row]
            for c = 1, m do
                cols[c][row] = cols[c][row] / pivot
            end

            -- eliminate column in other rows
            for r = 1, n do
                if r ~= row then
                    local factor = cols[col][r]
                    if math.abs(factor) > eps then
                        for c = 1, m do
                            cols[c][r] = cols[c][r] - factor * cols[c][row]
                        end
                        cols[col][r] = 0
                    end
                end
            end

            pivot_cols[#pivot_cols+1] = col
            rank = rank + 1
            row = row + 1
            if row > n then break end
        end
    end

    -- check for inconsistency: [0 0 ... | b] with b≠0
    for r = rank+1, n do
        local all_zero = true
        for c = 1, vars do
            if math.abs(cols[c][r]) > eps then
                all_zero = false
                break
            end
        end
        if all_zero and math.abs(cols[m][r]) > eps then
            return nil -- inconsistent
        end
    end

    -- Identify free variables
    local freevars = {}
    local pivotset = {}
    for _,c in ipairs(pivot_cols) do pivotset[c] = true end
    for c = 1, vars do
        if not pivotset[c] then
            freevars[#freevars+1] = c
        end
    end

    -- Extract one particular solution
    local sol = {}
    for i = 1, vars do sol[i] = 0 end
    for k,pcol in ipairs(pivot_cols) do
        sol[pcol] = cols[m][k]  -- solution from reduced form
    end

    return {solution = sol, freevars = freevars}
end

--- The orthogonal vector projection.
--- @param u table<table<number>> the vector being projected onto
--- @param v table<table<number>> the vector being projected
--- @return table<table<number>> the resulting projection
local function orthogonal_vector_projection(u, v)
    return scalar_multiplication(
        dot_product(v, u) / dot_product(u, u),
        u
    )
end

--- vector addition
--- @param u table<table<number>> the first vector
--- @param v table<table<number>> the second vector
--- @return table<table<number>> their sum
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], U[4]} }
end

--- determined whether a point and a line segment intersect
--- @param P table<table<number>> the point
--- @param L table<table<number>> the line segment
--- @return bool whether the point and line intersect
local function point_line_segment_intersecting(P, L)
    local LO = {L[1]}
    local LU = vector_subtraction({L[2]}, {L[1]})
    local LOP = vector_subtraction(P, LO)
    -- Check orthogonal projection distance
    local orth = orthogonal_vector_projection(LU, LOP)
    local real_orth = vector_addition(LO, orth)
    if distance(real_orth, P) >= eps then 
        return false 
    end
    local rhs = vector_subtraction(real_orth, LO)
    -- LO + (t) * LU = orth
    -- (t) * LU = orth - LO
    local augmented_matrix = {
        {LU[1][1], LU[1][2], LU[1][3]}
        ,{rhs[1][1], rhs[1][2], rhs[1][3]}
    }
    local sol = gauss_jordan(augmented_matrix)
    if not sol then return false end 
    if #sol.solution < 1 then return false end
    local t = sol.solution[1]
    return 0+eps<t and t<1-eps
end

--- the cross product
--- @param u table<table<number>> the first vector
--- @param v table<table<number>> the second vector
--- @return table<table<number>> their cross product
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]
    return { {x, y, z, u[1][4]} }
end

--- Projects a point onto a plane (in a given affine basis)
--- @param P table<table<number>> the point
--- @param T table<table<number>> the AFFINE basis of the plane NOT ITS TRIANGLE
--- @return table<table<number>> the literal intersection and not the coordinates
local function orthogonal_projection_onto_plane(P, T)
    local TO = {T[1]}
    local TU = {T[2]}
    local TV = {T[3]}
    local TW = cross_product(TU, TV)
    local TOP = vector_subtraction(P, TO)
    local proj_TW_TOP = orthogonal_vector_projection(TW, TOP)
    return vector_subtraction(P, proj_TW_TOP)
end

--- determines whether a point that is known to be in the triangle's plane intersects it
--- @param P table<table<number>> the point which is known to be coplannar with the triangle
--- @param T table<table<number>> the triangle
--- @return bool whether they intersect
local function point_in_triangle(P, T)
    local eps = 0.001
    local A, B, C = T[1], T[2], T[3]

    -- vectors relative to A, wrapped as table of tables
    local v0 = vector_subtraction({C}, {A})
    local v1 = vector_subtraction({B}, {A})
    local v2 = vector_subtraction(P, {A})

    -- dot products
    local d00 = dot_product(v0, v0)
    local d01 = dot_product(v0, v1)
    local d11 = dot_product(v1, v1)
    local d20 = dot_product(v2, v0)
    local d21 = dot_product(v2, v1)

    local denom = d00 * d11 - d01 * d01
    if math.abs(denom) < eps then
        return false -- degenerate triangle
    end

    -- barycentric coordinates
    local v = (d11 * d20 - d01 * d21) / denom
    local w = (d00 * d21 - d01 * d20) / denom
    local u = 1 - v - w

    -- inside triangle with epsilon buffer
    return u >= -3*eps and v >= -3*eps and w >= -3*eps
end


--- detects whether a point and a triangle intersect
--- @param P table<table<number>> the point
--- @param T table<table<number>> the triangle
--- @return bool whether they intersect
local function point_triangle_intersecting(P, T)
    local TO = {T[1]}
    local TU = vector_subtraction({T[2]}, TO)
    local TV = vector_subtraction({T[3]}, TO)
    local affine_basis = {TO[1], TU[1], TV[1]}
    local S = orthogonal_projection_onto_plane(P, affine_basis)
    return distance(S, P) < eps and point_in_triangle(P, T)
end

--- length of a vector
--- @param u table<table<number>> the vector to be measured
--- @return number the vectors length
local function length(u)
    local result = math.sqrt((u[1][1])^2 + (u[1][2])^2 + (u[1][3])^2)
    return result
end

--- return coordinates and free variables for line-line intersection
--- @param L1 table<table<number>> affine basis of a line
--- @param L1 table<table<number>> another 1D affine basis
--- @return table<table<number>> the solution and free variables
local function line_line_intersection(L1, L2)
    local L1O = {L1[1]}
    local L1U = {L1[2]}
    local L2O = {L2[1]}
    local L2U = {L2[2]}
    -- L1O + (t) * L1U = L2O + (s) * L2U
    -- (t) * L1U - (s) * L2U = L1O - L2O
    local rhs = vector_subtraction(L2O, L1O)
    local augmented_matrix = {
        {L1U[1][1], L1U[1][2], L1U[1][3]},
        {-L2U[1][1], -L2U[1][2], -L2U[1][3]},
        {rhs[1][1], rhs[1][2], rhs[1][3]}
    }
    return gauss_jordan(augmented_matrix)
end

--- determines the solution coefficients and solution set of two line segments, or produces nil
--- @param L1 table<table<number>> the first line segment
--- @param L2 table<table<number>> the second line segment
--- @return table<table<table<number>>,table<table<number>>>|nil the solution coefficients and solution set
local function line_segment_line_segment_intersection(L1, L2)
    local num = 0
    for _, P1 in ipairs(L1) do 
        for _, P2 in ipairs(L2) do 
            if distance({P1}, {P2}) < eps then 
                num = num + 1
            end
        end
    end
    if num ~= 0 then return nil end
    local L1O = {L1[1]}
    local L1U = vector_subtraction({L1[2]}, L1O)
    local L1A = {L1O[1], L1U[1]}
    local L2O = {L2[1]}
    local L2U = vector_subtraction({L2[2]}, L2O)
    local L2A = {L2O[1], L2U[1]}
    local coeffs = line_line_intersection(L1A, L2A)
    if coeffs == nil then return nil end
    if #coeffs.solution == 0 then return nil end
    local t, s = coeffs.solution[1], coeffs.solution[2]
    if 0 < t - eps and t < 1 - eps and 0 < s - eps and s < 1 - eps then
        return {
            coefficients = coeffs,
            intersection = vector_addition(
                L1O,
                scalar_multiplication(t, L1U)
            )
        }
    end
    return nil
end

--- obtains the coordinates and free variables of the intersection 
--- between a line and a plane, each defined by their affine bases.
--- @param L table<table<number>> an affine basis of a line
--- @param T table<table<number>> an affine basis of a plane
--- @return table<table<number>> the solution and free variables
local function line_plane_intersection(L, T)
    local LO = {L[1]}
    local LU = {L[2]}
    local TO = {T[1]}
    local TU = {T[2]}
    local TV = {T[3]}
    --- LO + (t) * LU = TO + (s) * TU + (w) * TV
    -- (t) * LU - (s) * TU - (w) * TV = TO - LO
    local rhs = vector_subtraction(TO, LO)

    local augmented_matrix = {
        {LU[1][1], LU[1][2], LU[1][3]},
        {-TU[1][1], -TU[1][2], -TU[1][3]},
        {-TV[1][1], -TV[1][2], -TV[1][3]},
        {rhs[1][1], rhs[1][2], rhs[1][3]}
    }
    local sol = gauss_jordan(augmented_matrix)
    if not sol then
        -- print("No solution: line-plane system inconsistent")
        return nil
    end
    if #sol.freevars > 0 then
        -- print("Line lies in the plane (coplanar case). Free vars:", #sol.freevars)
    end
    
    return sol
end

--- determines the intersection coefficients and intersection point of a line segment and triangle, if it exists, or returns nil
--- @param L table<table<number>> line segment defined by endpoints
--- @param T table<table<number>> triangle defined by vertices
--- @return table<table<table<number>>,table<table<number>>>|nil the solution coefficients and literal R3 intersection point
local function line_segment_triangle_intersection(L, T)
    local eps = 0.0000001
    local num = 0
    for _, P1 in ipairs(L) do 
        for _, P2 in ipairs(T) do 
            if distance({P1}, {P2}) < eps then 
                num = num + 1
            end
        end
    end
    if num > 1 then return nil end
    local LO = {L[1]}
    local LU = vector_subtraction({L[2]}, LO)
    local LA = {LO[1], LU[1]}
    local TO = {T[1]}
    local TU = vector_subtraction({T[2]}, TO)
    local TUA = {TO[1], TU[1]}
    local TV = vector_subtraction({T[3]}, TO)
    local TVA = {TO[1], TV[1]}
    local TA = {TO[1], TU[1], TV[1]}
    local TUVA = {
        vector_addition(TO, TU)[1]
        ,vector_subtraction(TV, TU)[1]
    }
    local coeffs = line_plane_intersection(LA, TA)
    if coeffs == nil then return nil end
    if #coeffs.freevars == 0  then 
        local t = coeffs.solution[1]
        if 0+eps<=t and t<=1-eps then
            local I = vector_addition(
                LO,
                scalar_multiplication(t, LU)
            )
            if point_in_triangle(I, T) then 
                return {
                    solution = coeffs,
                    intersection = I
                }
            end
        end
    end
    return nil
end

  local function fmt_vec(v)
    if type(v) ~= "table" then return tostring(v) end
    local parts = {}
    for i = 1, #v do parts[#parts+1] = tostring(v[i]) end
    return "{" .. table.concat(parts, ", ") .. "}"
  end

local function fmt_matrix(m)
if type(m) ~= "table" then return tostring(m) end
local rows = {}
for r = 1, #m do rows[#rows+1] = fmt_vec(m[r]) end
return "[" .. table.concat(rows, " | ") .. "]"
end

--- produces exactly zero, or exactly two intersection points between two triangles
--- @param T1 table<table<number>> a triangle defined by its vertices
--- @param T2 table<table<number>> another triangle defined by its vertices
--- @return table<table<number>> a matrix of exactly zero or exactly two points nothing else
local function triangle_triangle_intersections(T1, T2)
    local edges1 = { {T1[1], T1[2]}, {T1[2], T1[3]}, {T1[3], T1[1]} }
    local edges2 = { {T2[1], T2[2]}, {T2[2], T2[3]}, {T2[3], T2[1]} }

    local points = {}

    --- appends a point to a list if it is unique
    --- @param P table<table<number>> a point
    local function add_unique(P)
        -- if not P then return nil end
        for _, Q in ipairs(points) do
            if point_point_intersecting(P, Q) then return nil end
        end
        table.insert(points, P)
    end

    -- Check all edge pairs
    for _, E1 in ipairs(edges1) do
        local intersect = line_segment_triangle_intersection(E1, T2)
        if intersect ~= nil then
            add_unique(intersect.intersection)
        end
    end
    for _, E2 in ipairs(edges2) do
        local intersect = line_segment_triangle_intersection(E2, T1)
        if intersect ~= nil then
            add_unique(intersect.intersection)
        end
    end

    if #points ~= 2 then
        -- print(("two triangles intersected at %f points"):format(#points))
        -- print(fmt_matrix(T1))
        -- print(fmt_matrix(T2))
        return nil 
        --assert(false, ("two triangles intersected at %f points"):format(#points)) 
    end
    return points
end

--- returns a table of line segments - the partition of L into two - if L itnersects a point, else returns nil
--- @param P table<table<number>> the point
--- @param L table<table<number>> the line segment
--- @return table<table<table<number>>,table<table<number>>>|nil the intersection points if they exist
local function point_line_segment_partition(P, L)
    if point_line_segment_intersecting(P, L) then 
        return {
            point1 = {L[1], P[1]}, 
            point2 = {L[2], P[1]}
        }
    end 
    return nil
end

--- if two line segments intersect, returns the pieces of the first 
--- one, after being partitioned by the intersection with the other
--- @param L1 table<table<number>> the first line segment
--- @param L2 table<table<number>> the second line segment
--- @return table<table<table<number>>,table<table<number>>> the fragmented pieces of the first segment
local function line_segment_line_segment_partition(L1, L2)
    local eps = 0.0000001
    local intersect = line_segment_line_segment_intersection(L1, L2)
    if intersect == nil then return nil end
    local I = intersect.intersection
    -- if distance(I, L1) < eps then 
    --     return {
    --         line_segment1 = {L2[1], I[1]},
    --         line_segment2 = line_segment1
    --     }
    -- end
    -- if distance(I, L2) < eps then 
    --     return {
    --         line_segment1 = {L1[1], I[1]},
    --         line_segment2 = line_segment1
    --     }
    -- end
    return {
        line_segment1 = {L1[1], I[1]},
        line_segment2 = {L1[2], I[1]}
    }
end

--- partitions a line segment by its intersection with a triangle, or returns nil
--- @param L table<table<number>> the line segment
--- @param T table<table<number>> the triangle
--- @return table<table<table<number>>,table<table<number>>>|nil the partitioned line segment or nil
local function line_segment_triangle_clip(L, T)
    local intersect = line_segment_triangle_intersection(L, T)
    if intersect == nil then return nil end
    local I = intersect.intersection
    return {
        line_segment1 = {L[1], I[1]},
        line_segment2 = {L[2], I[1]}
    }
end

--- ChatGPT written centroid sorting function, for the purpose of convex polygon ordering
--- @param points table<table<number>> the not necessarily convex matrix of polygon points
--- @return table<table<number>> the definitively convex matrix of polygon points
local function centroid_sort(points)
    local num = #points
    assert(num >= 3, "Need at least 3 points to sort by centroid.")

    -- Compute centroid
    local sum = { {0, 0, 0, 1} }
    for i = 1, num do
        sum = vector_addition(sum, {points[i]})
    end
    local centroid = scalar_multiplication(1/num, sum)

    -- Build local basis for projection
    local v1 = vector_subtraction({points[2]}, {points[1]})
    local v2 = vector_subtraction({points[3]}, {points[1]})
    local normal = cross_product(v1, v2)

    -- Pick axes in the plane
    local x_axis = v1
    local y_axis = cross_product(normal, x_axis)

    x_axis = normalize(x_axis)
    y_axis = normalize(y_axis)

    -- Compute angle of each point relative to centroid
    local angles = {}
    for i, P in ipairs(points) do
        local rel = vector_subtraction({P}, centroid)[1]
        local x = rel[1]*x_axis[1][1] + rel[2]*x_axis[1][2] + rel[3]*x_axis[1][3]
        local y = rel[1]*y_axis[1][1] + rel[2]*y_axis[1][2] + rel[3]*y_axis[1][3]
        local theta = math.atan2(y, x)
        table.insert(angles, {point = P, angle = theta})
    end

    -- Sort by angle
    table.sort(angles, function(a, b) return a.angle < b.angle end)

    -- Extract sorted points
    local sorted = {}
    for _, entry in ipairs(angles) do
        table.insert(sorted, entry.point)
    end

    return sorted
end

--- returns the partitio of the first triangle into three subtriangles, 
--- if it intersects the second, otherwise produces nil
--- @param T1 table<table<number>> the first triangle
--- @param T2 table<table<number>> the second triangle
--- @return table<table<table<number>>,table<table<number>>,table<table<number>>> table of sub triangles
local function triangle_triangle_partition(T1, T2)
    local I = triangle_triangle_intersections(T1, T2)
    if I == nil then return nil end
    if #I == 0 then return nil end
    if #I == 1 then return nil end
    if #I ~= 2 then  assert(false, ("I is not 2, it is instead: %f"):format(#I)) end
    
    local IO = I[1]
    local IU = vector_subtraction(I[2], IO)
    local I_basis = {IO[1], IU[1]}
    local T1A = {T1[1], T1[2]}
    local T1AU = vector_subtraction({T1A[2]}, {T1A[1]})
    local T1A_basis = {T1A[1], T1AU[1]}
    local T1B = {T1[2], T1[3]}
    local T1BU = vector_subtraction({T1B[2]}, {T1B[1]})
    local T1B_basis = {T1B[1], T1BU[1]}
    local T1C = {T1[3], T1[1]}
    local T1CU = vector_subtraction({T1C[2]}, {T1C[1]})
    local T1C_basis = {T1C[1], T1CU[1]}
    local T2A = {T2[1], T2[2]}
    local T2AU = vector_subtraction({T2A[2]}, {T2A[1]})
    local T2A_basis = {T2A[1], T2AU[1]}
    local T2B = {T2[2], T2[3]}
    local T2BU = vector_subtraction({T2B[2]}, {T2B[1]})
    local T2B_basis = {T2B[1], T2BU[1]}
    local T2C = {T2[3], T2[1]}
    local T2CU = vector_subtraction({T2C[2]}, {T2C[1]})
    local T2C_basis = {T2C[1], T2CU[1]}

    local points = {}
    local non_intersecting = nil

    local int1 = line_line_intersection(I_basis, T1A_basis)
    if int1 == nil then 
        int1 = {solution = {}}
    end
    if #int1.solution ~= 0 then 
        local t = int1.solution[1]
        local intersect = vector_addition(
            IO,
            scalar_multiplication(t, IU)
        )
        if point_line_segment_intersecting(intersect, T1A) then
            table.insert(points, intersect)
        else 
            non_intersecting = "T1A" 
        end
    else 
        non_intersecting = "T1A"
    end

    local int2 = line_line_intersection(I_basis, T1B_basis)
    if int2 == nil then 
        int2 = {solution = {}}
    end
    if #int2.solution ~= 0 then 
        local t = int2.solution[1]
        local intersect = vector_addition(
            IO,
            scalar_multiplication(t, IU)
        )
        if point_line_segment_intersecting(intersect, T1B) then
            table.insert(points, intersect)
        else 
            non_intersecting = "T1B" 
        end
    else 
        non_intersecting = "T1B"
    end

    local int3 = line_line_intersection(I_basis, T1C_basis)
    if int3 == nil then 
        int3 = {solution = {}}
    end
    if #int3.solution ~= 0 then 
        local t = int3.solution[1]
        local intersect = vector_addition(
            IO,
            scalar_multiplication(t, IU)
        )
        if point_line_segment_intersecting(intersect, T1C) then
            table.insert(points, intersect)
        else 
            non_intersecting = "T1C" 
        end
    else 
        non_intersecting = "T1C"
    end

    -- if #points == 3 then return nil end 
    -- if #points == 1 then return nil end 
    if #points ~= 2 then
        -- print("Partition failure: got", #points, "points")
        -- print("Triangle 1:", T1)
        -- print("Triangle 2:", T2)
        -- print("Non-intersecting edge:", non_intersecting)
        -- for i, p in ipairs(points) do
        --     print("Point", i, p[1][1], p[1][2], p[1][3])
        -- end
        return nil
        --assert(false, "triangle_triangle_partition doesn't have exactly two points")
    end
    
    local quad = {}
    local tri1
    local A, B = points[1], points[2]
    table.insert(quad, A[1])
    table.insert(quad, B[1])
    if non_intersecting == "T1A" then 
        table.insert(quad, T1A[1])
        table.insert(quad, T1A[2])
        tri1 = {A[1], B[1], T1B[2]}
    elseif non_intersecting == "T1B" then 
        table.insert(quad, T1B[1])
        table.insert(quad, T1B[2])
        tri1 = {A[1], B[1], T1C[2]}
    elseif non_intersecting == "T1C" then
        table.insert(quad, T1C[1])
        table.insert(quad, T1C[2])
        tri1 = {A[1], B[1], T1A[2]}
    end
    
    quad = centroid_sort(quad)
    return {
        tri1 = tri1,
        tri2 = {quad[1], quad[2], quad[3]},
        tri3 = {quad[3], quad[4], quad[1]}
    }
end
























-- https://tex.stackexchange.com/a/747040
--- Creates a TeX command that evaluates a Lua function
---
--- @param name string The name of the `\csname` to define
--- @param func function
--- @param args table<string> The TeX types of the function arguments
--- @param protected boolean|nil Define the command as `\protected`
--- @return nil
local function register_tex_cmd(name, func, args, protected)
    -- The extended version of this function uses `N` and `w` where appropriate,
    -- but only using `n` is good enough for exposition purposes.
    name = "__lua_tikztdtools_" .. name .. ":" .. ("n"):rep(#args)

    -- Push the appropriate scanner functions onto the scanning stack.
    local scanners = {}
    for _, arg in ipairs(args) do
        scanners[#scanners+1] = token['scan_' .. arg]
    end

    -- An intermediate function that properly "scans" for its arguments
    -- in the TeX side.
    local scanning_func = function()
        local values = {}
        for _, scanner in ipairs(scanners) do
            values[#values+1] = scanner()
        end

        func(table.unpack(values))
    end

    local index = luatexbase.new_luafunction(name)
    lua.get_functions_table()[index] = scanning_func

    if protected then
        token.set_lua(name, index, "protected")
    else
        token.set_lua(name, index)
    end
end



local segments = {}









































-- lua-tikz3dtools-matrix.lua

local mm = {}
mm.tau = 2*math.pi 
local pi = math.pi



local cos, sin = math.cos, math.sin

--- matrix multiplication
---
--- @param A table<table<number>> left matrix
--- @param B table<table<number>> right matrix
--- @return table<table<number>> their product
function mm.matrix_multiply(A, B)
    local rows_A = #A
    local columns_A = #A[1]
    local rows_B = #B
    local columns_B = #B[1]
    assert(
        columns_A == rows_B,
        string.format(
            [[
                Wrong size matrices for multiplication.
                Size A: %d,%d Size B: %d,%d
            ]],
            rows_A, columns_A,
            rows_B, columns_B
        )
    )
    local product = {}
    for row = 1, rows_A do
        product[row] = {}
        for column = 1, columns_B do
            product[row][column] = 0
            for dot_product_step = 1, columns_A do
                local a = A[row][dot_product_step]
                local b = B[dot_product_step][column]
                assert(type(a) == "number", 
                    string.format("Expected number but got %s in A[%d][%d]", type(a), row, dot_product_step))
                assert(type(b) == "number", 
                    string.format("Expected number but got %s in B[%d][%d]", type(b), dot_product_step, column))
                product[row][column] = product[row][column] + a * b
            end
        end
    end
    return product
end

function mm.reciprocate_by_homogenous(vector)
    local result = {}
    for i = 1, #vector do
        local row = vector[i]
        local w = row[4]
        if w == 0 then
            error("Cannot reciprocate row " .. i .. ": homogeneous coordinate w = 0")
        end
        --if w<0 then w=-w end
        result[i] = {
            row[1]/w,
            row[2]/w,
            row[3]/w,
            1
        }
    end
    return result
end

function  mm.transpose(A)
    local rows_A = #A
    local columns_A = #A[1]
    local result = {}
    for row = 1, columns_A, 1 do
        result[row] = {}
        for column = 1, rows_A, 1 do
            result[row][column] = A[column][row]
        end
    end
    return result
end

function mm.matrix_inverse(matrix)
    local rows = #matrix
    local columns = #matrix[1]
    assert(rows == columns, "You can only take the inverse of a square matrix.")
    local det = mm.det(matrix)
    assert(math.abs(math.abs(det)) > 0.00001, "You cannot take the inverse of a singular matrix.")

    local n = rows
    -- Build an augmented matrix [A | I]
    local augment = {}
    for i = 1, n do
        augment[i] = {}
        -- copy row i of A
        for j = 1, n do
            augment[i][j] = matrix[i][j]
        end
        -- append row i of I
        for j = 1, n do
            augment[i][n + j] = (i == j) and 1 or 0
        end
    end

    -- Gauss-Jordan elimination
    for i = 1, n do
        -- If pivot is zero (or very close), swap with a lower row that has a nonzero pivot
        if math.abs(augment[i][i]) < 1e-12 then
            local swapRow = nil
            for r = i + 1, n do
                if math.abs(augment[r][i]) > 1e-12 then
                    swapRow = r
                    break
                end
            end
            assert(swapRow, "Matrix is singular (zero pivot encountered).")
            augment[i], augment[swapRow] = augment[swapRow], augment[i]
        end

        -- Normalize row i so that augment[i][i] == 1
        local pivot = augment[i][i]
        for col = 1, 2 * n do
            augment[i][col] = augment[i][col] / pivot
        end

        -- Eliminate column i in all other rows
        for r = 1, n do
            if r ~= i then
                local factor = augment[r][i]
                for col = 1, 2 * n do
                    augment[r][col] = augment[r][col] - factor * augment[i][col]
                end
            end
        end
    end

    -- Extract the inverse matrix from the augmented result
    local inv = {}
    for i = 1, n do
        inv[i] = {}
        for j = 1, n do
            inv[i][j] = augment[i][n + j]
        end
    end

    return inv
end




function mm.det(matrix)
    local rows = #matrix
    local columns = #matrix[1]
    assert(rows > 0, "Matrix must have at least one row to take determinant.")
    assert(columns > 0, "Matrix must have at least one column to take determinant.")
    assert(rows == columns, "You can only take the determinant of a square matrix.")
    if rows == 1 then
        return matrix[1][1]
    elseif rows == 2 then
        -- return a*d - b*c
        return matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]
    end
    -- We will do a cofactor expansion on the first row.
    local det = 0
    local minor
    local new_row
    for element = 1, columns, 1 do
        minor = {}
        for row = 2, rows, 1 do
            new_row = {}
            for column = 1, columns, 1 do
                if column ~= element then
                    table.insert(new_row, matrix[row][column])
                end
            end
            table.insert(minor,new_row)
        end
        det = det + matrix[1][element] * (-1)^(element+1) * mm.det(minor)
    end
    return det
end

function mm.yrotation(angle)
    local c = cos(angle)
    local s = sin(angle)
    return {
        {c,0,-s,0}
        ,{0,1,0,0}
        ,{s,0,c,0}
        ,{0,0,0,1}
    }
end

function mm.translate(x,y,z)
    return {
        {1,0,0,0}
        ,{0,1,0,0}
        ,{0,0,1,0}
        ,{x,y,z,1}
    }
end

function mm.xscale(scale)
    return {
        {scale,0,0,0}
        ,{0,1,0,0}
        ,{0,0,1,0}
        ,{0,0,0,1}
    }
end

function mm.yscale(scale)
    return {
        {1,0,0,0}
        ,{0,scale,0,0}
        ,{0,0,1,0}
        ,{0,0,0,1}
    }
end

function mm.zscale(scale)
    return {
        {1,0,0,0}
        ,{0,1,0,0}
        ,{0,0,scale,0}
        ,{0,0,0,1}
    }
end

function mm.scale(scale)
    return {
        {scale,0,0,0}
        ,{0,scale,0,0}
        ,{0,0,scale,0}
        ,{0,0,0,1}
    }
end

function mm.xrotation(angle)
    return {
        {1,0,0,0}
        ,{0,math.cos(angle),math.sin(angle),0}
        ,{0,-math.sin(angle),math.cos(angle),0}
        ,{0,0,0,1}
    }
end


function mm.zrotation(angle)
    local c = cos(angle)
    local s = sin(angle)
    return {
        {c,s,0,0}
        ,{-s,c,0,0}
        ,{0,0,1,0}
        ,{0,0,0,1}
    }
end

function mm.euler(alpha,beta,gamma)
    return mm.matrix_multiply(
        mm.zrotation(gamma)
        ,mm.matrix_multiply(
            mm.yrotation(beta)
            ,mm.zrotation(alpha)
        )
    )
end

function mm.sphere(longitude,latitude)
    local s = sin(latitude)
    return {{
        s * cos(longitude)
        ,s * sin(longitude)
        ,cos(latitude)
        ,1
    }}
end

function mm.matrix_add(A, B)
    local rows_A = #A
    local columns_A = #A[1]
    local rows_B = #B
    local columns_B = #B[1]
    assert(rows_A == rows_B and columns_A == columns_B, "Wrong size matrices for addition.")
    local sum = {}
    for row = 1, rows_A, 1 do
        sum[row] = {}
        for column = 1, columns_A, 1 do
            sum[row][column] = A[row][column] + B[row][column]
        end
    end
    return sum
end

function mm.matrix_subtract(A,B)
    local rows_A = #A
    local columns_A = #A[1]
    local rows_B = #B
    local columns_B = #B[1]
    assert(rows_A == rows_B and columns_A == columns_B, "Wrong size matrices for subtraction.")
    local sum = {}
    for row = 1, rows_A, 1 do
        sum[row] = {}
        for column = 1, columns_A, 1 do
            sum[row][column] = A[row][column] - B[row][column]
        end
    end
    return sum
end

function mm.matrix_scale(factor, A)
    local rows = #A
    local cols = #A[1]
    local result = {}
    for i = 1, rows do
        result[i] = {}
        for j = 1, cols do
            result[i][j] = A[i][j] * factor
        end
    end
    return result
end

function mm.sign(number)
    if number >= 0 then return "positive" end
    return "negative"
end

function mm.dot_product(u,v)
    local result = u[1][1]*v[1][1] + u[1][2]*v[1][2] + u[1][3]*v[1][3]
    return result
end

function mm.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,u[1][4]}}
    return result
end

function mm.norm(u)
    local result = math.sqrt((u[1][1])^2 + (u[1][2])^2 + (u[1][3])^2)
    return result
end

function mm.normalize(vector)
    local len = mm.norm(vector)
    return {{
        vector[1][1]/len
        ,vector[1][2]/len
        ,vector[1][3]/len
        ,1
    }}
end

function mm.identity_matrix()
    local I = {}
    for i = 1, 4 do
        I[i] = {}
        for j = 1, 4 do
            I[i][j] = (i == j) and 1 or 0
        end
    end
    return I
end

function mm.midpoint(triangle)
    local P,Q,R = table.unpack(triangle)
    local x = (P[1]+Q[1]+R[1])/3
    local y = (P[2]+Q[2]+R[2])/3
    local z = (P[3]+Q[3]+R[3])/3
    return {{x,y,z,1}}
end

function mm.orthogonal_vector(u)
    local v
    if (math.abs(u[1][1])>0.001 and math.abs(u[1][2])<0.001 and math.abs(u[1][3])<0.001) then
        v = mm.cross_product(u,{{0,1,0,1}})
    else
        v = mm.cross_product(u,{{1,0,0,1}})
    end
    return v
end

function mm.get_observer_plane_basis(observer)
    local origin = {{0,0,0,1}}
    local basis_i = mm.orthogonal_vector(observer)
    basis_i = mm.normalize(basis_i)
    local basis_j = mm.cross_product(observer,basis_i)
    basis_j = mm.normalize(basis_j)
    return {origin,basis_i,basis_j}
end

function mm.orthogonal_vector_projection(base_vector,projected_vector)
    local scale = (
        mm.dot_product(base_vector,projected_vector) / 
        mm.dot_product(base_vector,base_vector)
    )
    return {{base_vector[1][1]*scale,base_vector[1][2]*scale,base_vector[1][3]*scale,1}}
end

function mm.project_point_onto_basis(point,basis)
    local normal = mm.cross_product(basis[2],basis[3])
    normal = mm.normalize(normal)
    local vector_from_plane = mm.orthogonal_vector_projection(point,normal)
    local result = {{
        point[1][1]-vector_from_plane[1][1]
        ,point[1][2]-vector_from_plane[1][2]
        ,point[1][3]-vector_from_plane[1][3]
        ,1
    }}
    return result
end

function mm.transform_about(triplet_point,f)
    local p = triplet_point
    local a = triplet_angles
    return mm.matrix_multiply(
        mm.matrix_inverse(mm.translate(p[1],p[2],p[3]))
        ,mm.matrix_multiply(f,mm.translate(p[1],p[2],p[3]))
    )
end


function mm.scale3(x,y,z)
    return {
        {x,0,0,0}
        ,{0,y,0,0}
        ,{0,0,z,0}
        ,{0,0,0,1}
    }
end








































































local lua_tikz3dtools_parametric = {}
lua_tikz3dtools_parametric.math = {}
for k, v in pairs(_G) do lua_tikz3dtools_parametric.math[k] = v end
for k, v in pairs(mm) do lua_tikz3dtools_parametric.math[k] = v end
for k, v in pairs(math) do lua_tikz3dtools_parametric.math[k] = v end

local function single_string_expression(str)
    local func = load(([[return %s]]):format(str), "expression", "t", lua_tikz3dtools_parametric.math)()
    -- if not func then
    --     error("Failed to parse expression: " .. expr)
    --     return
    -- end
    -- local ok, result = pcall(func)
    -- if not ok then
    --     error("Error evaluating expression: " .. result)
    --     return
    -- end
    -- if type(result) ~= "function" then
    --     error("Expected a function, got: " .. type(result))
    --     return
    -- end

    return func--result
end

local function single_string_function(str)
    return load(("return function(u) return %s end"):format(str), "expression", "t", lua_tikz3dtools_parametric.math)()
end

local function double_string_function(str)
    return load(("return function(u,v) return %s end"):format(str), "expression", "t", lua_tikz3dtools_parametric.math)()
end

local function triple_string_function(str)
    return load(("return function(u,v,w) return %s end"):format(str), "expression", "t", lua_tikz3dtools_parametric.math)()
end


local function append_point(hash)
    local x              = hash.x
    local y              = hash.y
    local z              = hash.z
    local filloptions   = hash.filloptions
    local transformation = hash.transformation

    x = single_string_expression(x)
    y = single_string_expression(y)
    z = single_string_expression(z)
    transformation = single_string_expression(transformation)

    local A = { { x, y, z, 1 } }
    local the_segment = mm.matrix_multiply(A, transformation)
    table.insert(
        segments, 
        { 
            segment = the_segment, 
            filloptions = filloptions, 
            type = "point" 
        }
    )
end
register_tex_cmd(
    "appendpoint",
    function()
        append_point{
            x              = token.get_macro("luatikztdtools@p@p@x"),
            y              = token.get_macro("luatikztdtools@p@p@y"),
            z              = token.get_macro("luatikztdtools@p@p@z"),
            filloptions    = token.get_macro("luatikztdtools@p@p@filloptions"),
            transformation = token.get_macro("luatikztdtools@p@p@transformation")
        } 
    end,
    { }
)









local function append_surface(hash)
    local ustart         = hash.ustart
    local ustop          = hash.ustop
    local usamples       = hash.usamples
    local vstart         = hash.vstart
    local vstop          = hash.vstop
    local vsamples       = hash.vsamples
    local x              = hash.x
    local y              = hash.y
    local z              = hash.z
    local drawoptions    = hash.drawoptions
    local filloptions    = hash.filloptions
    local transformation = hash.transformation
    local filter         = hash.filter

    x = double_string_function(x)
    y = double_string_function(y)
    z = double_string_function(z)

    ustart = single_string_expression(ustart)
    ustop = single_string_expression(ustop)
    usamples = single_string_expression(usamples)
    vstart = single_string_expression(vstart)
    vstop = single_string_expression(vstop)
    vsamples = single_string_expression(vsamples)
    transformation = single_string_expression(transformation)

    local ustep = (ustop - ustart) / (usamples - 1)
    local vstep = (vstop - vstart) / (vsamples - 1)

    local function parametric_surface(u, v)
        return { x(u,v), y(u,v), z(u,v), 1 }
    end

    local function distance(P1, P2)
        return math.sqrt(
            (P1[1]-P2[1])^2 +
            (P1[2]-P2[2])^2 +
            (P1[3]-P2[3])^2
        )
    end
    local eps = 0.0000001

    for i = 0, usamples - 2 do
        local u = ustart + i * ustep
        for j = 0, vsamples - 2 do
            local v = vstart + j * vstep
            local A = parametric_surface(u, v)
            local B = parametric_surface(u + ustep, v)
            local C = parametric_surface(u, v + vstep)
            local D = parametric_surface(u + ustep, v + vstep)
            local the_segment1 = mm.matrix_multiply({ A, B, D },transformation)
            local the_segment2 = mm.matrix_multiply({ C, A, D },transformation)
            if not (
                distance(A, B) < eps or
                distance(B, D) < eps or
                distance(A, D) < eps
            ) then
                table.insert(
                    segments, 
                    { 
                        segment      = the_segment1, 
                        filloptions  = filloptions, 
                        type         = "triangle" 
                        ,filter      = filter
                    }
                )
            end
            if not (
                distance(A, C) < eps or
                distance(C, D) < eps or
                distance(A, D) < eps
            ) then
                table.insert(
                    segments, 
                    { 
                        segment      = the_segment2, 
                        filloptions = filloptions, 
                        type         = "triangle",
                        filter       = filter
                    }
                )
            end
        end
    end
end



register_tex_cmd(
    "appendsurface", 
    function()
        append_surface{
            ustart         = token.get_macro("luatikztdtools@p@s@ustart"),
            ustop          = token.get_macro("luatikztdtools@p@s@ustop"),
            usamples       = token.get_macro("luatikztdtools@p@s@usamples"),
            vstart         = token.get_macro("luatikztdtools@p@s@vstart"),
            vstop          = token.get_macro("luatikztdtools@p@s@vstop"),
            vsamples       = token.get_macro("luatikztdtools@p@s@vsamples"),
            x              = token.get_macro("luatikztdtools@p@s@x"),
            y              = token.get_macro("luatikztdtools@p@s@y"),
            z              = token.get_macro("luatikztdtools@p@s@z"),
            transformation = token.get_macro("luatikztdtools@p@s@transformation"),
            filloptions    = token.get_macro("luatikztdtools@p@s@filloptions"),
            filter         = token.get_macro("luatikztdtools@p@s@filter"),
        } 
    end,
    { }
)

local function append_label(hash)
    local x              = hash.x
    local y              = hash.y
    local z              = hash.z
    local name           = hash.name
    local transformation = hash.transformation

    x = single_string_expression(x)
    y = single_string_expression(y)
    z = single_string_expression(z)

    transformation = single_string_expression(transformation)

    local A = { { x, y, z, 1 } }
    local the_segment = mm.matrix_multiply(A, transformation)
    table.insert(
        segments, 
        { 
            segment = the_segment, 
            type = "point",
            name = name
        }
    )
end
register_tex_cmd(
    "appendlabel",
    function()
        append_label{
            x              = token.get_macro("luatikztdtools@p@l@x"),
            y              = token.get_macro("luatikztdtools@p@l@y"),
            z              = token.get_macro("luatikztdtools@p@l@z"),
            name           = token.get_macro("luatikztdtools@p@l@name"),
            transformation = token.get_macro("luatikztdtools@p@l@transformation")
        } 
    end,
    { }
)


local function append_curve(hash)
    local ustart         = hash.ustart
    local ustop          = hash.ustop
    local usamples       = hash.usamples
    local x              = hash.x
    local y              = hash.y
    local z              = hash.z
    local transformation = hash.transformation
    local drawoptions    = hash.drawoptions 
    local arrowtip       = single_string_expression(hash.arrowtip)
    local arrowtail      = single_string_expression(hash.arrowtail)
    local arrowoptions   = hash.arrowoptions

    ustart = single_string_expression(ustart)
    ustop = single_string_expression(ustop)
    usamples = single_string_expression(usamples)
    x = single_string_function(x)
    y = single_string_function(y)
    z = single_string_function(z)
    transformation = single_string_expression(transformation)
    drawoptions = drawoptions

    local ustep = (ustop - ustart) / (usamples - 1)

    local function parametric_curve(u)
        return { { x(u), y(u), z(u), 1 } }
    end

    for i = 0, usamples - 2 do
        local u = ustart + i * ustep
        local A = parametric_curve(u)
        local B = parametric_curve(u+ustep)
        local thesegment = mm.matrix_multiply({ A[1], B[1] }, transformation)
        table.insert(
            segments, 
            { 
                segment = thesegment, 
                drawoptions = drawoptions, 
                type = "line segment" 
            }
        )            
        if i == 0 and arrowtail == true then
            -- append_surface{
                
            -- }
        elseif i == usamples - 2 and arrowtip == true then
            local P = parametric_curve(ustop)
            local NP = parametric_curve(ustop-ustep)
            NP = mm.matrix_subtract(P, NP)
            NP = mm.normalize(NP)
            local U = mm.cross_product(NP, mm.orthogonal_vector(NP))
            local V = mm.cross_product(NP, U)
            append_surface{
                ustart = 0
                ,ustop = 0.1
                ,usamples = 2
                ,vstart = 0
                ,vstop = 1
                ,vsamples = 4
                ,x = "u*cos(v*tau)"
                ,y = "u*sin(v*tau)"
                ,z = "-u"
                ,filloptions = arrowoptions
                ,transformation = ([[
                    matrix_multiply(
                        {
                            {%f,%f,%f,0}
                            ,{%f,%f,%f,0}
                            ,{%f,%f,%f,0}
                            ,{%f,%f,%f,1}
                        }
                        ,%s
                    )
                ]]):format(
                    U[1][1],U[1][2],U[1][3]
                    ,V[1][1],V[1][2],V[1][3]
                    ,NP[1][1],NP[1][2],NP[1][3]
                    ,P[1][1],P[1][2],P[1][3]
                    ,hash.transformation
                )
            }
        end
    end
end
register_tex_cmd(
    "appendcurve",
    function()
        append_curve{
            ustart         = token.get_macro("luatikztdtools@p@c@ustart"),
            ustop          = token.get_macro("luatikztdtools@p@c@ustop"),
            usamples       = token.get_macro("luatikztdtools@p@c@usamples"),
            x              = token.get_macro("luatikztdtools@p@c@x"),
            y              = token.get_macro("luatikztdtools@p@c@y"),
            z              = token.get_macro("luatikztdtools@p@c@z"),
            transformation = token.get_macro("luatikztdtools@p@c@transformation"),
            drawoptions    = token.get_macro("luatikztdtools@p@c@drawoptions"),
            arrowtip       = token.get_macro("luatikztdtools@p@c@arrowtip"),
            arrowtail      = token.get_macro("luatikztdtools@p@c@arrowtail"),
            arrowoptions   = token.get_macro("luatikztdtools@p@c@arrowtipoptions")
        } 
    end,
    { }
)









local function append_solid(hash)
    -- Evaluate bounds and samples
    local ustart   = single_string_expression(hash.ustart)
    local ustop    = single_string_expression(hash.ustop)
    local usamples = single_string_expression(hash.usamples)
    local vstart   = single_string_expression(hash.vstart)
    local vstop    = single_string_expression(hash.vstop)
    local vsamples = single_string_expression(hash.vsamples)
    local wstart   = single_string_expression(hash.wstart)
    local wstop    = single_string_expression(hash.wstop)
    local wsamples = single_string_expression(hash.wsamples)
    local filloptions    = hash.filloptions
    local transformation = single_string_expression(hash.transformation)
    local filter = hash.filter

    -- Compile x,y,z as triple-string functions
    local x = triple_string_function(hash.x)
    local y = triple_string_function(hash.y)
    local z = triple_string_function(hash.z)

    -- Parametric solid function
    local function parametric_solid(u,v,w)
        return { x(u,v,w), y(u,v,w), z(u,v,w), 1 }
    end

    -- Steps for tessellation
    local ustep = (ustop - ustart)/(usamples-1)
    local vstep = (vstop - vstart)/(vsamples-1)
    local wstep = (wstop - wstart)/(wsamples-1)

    -- Helper to tessellate a single face
    local function tessellate_face(fixed_var, fixed_val, s1_start, s1_step, s1_count, s2_start, s2_step, s2_count)
        for i = 0, s1_count-2 do
            local s1a = s1_start + i*s1_step
            local s1b = s1_start + (i+1)*s1_step
            for j = 0, s2_count-2 do
                local s2a = s2_start + j*s2_step
                local s2b = s2_start + (j+1)*s2_step

                -- Compute 4 corners of the quad
                local A,B,C,D
                if fixed_var == "u" then
                    A = parametric_solid(fixed_val, s1a, s2a)
                    B = parametric_solid(fixed_val, s1b, s2a)
                    C = parametric_solid(fixed_val, s1a, s2b)
                    D = parametric_solid(fixed_val, s1b, s2b)
                elseif fixed_var == "v" then
                    A = parametric_solid(s1a, fixed_val, s2a)
                    B = parametric_solid(s1b, fixed_val, s2a)
                    C = parametric_solid(s1a, fixed_val, s2b)
                    D = parametric_solid(s1b, fixed_val, s2b)
                else -- fixed_var == "w"
                    A = parametric_solid(s1a, s2a, fixed_val)
                    B = parametric_solid(s1b, s2a, fixed_val)
                    C = parametric_solid(s1a, s2b, fixed_val)
                    D = parametric_solid(s1b, s2b, fixed_val)
                end

                -- Insert two triangles for the quad
                table.insert(segments, { segment = mm.matrix_multiply({A,B,D}, transformation), filloptions=filloptions, type="triangle", filter = filter })
                table.insert(segments, { segment = mm.matrix_multiply({A,C,D}, transformation), filloptions=filloptions, type="triangle", filter = filter })
            end
        end
    end

    -- Tessellate all 6 faces
    tessellate_face("w", wstart, ustart, ustep, usamples, vstart, vstep, vsamples) -- front
    tessellate_face("w", wstop,  ustart, ustep, usamples, vstart, vstep, vsamples) -- back
    tessellate_face("u", ustart, vstart, vstep, vsamples, wstart, wstep, wsamples) -- left
    tessellate_face("u", ustop,  vstart, vstep, vsamples, wstart, wstep, wsamples) -- right
    tessellate_face("v", vstart, ustart, ustep, usamples, wstart, wstep, wsamples) -- bottom
    tessellate_face("v", vstop,  ustart, ustep, usamples, wstart, wstep, wsamples) -- top
end


register_tex_cmd(
    "appendsolid", 
    function()
        append_solid{
            ustart         = token.get_macro("luatikztdtools@p@solid@ustart"),
            ustop          = token.get_macro("luatikztdtools@p@solid@ustop"),
            usamples       = token.get_macro("luatikztdtools@p@solid@usamples"),
            vstart         = token.get_macro("luatikztdtools@p@solid@vstart"),
            vstop          = token.get_macro("luatikztdtools@p@solid@vstop"),
            vsamples       = token.get_macro("luatikztdtools@p@solid@vsamples"),
            wstart         = token.get_macro("luatikztdtools@p@solid@wstart"),
            wstop          = token.get_macro("luatikztdtools@p@solid@wstop"),
            wsamples       = token.get_macro("luatikztdtools@p@solid@wsamples"),
            x              = token.get_macro("luatikztdtools@p@solid@x"),
            y              = token.get_macro("luatikztdtools@p@solid@y"),
            z              = token.get_macro("luatikztdtools@p@solid@z"),
            transformation = token.get_macro("luatikztdtools@p@solid@transformation"),
            filloptions    = token.get_macro("luatikztdtools@p@solid@filloptions"),
            filter         = token.get_macro("luatikztdtools@p@solid@filter")
        } 
    end,
    { }
)

--- Create a unique signature string for any fragment
--- Ensures line segments are order-independent
--- @param frag table a fragment (point, line segment, or triangle)
--- @param eps number tolerance
--- @return string the signature
local function fragment_signature(frag, eps)
  eps = eps or 1e-7

  local function key_point(p)
    -- round to tolerance to avoid floating noise
    return string.format("%.7f,%.7f,%.7f,%.7f", p[1], p[2], p[3],p[4])
  end

  if frag.type == "point" then
    return "P:" .. key_point(frag.segment[1])

  elseif frag.type == "line segment" then
    local a = key_point(frag.segment[1])
    local b = key_point(frag.segment[2])
    -- unordered: sort lexicographically
    if a < b then
      return "L:" .. a .. "|" .. b
    else
      return "L:" .. b .. "|" .. a
    end

  elseif frag.type == "triangle" then
    local pts = {}
    for i = 1,3 do pts[i] = key_point(frag.segment[i]) end
    table.sort(pts)
    return "T:" .. table.concat(pts, "|")
  end

  return "?"
end



-----------------------

-- ========================================================================
-- Topological sort with partitioning and cycle handling
-- ========================================================================

local function apply_filters(segments)
        local DEBUG = true
        local function dprintf(tag, msg)
            if DEBUG then texio.write_nl(("lua | [%s] %s"):format(tag, msg)) end
        end

        local function fmt_vec(v)
            if type(v) ~= "table" then return tostring(v) end
            local parts = {}
            for i = 1, #v do parts[#parts+1] = tostring(v[i]) end
            return "{" .. table.concat(parts, ", ") .. "}"
        end

        local function fmt_matrix(m)
            if type(m) ~= "table" then return tostring(m) end
            local rows = {}
            for r = 1, #m do rows[#rows+1] = fmt_vec(m[r]) end
            return "[" .. table.concat(rows, " | ") .. "]"
        end

        --   dprintf("FILTER", ("starting filter pass; total segments=%d"):format(#segments))

        local new_segments = {}
        local kept, culled, errors = 0, 0, 0

        for idx, frag in ipairs(segments) do
            local tag = ("frag #%d"):format(idx)

            if frag.type ~= "triangle" or not frag.filter then
            new_segments[#new_segments+1] = frag
            --   dprintf("KEEP", ("%s type=%s no filter -> keep"):format(tag, tostring(frag.type)))
            else
            --   dprintf("FILTER", ("%s raw filter string: %s"):format(tag, tostring(frag.filter)))

            -- homogenize vertices as 1x4 matrices
            local A = { { table.unpack(frag.segment[1]) } }
            local B = { { table.unpack(frag.segment[2]) } }
            local C = { { table.unpack(frag.segment[3]) } }

            --   dprintf("ENV", ("%s vertices: A=%s B=%s C=%s"):format(tag, fmt_matrix(A), fmt_matrix(B), fmt_matrix(C)))
            --   dprintf("ENV", ("%s segment matrix: %s"):format(tag, fmt_matrix(frag.segment)))

            -- build environment
            local env = { A = A, B = B, C = C }
            for k,v in pairs(lua_tikz3dtools_parametric.math) do
                env[k] = v
            end

            local chunk, load_err = load("return (" .. frag.filter .. ")", "filter-chunk", "t", env)
            if not chunk then
                dprintf("ERROR", ("%s load failed: %s"):format(tag, tostring(load_err)))
                new_segments[#new_segments+1] = frag
                errors = errors + 1
            else
                local ok, result = pcall(chunk)
                -- dprintf("EVAL", ("%s ok=%s result(type=%s)=%s"):format(tag, tostring(ok), type(result), tostring(result)))
                if not ok then
                dprintf("ERROR", ("%s runtime error: %s"):format(tag, tostring(result)))
                new_segments[#new_segments+1] = frag
                errors = errors + 1
                elseif result then
                new_segments[#new_segments+1] = frag
                kept = kept + 1
                else
                culled = culled + 1
                end
            end
            end
        end

        --   dprintf("FILTER", ("filtering done; kept=%d culled=%d errors=%d"):format(kept, culled, errors))
        return new_segments
end


-- local function topo_sort_with_cycles(items, cmp)
--     max_depth = 16
--     local eps_local = 1e-7   -- less strict tolerance
--     local DEBUG = false

--     local global_seen = {}   -- local deduplication per call

--     -- ---------------------------------------------------------------
--     -- type checks
--     -- ---------------------------------------------------------------
--     local function is_line_segment(p) return p.type == "line segment" or p.type == "line" end
--     local function is_triangle(p) return p.type == "triangle" end
--     local function is_point(p) return p.type == "point" end

--     local function get_bbox(p)
--         local minX, maxX = math.huge, -math.huge
--         local minY, maxY = math.huge, -math.huge
--         for _, P in ipairs(p.segment or {}) do
--         local x, y = P[1] or 0, P[2] or 0
--         minX, maxX = math.min(minX, x), math.max(maxX, x)
--         minY, maxY = math.min(minY, y), math.max(maxY, y)
--         end
--         return minX, maxX, minY, maxY
--     end

--     local function bboxes_overlap(b1, b2)
--         local minX1, maxX1, minY1, maxY1 = table.unpack(b1)
--         local minX2, maxX2, minY2, maxY2 = table.unpack(b2)
--         return not (maxX1 < minX2+eps or maxX2 < minX1 or maxY1 < minY2 or maxY2 < minY1)
--     end

--     -- ---------------------------------------------------------------
--     -- cloning / projection helpers
--     -- ---------------------------------------------------------------
--     local next_frag_id = 1
--         local function clone_fragment(orig)
--         local copy = {}
--         for k,v in pairs(orig) do if k ~= "segment" then copy[k] = v end end
--         copy.segment = {}
--         for _, P in ipairs(orig.segment or {}) do
--             copy.segment[#copy.segment+1] = {P[1], P[2], P[3], P[4]}
--         end
--         copy.__id = next_frag_id; next_frag_id = next_frag_id + 1
--         return copy
--         end




--         -- Normalize caller-supplied arrays of points but preserve original w exactly.
--         -- We ensure x,y,z are numbers (default 0 if missing) but do NOT change w.
--         local function normalize_points(arr)
--         local out = {}
--         for _, P in ipairs(arr or {}) do
--             if #P >= 3 then
--             local x = P[1]
--             local y = P[2]
--             local z = P[3]
--             local w = P[4]
--             out[#out+1] = { x, y, z, w }
--             end
--         end
--         return out
--         end


--         local function make_fragment_from_points(pts, orig, parent)
--         local typ = (#pts==1 and "point") or (#pts==2 and "line segment") or "triangle"
--         local frag = {}
--         for k,v in pairs(orig or {}) do frag[k] = v end

--         frag.segment = {}
--         for _, P in ipairs(pts) do
--             -- preserve w exactly (do not default to 1)
--             local x = P[1]
--             local y = P[2]
--             local z = P[3]
--             local w = P[4]
--             frag.segment[#frag.segment+1] = { x, y, z, w }
--         end

--         frag.type = typ
--         frag.__id = next_frag_id; next_frag_id = next_frag_id + 1

--         -- merge parent metadata into frag without a whitelist
--         if parent then
--         for k, v in pairs(parent) do
--             if k ~= "segment" and frag[k] == nil then
--             if type(v) == "table" then
--                 -- shallow copy to avoid reference sharing
--                 local cp = {}
--                 for i = 1, #v do cp[i] = v[i] end
--                 frag[k] = cp
--             else
--                 frag[k] = v
--             end
--             end
--         end
--         end


--         return frag
--         end


--         local function fragments_from_partition_result(result, orig, parent)
--         local frags = {}
--         if not result then return frags end
--         for _, pts in pairs(result) do
--             local eu = normalize_points(pts)  -- will preserve w
--             if #eu > 0 then
--             frags[#frags+1] = make_fragment_from_points(eu, orig, parent)
--             end
--         end
--         return frags
--         end

--         local working = {}
--         for i=1,#items do
--         local c = clone_fragment(items[i])
--         working[#working+1] = c
--         end

--         -- -- Copy originals to act as parents (never reinserted)
--         local parents = {}
--         for i = 1, #working do
--             parents[i] = working[i]
--         end

--         -- Cut a single target using a single parent
--         -- Returns: changed:boolean, kept:table|nil
--         local function cut_once(target, parent)
--         if target == parent then return false, nil end
--         if not bboxes_overlap({get_bbox(target)}, {get_bbox(parent)}) then
--             return false, nil
--         end

--         local partA, partB
--         if is_triangle(target) and is_triangle(parent) then
--             partA = triangle_triangle_partition(target.segment, parent.segment)

--         elseif is_line_segment(target) and is_triangle(parent) then
--             partA = line_segment_triangle_clip(target.segment, parent.segment)

--         elseif is_triangle(target) and is_line_segment(parent) then
--             -- no triangle-by-segment partitioner available, so don’t consume the triangle
--             partB = nil

--         elseif is_line_segment(target) and is_line_segment(parent) then
--             partA = line_segment_line_segment_partition(target.segment, parent.segment)
--         end

--         local frags
--         if partA then
--             frags = fragments_from_partition_result(partA, target, parent)
--         elseif partB then
--             frags = fragments_from_partition_result(partB, target, parent)
--         end
--         if not frags or #frags == 0 then
--             return false, nil
--         end

--         -- ✅ strict rule: if any fragment is of a different type, keep the original
--         for i = 1, #frags do
--             if frags[i].type ~= target.type then
--             return false, nil
--             end
--         end

--         -- deduplication
--         local kept, sig_target, any_new = {}, fragment_signature(target, eps), false
--         for _, f in ipairs(frags) do
--             local s = fragment_signature(f, eps)
--             if s ~= sig_target and not global_seen[s] then
--             global_seen[s] = true
--             kept[#kept+1] = f
--             any_new = true
--             end
--         end

--         if not any_new then
--             return false, nil
--         end
--         return true, kept
--         end



--         local changed = true
--         local iter = 0

--         while changed and iter < 50 do
--         iter = iter + 1
--         changed = false

--         local new_working = {}

--         for _, target in ipairs(working) do
--             local fragments = { target }
--             local cut_history = {}

--             while #fragments > 0 do
--             local frag = table.remove(fragments)
--             local sig = fragment_signature(frag, eps)
--             cut_history[sig] = cut_history[sig] or {}

--             local cut_any = false

--             for _, parent in ipairs(parents) do
--                 if not cut_history[sig][parent] then
--                 local did_change, cut = cut_once(frag, parent)
--                 cut_history[sig][parent] = true

--                 if did_change and cut then
--                     for _, f in ipairs(cut) do
--                     fragments[#fragments + 1] = f
--                     end
--                     cut_any = true
--                     changed = true
--                     break -- process newly created fragments immediately
--                 end
--                 end
--             end

--             if not cut_any then
--                 new_working[#new_working + 1] = frag
--             end
--             end
--         end

--         working = new_working
--     end



--     working = apply_filters(working)

--       for _, segment in ipairs(working) do
--         local pts = {}
        
--         -- Reciprocate each point and skip those that are out of bounds
--         for _, P in ipairs(segment.segment) do
--             local R = mm.reciprocate_by_homogenous({P})[1]
            
--             -- Skip points outside the threshold range (culling logic)
--             if math.abs(R[1]) > 150 or math.abs(R[2]) > 150 or R[3] > 0 then
--                 -- Mark segment to be skipped
--                 segment.skip = true
--                 break
--             end
            
--             -- Add reciprocated point to the list
--             table.insert(pts, R)
--         end
        
--         -- Update the segment with the reciprocated points
--         segment.segment = pts
--     end



--   -- ---------------------------------------------------------------
--   -- replace items with working
--   -- ---------------------------------------------------------------
--   for k=#items,1,-1 do items[k]=nil end
--   for _,v in ipairs(working) do items[#items+1]=v end



--   -- ---------------------------------------------------------------
--   -- occlusion graph + SCC topo sort
--   -- ---------------------------------------------------------------
--   local n=#items
--   local bboxes,graph={},{}
--   for i=1,n do bboxes[i]={get_bbox(items[i])}; graph[i]={} end
--   for i=1,n-1 do
--     for j=i+1,n do
--       if bboxes_overlap(bboxes[i],bboxes[j]) then
--         local r=cmp(items[i],items[j])
--         if r==true then table.insert(graph[i],j)
--         elseif r==false then table.insert(graph[j],i) end
--       end
--     end
--   end

--   local index,stack,indices,lowlink,onstack,sccs=0,{}, {}, {}, {}, {}
--   for i=1,n do indices[i],lowlink[i]=-1,-1 end
--   local function dfs(v)
--     indices[v],lowlink[v]=index,index; index=index+1
--     stack[#stack+1]=v; onstack[v]=true
--     for _,w in ipairs(graph[v]) do
--       if indices[w]==-1 then dfs(w); lowlink[v]=math.min(lowlink[v],lowlink[w])
--       elseif onstack[w] then lowlink[v]=math.min(lowlink[v],indices[w]) end
--     end
--     if lowlink[v]==indices[v] then
--       local scc={}
--       while true do
--         local w=table.remove(stack); onstack[w]=false
--         scc[#scc+1]=w
--         if w==v then break end
--       end
--       sccs[#sccs+1]=scc
--     end
--   end
--   for v=1,n do if indices[v]==-1 then dfs(v) end end

--   local scc_index,scc_graph,indeg={}, {}, {}
--   for i,comp in ipairs(sccs) do
--     for _,v in ipairs(comp) do scc_index[v]=i end
--     scc_graph[i],indeg[i]={},0
--   end
--   for v=1,n do
--     for _,w in ipairs(graph[v]) do
--       local si,sj=scc_index[v],scc_index[w]
--       if si~=sj then table.insert(scc_graph[si],sj); indeg[sj]=indeg[sj]+1 end
--     end
--   end

--   local queue,sorted={},{}
--   for i=1,#sccs do if indeg[i]==0 then queue[#queue+1]=i end end
--   while #queue>0 do
--     local i=table.remove(queue,1)
--     for _,v in ipairs(sccs[i]) do sorted[#sorted+1]=items[v] end
--     for _,j in ipairs(scc_graph[i]) do
--       indeg[j]=indeg[j]-1
--       if indeg[j]==0 then queue[#queue+1]=j end
--     end
--   end

--   return sorted
-- end




local function topo_sort_with_cycles(items, cmp)
    local max_depth = 16
    local eps = 1e-7   -- consistent tolerance
    local DEBUG = false

    local global_seen = {}

    -- type checks
    local function is_line_segment(p) return p.type == "line segment" or p.type == "line" end
    local function is_triangle(p) return p.type == "triangle" end
    local function is_point(p) return p.type == "point" end

    local function nearly_equal(a, b, tol) return math.abs((a or 0) - (b or 0)) <= (tol or eps) end

    -- get_bbox expects Cartesian points (after divide-by-w)
    local function get_bbox(p)
        local minX, maxX = math.huge, -math.huge
        local minY, maxY = math.huge, -math.huge
        for _, P in ipairs(p.segment or {}) do
            local x, y = P[1] or 0, P[2] or 0
            minX, maxX = math.min(minX, x), math.max(maxX, x)
            minY, maxY = math.min(minY, y), math.max(maxY, y)
        end
        return minX, maxX, minY, maxY
    end

    local function bboxes_overlap(b1, b2)
        local minX1, maxX1, minY1, maxY1 = table.unpack(b1)
        local minX2, maxX2, minY2, maxY2 = table.unpack(b2)
        return not (maxX1 < minX2 + eps or maxX2 < minX1 or maxY1 < minY2 or maxY2 < minY1)
    end

    -- clone fragment (preserve original w exactly; default nil->1 not forced here)
    local next_frag_id = 1
    local function clone_fragment(orig)
        local copy = {}
        for k,v in pairs(orig) do if k ~= "segment" then copy[k] = v end end
        copy.segment = {}
        for _, P in ipairs(orig.segment or {}) do
            copy.segment[#copy.segment+1] = { P[1], P[2], P[3], P[4] }
        end
        copy.__id = next_frag_id; next_frag_id = next_frag_id + 1
        return copy
    end



    -- create fragment from points, never replace w
    local function make_fragment_from_points(pts, orig, parent)
        local typ = (#pts==1 and "point") or (#pts==2 and "line segment") or "triangle"
        local frag = {}
        for k,v in pairs(orig or {}) do frag[k] = v end
        frag.segment = {}
        for _, P in ipairs(pts) do
            -- preserve the w in P exactly
            frag.segment[#frag.segment+1] = { P[1], P[2], P[3], P[4] }
        end
        frag.type = typ
        frag.__id = next_frag_id; next_frag_id = next_frag_id + 1
        -- copy non-segment parent attributes if missing
        if parent then
            for k,v in pairs(parent) do
                if k ~= "segment" and frag[k] == nil then
                    if type(v) == "table" then
                        local cp = {}
                        for i = 1, #v do cp[i] = v[i] end
                        frag[k] = cp
                    else
                        frag[k] = v
                    end
                end
            end
        end
        return frag
    end

    -- convert partition results to fragments, keeping the exact w of generated points
    local function fragments_from_partition_result(result, orig, parent)
        local frags = {}
        if not result then return frags end
        for _, pts in pairs(result) do
            -- pts already contain correct w, do not modify
            frags[#frags+1] = make_fragment_from_points(pts, orig, parent)
        end
        return frags
    end


    -- prepare working copy with raw Cartesian coords (x,y,z,w if present)
    local working = {}
    for i = 1, #items do working[#working+1] = clone_fragment(items[i]) end
    local parents = {}
    for i = 1, #working do parents[i] = working[i] end


    local function is_degenerate(f)
        if f.type == "triangle" then
            local A, B, C = f.segment[1], f.segment[2], f.segment[3]
            -- compute area in 2D or 3D Cartesian coords
            local ab = {B[1]-A[1], B[2]-A[2], (B[3] or 0)-(A[3] or 0)}
            local ac = {C[1]-A[1], C[2]-A[2], (C[3] or 0)-(A[3] or 0)}
            local area = 0.5 * math.sqrt(
                (ab[2]*ac[3]-ab[3]*ac[2])^2 +
                (ab[3]*ac[1]-ab[1]*ac[3])^2 +
                (ab[1]*ac[2]-ab[2]*ac[1])^2
            )
            return area < eps
        elseif f.type == "line segment" then
            local A, B = f.segment[1], f.segment[2]
            local dx, dy = B[1]-A[1], B[2]-A[2]
            return (dx*dx + dy*dy) < eps*eps
        else
            return false
        end
    end

    local function cut_once(target, parent)
    if target == parent then return false, nil end
    if not bboxes_overlap({ get_bbox(target) }, { get_bbox(parent) }) then return false, nil end

    local partA, partB
    if is_triangle(target) and is_triangle(parent) then
        partA = triangle_triangle_partition(target.segment, parent.segment)
    elseif is_line_segment(target) and is_triangle(parent) then
        partA = line_segment_triangle_clip(target.segment, parent.segment)
    elseif is_triangle(target) and is_line_segment(parent) then
        partB = nil
    elseif is_line_segment(target) and is_line_segment(parent) then
        partA = line_segment_line_segment_partition(target.segment, parent.segment)
    end

    local frags = partA and fragments_from_partition_result(partA, target, parent)
                 or partB and fragments_from_partition_result(partB, target, parent)
                 or {}

    if #frags == 0 then return false, nil end

    local kept = {}
    local sig_target = fragment_signature(target, eps)
    local any_new = false
    local saw_degenerate = false

    for _, f in ipairs(frags) do
        if f.type == "triangle" and #f.segment < 3 then
            saw_degenerate = true  -- mark degenerate
        else
            local s = fragment_signature(f, eps)
            if s ~= sig_target and not global_seen[s] then
                global_seen[s] = true
                kept[#kept+1] = f
                any_new = true
            end
        end
    end

    -- If degenerate was skipped, allow the original-sized fragment
    if saw_degenerate and not global_seen[sig_target] then
        global_seen[sig_target] = true
        kept[#kept+1] = target
        any_new = true
    end

    if not any_new then return false, nil end
    return true, kept
end

    -- iterative cutting (raw coords)
    local changed = true
    local iter = 0
    while changed and iter < 50 do
        iter = iter + 1
        changed = false
        local new_working = {}
        for _, target in ipairs(working) do
            local fragments = { target }
            local cut_history = {}
            while #fragments > 0 do
                local frag = table.remove(fragments)
                local sig = fragment_signature(frag, eps)
                cut_history[sig] = cut_history[sig] or {}
                local cut_any = false
                for _, parent in ipairs(parents) do
                    if not cut_history[sig][parent] then
                        local did_change, cut = cut_once(frag, parent)
                        cut_history[sig][parent] = true
                        if did_change and cut then
                            for _, f in ipairs(cut) do fragments[#fragments+1] = f end
                            cut_any = true
                            changed = true
                            break
                        end
                    end
                end
                if not cut_any then new_working[#new_working+1] = frag end
            end
        end
        working = new_working
    end

    -- SAFE wrapper: ensure apply_filters sees raw Cartesian coordinates and map results back
    local function safe_apply_filters(segments)
        local DEBUG = false
        local function dprintf(tag, msg)
            if DEBUG then texio.write_nl(("lua | [%s] %s"):format(tag, msg)) end
        end

        local function fmt_vec(v)
            if type(v) ~= "table" then return tostring(v) end
            local parts = {}
            for i = 1, #v do parts[#parts+1] = tostring(v[i]) end
            return "{" .. table.concat(parts, ", ") .. "}"
        end

        local function fmt_matrix(m)
            if type(m) ~= "table" then return tostring(m) end
            local rows = {}
            for r = 1, #m do rows[#rows+1] = fmt_vec(m[r]) end
            return "[" .. table.concat(rows, " | ") .. "]"
        end

        local new_segments = {}
        local kept, culled, errors = 0, 0, 0

        for idx, frag in ipairs(segments) do
            local tag = ("frag #%d"):format(idx)

            -- Skip non-triangles or segments without filters
            if frag.type ~= "triangle" or not frag.filter then
                new_segments[#new_segments+1] = frag
            else
                -- homogenize vertices (w preserved)
                local A = {frag.segment[1]}
                local B = {frag.segment[2]}
                local C = {frag.segment[3]}

                -- Build environment for filter
                local env = { A = A, B = B, C = C }
                for k,v in pairs(lua_tikz3dtools_parametric.math) do
                    env[k] = v
                end

                -- Always evaluate filter in Euclidean coordinates
                local chunk, load_err = load("return (" .. frag.filter .. ")", "filter-chunk", "t", env)
                if not chunk then
                    dprintf("ERROR", ("%s load failed: %s"):format(tag, tostring(load_err)))
                    new_segments[#new_segments+1] = frag
                    errors = errors + 1
                else
                    local ok, result = pcall(chunk)
                    if not ok then
                        dprintf("ERROR", ("%s runtime error: %s"):format(tag, tostring(result)))
                        new_segments[#new_segments+1] = frag
                        errors = errors + 1
                    elseif result then
                        new_segments[#new_segments+1] = frag
                        kept = kept + 1
                    else
                        culled = culled + 1
                    end
                end
            end
        end

        -- Return filtered segments (still Euclidean, w untouched)
        return new_segments
    end


    working = safe_apply_filters(working)

        -- NOW do the projective divide (post-filter): convert every segment point to Cartesian (x/w,y/w,z/w)
    for _, segment in ipairs(working) do
        local cart = {}
        for _, P in ipairs(segment.segment or {}) do
            local R = mm.reciprocate_by_homogenous({P})
            cart[#cart+1] = R[1]
        end
        -- assign back the Cartesian coordinates
        segment.segment = cart
    end

    -- prune skipped segments
    local pruned = {}
    for _, s in ipairs(working) do if not s.skip and s.segment and #s.segment > 0 then pruned[#pruned+1] = s end end
    working = pruned

    -- replace items with working (now Cartesian)
    for k = #items, 1, -1 do items[k] = nil end
    for _, v in ipairs(working) do items[#items+1] = v end

    -- occlusion graph + SCC topo sort (bboxes computed on Cartesian coords)
    local n = #items
    local bboxes, graph = {}, {}
    for i = 1, n do bboxes[i] = { get_bbox(items[i]) }; graph[i] = {} end
    for i = 1, n - 1 do
        for j = i + 1, n do
            if bboxes_overlap(bboxes[i], bboxes[j]) then
                local r = cmp(items[i], items[j])
                if r == true then table.insert(graph[i], j)
                elseif r == false then table.insert(graph[j], i) end
            end
        end
    end

    local index, stack, indices, lowlink, onstack, sccs = 0, {}, {}, {}, {}, {}
    for i = 1, n do indices[i], lowlink[i] = -1, -1 end
    local function dfs(v)
        indices[v], lowlink[v] = index, index; index = index + 1
        stack[#stack+1] = v; onstack[v] = true
        for _, w in ipairs(graph[v]) do
            if indices[w] == -1 then dfs(w); lowlink[v] = math.min(lowlink[v], lowlink[w])
            elseif onstack[w] then lowlink[v] = math.min(lowlink[v], indices[w]) end
        end
        if lowlink[v] == indices[v] then
            local scc = {}
            while true do
                local w = table.remove(stack); onstack[w] = false
                scc[#scc+1] = w
                if w == v then break end
            end
            sccs[#sccs+1] = scc
        end
    end
    for v = 1, n do if indices[v] == -1 then dfs(v) end end

    local scc_index, scc_graph, indeg = {}, {}, {}
    for i, comp in ipairs(sccs) do
        for _, v in ipairs(comp) do scc_index[v] = i end
        scc_graph[i], indeg[i] = {}, 0
    end
    for v = 1, n do
        for _, w in ipairs(graph[v]) do
            local si, sj = scc_index[v], scc_index[w]
            if si ~= sj then table.insert(scc_graph[si], sj); indeg[sj] = indeg[sj] + 1 end
        end
    end

    local queue, sorted = {}, {}
    for i = 1, #sccs do if indeg[i] == 0 then queue[#queue+1] = i end end
    while #queue > 0 do
        local i = table.remove(queue, 1)
        for _, v in ipairs(sccs[i]) do sorted[#sorted+1] = items[v] end
        for _, j in ipairs(scc_graph[i]) do indeg[j] = indeg[j] - 1 if indeg[j] == 0 then queue[#queue+1] = j end end
    end

    return sorted
end






local function reverse_inplace(t)
    local n = #t
    for i = 1, math.floor(n/2) do
        t[i], t[n-i+1] = t[n-i+1], t[i]
    end
end









local function sign(number)
    if number > -0.0000001 then return "positive" end
    return "negative"
end


--- Occlusion sorts two points.
--- @param P1 table<table<number>> a point
--- @param P2 table<table<number>> 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.0000001 then 
        if distance({{P1[1][1], P1[1][2], 0, 1}}, {{P2[1][1], P2[1][2], 0, 1}}) < 0.0000001 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 = 0.0000001

    -- 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<table<number>> the point
--- @param L table<table<number>> 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.0000001 then
        local tmp
        if (
            sign(projection[1][1]) == sign(line_direction_vector[1][1]) and
            sign(projection[1][2]) == sign(line_direction_vector[1][2]) and
            sign(projection[1][3]) == sign(line_direction_vector[1][3])
        ) then 
            tmp = 1
        else 
            tmp = -1
        end
        local test = tmp * length(projection) / length(line_direction_vector)
        if (0+eps<=test and test<=1-eps) then
            return point_point_occlusion_sort(P, true_projection)
        end
    end
    return nil
end



--- point versus triangle occlusion comparator
--- @param P table<table<number>> the point
--- @param T table<table<number>> the triangle
--- @return boolean|nil the result of the occlusion comparison
local function point_triangle_occlusion_sort(P, T)

    --[[
        First we project the point and the triangle onto the
        viewing plane (xy). We do this to determine if they occluse 
        using a cross product test. If they do, then we vertically 
        project the point onto the triangle and compare by that.
    ]]
    local point_projection = { {P[1][1], P[1][2], 0, 1} }
    local T1_projection = { {T[1][1], T[1][2], 0, 1} }
    local T2_projection = { {T[2][1], T[2][2], 0, 1} }
    local T3_projection = { {T[3][1], T[3][2], 0, 1} }

    local dist1 = distance(P, {T[1]})
    if dist1 > 0.0000001 == false then return nil end
    local dist2 = distance(P, {T[2]})
    if dist2 > 0.0000001 == false then return nil end
    local dist3 = distance(P, {T[3]})
    if dist3 > 0.0000001 == false then return nil end

    --[[
        For reasons of space, we will invoke smaller names
        for our points.
    ]]
    local P1 = point_projection
    local T1 = T1_projection
    local T2 = T2_projection
    local T3 = T3_projection

    local T1T2 = vector_subtraction(T2, T1)
    local T1P1 = vector_subtraction(P1, T1)
    local T1T2xT1P1 = cross_product(T1T2, T1P1)
    local sign_T1T2xT1P1 = sign(T1T2xT1P1[1][3])

    local T2T3 = vector_subtraction(T3, T2)
    local T2P1 = vector_subtraction(P1, T2)
    local T2T3xT2P1 = cross_product(T2T3, T2P1)
    local sign_T2T3xT2P1 = sign(T2T3xT2P1[1][3])

    if not (sign_T1T2xT1P1 == sign_T2T3xT2P1) then return nil end

    local T3T1 = vector_subtraction(T1, T3)
    local T3P1 = vector_subtraction(P1, T3)
    local T3T1xT3P1 = cross_product(T3T1, T3P1)
    local sign_T3T1xT3P1 = sign(T3T1xT3P1[1][3])

    
    if sign_T2T3xT2P1 == sign_T3T1xT3P1 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+eps<=t and t<=1-eps and 0+eps<=s and s<=1-eps then 
            local a = vector_addition(o, scalar_multiplication(t, u))
            local b = vector_addition(a, scalar_multiplication(s, v))
            return point_point_occlusion_sort(P, b)
        end
    end
    return nil
end




--- line segment versus line segment occlusion sort
--- @param L1 table<table<number>> the first line segment
--- @param L2 table<table<number>> the second line segment
--- @return boolean|nil the result of the occlusion comparison
local function line_segment_line_segment_occlusion_sort(L1, L2)
    local L1A, L1B = {{L1[1][1], L1[1][2], 0, 1}}, {{L1[2][1], L1[2][2], 0, 1}}
    local L2A, L2B = {{L2[1][1], L2[1][2], 0, 1}}, {{L2[2][1], L2[2][2], 0, 1}}
    local L1_dir = vector_subtraction(L1B, L1A)
    local L2_dir = vector_subtraction(L2B, L2A)
    if length(cross_product(L1_dir,L2_dir)) > 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-eps <= t and t<=1+eps and 0-eps <= s and s <= 1+eps) 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


--- line segment versus triangle test
--- @param L table<table<number>> the line segment
--- @param LT table<table<number>> the triangle
--- @return boolean|nil the result of the occlusion comparison
local function line_segment_triangle_occlusion_sort(L, T)
    local A = point_triangle_occlusion_sort({L[1]}, T)
    if A ~= nil then return A end
    local B = point_triangle_occlusion_sort({L[2]}, T)
    if B ~= nil then return B end
    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)
    if C ~= nil then return C end
    local D = line_segment_line_segment_occlusion_sort(L, T2)
    if D ~= nil then return D end
    local E = line_segment_line_segment_occlusion_sort(L, T3)
    if E ~= nil then return E end
    return nil
end


--- triangle versus triangle occlusion sort
--- @param L1 table<table<number>> the first triangle
--- @param L2 table<table<number>> the second triangle
--- @return boolean|nil the result of the occlusion comparison
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)
    if A ~= nil then return A end
    local B = line_segment_triangle_occlusion_sort(T1BC, T2)
    if B ~= nil then return B end
    local C = line_segment_triangle_occlusion_sort(T1CA, T2)
    if C ~= nil then return C end

    local D = point_triangle_occlusion_sort({T1[1]}, T2)
    if D ~= nil then return D end
    local E = point_triangle_occlusion_sort({T1[2]}, T2)
    if E ~= nil then return E end
    local F = point_triangle_occlusion_sort({T1[3]}, T2)
    if F ~= nil then return F end

    local G = point_triangle_occlusion_sort({T2[1]}, T1)
    if G ~= nil then return not G end
    local H = point_triangle_occlusion_sort({T2[2]}, T1)
    if H ~= nil then return not H end
    local I = point_triangle_occlusion_sort({T2[3]}, T1)
    if I ~= nil then return not I end

    return nil
end


--- Occlusion sorts an arbitrary set of segments.
--- @param S1 table<table<number>> a segment
--- @param S2 table<table<number>> 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


local function display_segments()
    
    filtered_segments = topo_sort_with_cycles(segments, occlusion_sort_segments, 16)
    
    -- Reverse the order after sorting
    reverse_inplace(filtered_segments)
local function prune_far_segments(filtered_segments, max_dist)
    max_dist = max_dist or 150
    local pruned = {}
    
    for _, seg in ipairs(filtered_segments) do
        local keep = true
        for _, P in ipairs(seg.segment or {}) do
            local d = math.sqrt(P[1]^2 + P[2]^2 + P[3]^2)
            if d > max_dist then
                keep = false
                break
            end
        end
        if keep then
            pruned[#pruned+1] = seg
        end
    end

    return pruned
end

-- Usage:
filtered_segments = prune_far_segments(filtered_segments, 150)


    -- Iterate through the sorted segments for rendering
    for _, segment in ipairs(filtered_segments) do
        -- Handle point segments
        if segment.type == "point" and not segment.name then
            local P = segment.segment[1]
            tex.sprint(string.format(
                "\\path[%s] (%f,%f) circle[radius = 0.06];",
                segment.filloptions,
                P[1], P[2]
            ))
        end

        if segment.type == "point" and segment.name then
            local P = segment.segment[1]
            tex.sprint(string.format(
                "\\node at (%f,%f) {%s};",
                P[1], P[2]
                ,segment.name
            ))
        end
        
        -- Handle line segments
        if segment.type == "line segment" then 
            local S, E = segment.segment[1], segment.segment[2]
            tex.sprint(string.format(
                "\\path[%s] (%f,%f) -- (%f,%f);",
                segment.drawoptions, 
                S[1], S[2], 
                E[1], E[2]
            ))
        end
        
        -- Handle triangle segments
        if segment.type == "triangle" then 
            local P, Q, R = segment.segment[1], segment.segment[2], segment.segment[3]
            tex.sprint(string.format(
                "\\path[%s] (%f,%f) -- (%f,%f) -- (%f,%f) -- cycle;",
                segment.filloptions,
                P[1], P[2], 
                Q[1], Q[2], 
                R[1], R[2]
            ))
        end
    end
    
    -- Clear the segments array after processing
    segments = {}
end

register_tex_cmd("displaysegments", function() display_segments() end, { })


local function make_object(name, tbl)
    local obj = {}
    obj[name] = tbl
    return obj
end

local function set_matrix(hash)
    local transformation = single_string_expression(hash.transformation)
    local name = hash.name
    local tbl = make_object(name, transformation)
    for k, v in pairs(tbl) do lua_tikz3dtools_parametric.math[k] = v end
end
register_tex_cmd(
    "setobject", 
    function()
        set_matrix{
            name           = token.get_macro("luatikztdtools@p@m@name"),
            transformation = token.get_macro("luatikztdtools@p@m@transformation"),
        } 
    end,
    { }
)