I have a surface composed of triangles. I want to recursively subdivide these triangles into 6 sub-triangles, which are then themselves further subdivided into 6 more and so on.

I want for this subdivision to carry on only until enough detail is captured, and no more after that.

I do not know what the heck I should measure for this. I tried using the laplacian, but it didn't seem to work. Maybe there is a better measurement, or perhaps either my laplaican function could be wrong.

There are two functions called append_surface in the MWE. swap their order to see each's effect. In particular, I've included the current function, along with my experimental function which I made in an attempt to solve this question.

The current function produces:

The experimental function produces:

mwe.tex

\documentclass[ tikz ,border = 1cm ]{standalone} \makeatletter \tikzset{ /lua-tikz3dtools/.is family % abbreviated "td" ,/lua-tikz3dtools/.cd ,/lua-tikz3dtools/parametric/.cd % appreviated "p" ,/lua-tikz3dtools/parametric/matrix/.cd % abbreviated "m" ,object/.code = {\protected@edef\luatikztdtools@p@m@transformation{#1}} ,object = {identity_matrix()} ,name/.code = {\protected@edef\luatikztdtools@p@m@name{#1}} ,name = {Change me!} ,/lua-tikz3dtools/parametric/surface/.cd % abbreviated "s" ,ustart/.estore in = \luatikztdtools@p@s@ustart ,ustop/.estore in = \luatikztdtools@p@s@ustop ,usamples/.estore in = \luatikztdtools@p@s@usamples ,vstart/.estore in = \luatikztdtools@p@s@vstart ,vstop/.estore in = \luatikztdtools@p@s@vstop ,vsamples/.estore in = \luatikztdtools@p@s@vsamples ,x/.code = {\protected@edef\luatikztdtools@p@s@x{#1}} ,y/.code = {\protected@edef\luatikztdtools@p@s@y{#1}} ,z/.code = {\protected@edef\luatikztdtools@p@s@z{#1}} ,transformation/.code = {\protected@edef\luatikztdtools@p@s@transformation{#1}} ,fill options/.code = {\protected@edef\luatikztdtools@p@s@filloptions{#1}} ,ustart = {0} ,ustop = {1} ,usamples = {10} ,vstart = {0} ,vstop = {1} ,vsamples = {10} ,x = {sphere(tau*u, pi*v)[1][1]} ,y = {sphere(tau*u, pi*v)[1][2]} ,z = {sphere(tau*u, pi*v)[1][3]} ,transformation = {identity_matrix()} ,fill options = {fill, fill opacity = 0.3} } \ExplSyntaxOn \lua_load_module:n { mwe-lua } \NewDocumentCommand { \appendsurface } {o} { \group_begin: \tikzset{ /lua-tikz3dtools/parametric/surface/.search~also = {/tikz} ,/lua-tikz3dtools/parametric/surface/.cd ,#1 } \__lua_tikztdtools_appendsurface: \group_end: } \NewDocumentCommand { \displaysegments } {} { \group_begin: \__lua_tikztdtools_displaysegments: \group_end: } \NewDocumentCommand { \setobject } {o} { \group_begin: \tikzset{ /lua-tikz3dtools/parametric/matrix/.search~also = {/tikz} ,/lua-tikz3dtools/parametric/matrix/.cd ,#1 } \__lua_tikztdtools_setobject: \group_end: } \ExplSyntaxOff \makeatother \begin{document} \begin{tikzpicture} \setobject[ name = {T} ,object = { matrix_multiply( euler(pi/2,pi/3,5*pi/6) ,translate(0,0,-15) ) } ] \appendsurface[ ustart = {-7} ,ustop = {7} ,usamples = {7} ,vstart = {-7} ,vstop = {7} ,vsamples = {7} ,transformation = {T} ,x = {u} ,y = {v} ,z = { 3 / ((u)^2+(v)^2+1/2) } ,fill options = { preaction = { fill = red!70!black ,fill opacity = 1 } ,postaction = { draw = black ,ultra thin ,line join = round ,line cap = round } } ] \displaysegments \end{tikzpicture} \end{document} 

mwe-lua.lua

 local eps = 0.0000001 --- 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], 1} } 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], 1} } 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], 1} } 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 --- 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 -- 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 = {} 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>> the 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.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.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.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 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) return load(("return %s"):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 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 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({ A, C, 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" } ) 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" } ) end end end end local function append_surface(hash) local ustart = hash.ustart local ustop = hash.ustop local vstart = hash.vstart local vstop = hash.vstop local x = hash.x local y = hash.y local z = hash.z local filloptions = hash.filloptions local transformation = hash.transformation 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) vstart = single_string_expression(vstart) vstop = single_string_expression(vstop) transformation = single_string_expression(transformation) --- parametric definition of our surface This function maps --- a point from the 2D domain onto the 3D codomain. --- @param u number the first coordinate component --- @param v number the second coordinate component --- @return table<number> the map of the given coordinates onto R3 local function parametric_surface(u, v) return { { x(u,v), y(u,v), z(u,v), 1 } } end --- distance between two DOMAIN points --- @param A table<table<number>> the first point in the DOMAIN --- @param B table<table<number>> the second point in the DOMAIN --- @return number their Euclidean distance in the DOMAIN local function dist2(A,B) local Ax, Ay = A[1], A[2] local Bx, By = B[1], B[2] return math.sqrt((Ax-Bx)^2+(Ay-By)^2) end --- gives the incentre of a triangle on the 2D domain --- @param T table<table<number>> --- @return table<number> the incentre (u,v) local function incentre2(T) local A,B,C = T[1],T[2],T[3] local DAB = dist2(A,B) local DBC = dist2(B,C) local DCA = dist2(C,A) local S = DAB + DBC + DCA if S == 0 then -- degenerate triangle, return centroid as fallback return { (A[1]+B[1]+C[1]) / 3, (A[2]+B[2]+C[2]) / 3 } end -- weighted average: (a*A + b*B + c*C) / (a+b+c) local dotx = (A[1]*DAB + B[1]*DBC + C[1]*DCA) / S local doty = (A[2]*DAB + B[2]*DBC + C[2]*DCA) / S return { dotx, doty } end --- partial derivatives of f(x,y) at the point (x=a,y=b) local function first_partial_derivatives(a,b) local eps = 1e-6 local dfdu = scalar_multiplication( 1/eps ,vector_subtraction( parametric_surface(a+eps, b) ,parametric_surface(a,b) ) ) local dfdv = scalar_multiplication( 1/eps ,vector_subtraction( parametric_surface(a, b+eps) ,parametric_surface(a,b) ) ) return {dfdu, dfdv} end --- second partial derivatives --- @param a number the first coordinate component --- @param b number the second coordinate component --- @return table<nnumber> the second partial derivates local function second_partial_derivatives(a, b) local eps = 1e-6 local FPD1, FPD2 = first_partial_derivatives(a, b)[1], first_partial_derivatives(a, b)[2] local FPD1h = first_partial_derivatives(a+eps, b)[1] local FPD2h = first_partial_derivatives(a, b+eps)[2] local d2fdu2 = scalar_multiplication( 1/eps ,vector_subtraction( FPD1h ,FPD1 ) ) local d2fdv2 = scalar_multiplication( 1/eps ,vector_subtraction( FPD2h ,FPD2 ) ) return {d2fdu2, d2fdv2} end local function laplacian(a,b) local P, Q = second_partial_derivatives(a, b)[1], second_partial_derivatives(a, b)[2] return scalar_multiplication(1/2,vector_addition(P, Q)) end -- local LL = parametric_surface(ustart, vstart) -- local LR = parametric_surface(ustop, vstart) -- local UL = parametric_surface(ustart, vstop) -- local UR = parametric_surface(ustop, vstop) local TRI1 = {{ustart, vstart}, {ustop, vstart}, {ustart, vstop}} local TRI2 = {{ustop, vstop}, {ustop, vstart}, {ustart, vstop}} local tris = {TRI1, TRI2} local max_depth = 4 -- you can tune this local iteration = 0 local min_DT = 0.05 local max_DT = 0.5 local tiny = 1e-12 while iteration < max_depth do local newtris = {} local refined = false for _, tri in ipairs(tris) do local v1, v2, v3 = tri[1], tri[2], tri[3] local P = incentre2(tri) -- Adaptive threshold based on Laplacian local PL = laplacian(P[1], P[2]) local px, py, pz = PL[1][1], PL[1][2], PL[1][3] local Lmag = math.sqrt(px*px + py*py + pz*pz) local invL = 1 / (Lmag + tiny) local DT = math.max(min_DT, math.min(invL, max_DT)) -- Distances to control refinement local D1 = dist2(P, v1) local D2 = dist2(P, v2) local D3 = dist2(P, v3) if (D1 < DT) or (D2 < DT) or (D3 < DT) then table.insert(newtris, tri) else -- Midpoints of edges (not P→vertex midpoints) local e12 = { (v1[1] + v2[1]) * 0.5, (v1[2] + v2[2]) * 0.5 } local e23 = { (v2[1] + v3[1]) * 0.5, (v2[2] + v3[2]) * 0.5 } local e31 = { (v3[1] + v1[1]) * 0.5, (v3[2] + v1[2]) * 0.5 } -- Six subtriangles around P that fill the parent: local t1 = { P, v1, e12 } local t2 = { P, e12, v2 } local t3 = { P, v2, e23 } local t4 = { P, e23, v3 } local t5 = { P, v3, e31 } local t6 = { P, e31, v1 } table.insert(newtris, t1) table.insert(newtris, t2) table.insert(newtris, t3) table.insert(newtris, t4) table.insert(newtris, t5) table.insert(newtris, t6) refined = true end end tris = newtris iteration = iteration + 1 if not refined then break end end for _, tri in ipairs(tris) do -- map UV to 3D local A = parametric_surface(tri[1][1], tri[1][2])[1] local B = parametric_surface(tri[2][1], tri[2][2])[1] local C = parametric_surface(tri[3][1], tri[3][2])[1] -- matrix-multiply by transformation local the_segment = mm.matrix_multiply({A, B, C}, transformation) if not ( distance({A}, {B}) < eps or distance({B}, {C}) < eps or distance({A}, {C}) < eps ) then table.insert( segments, { segment = the_segment, filloptions = filloptions, type = "triangle" } ) 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"), } 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, { } ) local function display_segments() local function comparator(a,b) return mm.midpoint(a.segment)[1][3] < mm.midpoint(b.segment)[1][3] end -- Perform sorting after filtering out skipped segments table.sort(segments, comparator) -- Iterate through the sorted segments for rendering for _, segment in ipairs(segments) do -- 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, { })