Lua: Baricentric trasformation for arbitrary 3d point (position relative to a1b1c1-triangle to a2b2c2-triangle)

local function clamp(value, low1, high1)
	if value<low1 then
		return low1
	end
	if value>high1 then
		return high1
	end
	return value
end

local function remap(value, low1, high1, low2, high2)
	return low2 + (value - low1) * (high2 - low2) / (high1 - low1)
end

local function lerp(low1, high1, value)
	return low1 + (high1 - low1) * value
end

function vec3_length(x1,y1,z1)
	return math.sqrt(x1 ^ 2 + y1 ^ 2 + z1 ^ 2)
end

function vec3_normalize(x1,y1,z1)
	local len = vec3_length(x1,y1,z1)
	x1 = x1 / len
	y1 = y1 / len
	z1 = z1 / len
	return x1,y1,z1
end

function vec3_dot(x1,y1,z1, x2,y2,z2)
	return x1*x2 + y1*y2 + z1*z2
end

function vec3_add(x1,y1,z1, x2,y2,z2)
	return x1+x2, y1+y2, z1+z2
end

function vec3_sub(x1,y1,z1, x2,y2,z2)
	return x1-x2, y1-y2, z1-z2
end

function vec3_scale(x1,y1,z1, f)
	return x1*f, y1*f, z1*f
end

function vec3_cross(x1,y1,z1, x2,y2,z2)
	-- https://gist.github.com/williame/3751097
	return y1*z2-z1*y2, z1*x2-x1*z2, x1*y2-y1*x2
end

function barycentric_transform(p1x,p1y,p1z, a1x,a1y,a1z, b1x,b1y,b1z, c1x,c1y,c1z, a2x,a2y,a2z, b2x,b2y,b2z, c2x,c2y,c2z)
	-- BL: BARICENTRIC calcs:  mathutils.geometry.barycentric_transform
	-- // https://github.com/dfelinto/blender/blob/c4ef90f5a0b1d05b16187eb6e32323defe6461c0/source/blender/blenlib/intern/math_geom.c
	-- a1b1c1 normal and area
	local b1a1x, b1a1y, b1a1z = vec3_sub(b1x,b1y,b1z, a1x,a1y,a1z);
	local c1a1x, c1a1y, c1a1z = vec3_sub(c1x,c1y,c1z, a1x,a1y,a1z);
	local a1b1c1n_x, a1b1c1n_y, a1b1c1n_z = vec3_cross(b1a1x, b1a1y, b1a1z, c1a1x, c1a1y, c1a1z) -- non-normalized Normal
	local a1b1c1n_area = vec3_length(a1b1c1n_x, a1b1c1n_y, a1b1c1n_z) / 2
	if a1b1c1n_area < 0.001 then
		-- normal not calcualtable
		-- error("-- bailing: 1")
		return p1x,p1y,p1z
	end
	local a1b1c1nn_x, a1b1c1nn_y, a1b1c1nn_z = vec3_normalize(a1b1c1n_x, a1b1c1n_y, a1b1c1n_z) -- normalized Normal

	-- a2b2c2 normal and area
	local b2a2x, b2a2y, b2a2z = vec3_sub(b2x,b2y,b2z, a2x,a2y,a2z);
	local c2a2x, c2a2y, c2a2z = vec3_sub(c2x,c2y,c2z, a2x,a2y,a2z);
	local a2b2c2n_x, a2b2c2n_y, a2b2c2n_z = vec3_cross(b2a2x, b2a2y, b2a2z, c2a2x, c2a2y, c2a2z) -- non-normalized Normal
	local a2b2c2n_area = vec3_length(a2b2c2n_x, a2b2c2n_y, a2b2c2n_z) / 2
	if a2b2c2n_area < 0.001 then
		-- normal not calcualtable
		-- error("-- bailing: 2")
		return p1x,p1y,p1z
	end
	local a2b2c2nn_x, a2b2c2nn_y, a2b2c2nn_z = vec3_normalize(a2b2c2n_x, a2b2c2n_y, a2b2c2n_z) -- normalized Normal

	-- projecting point on a1-b1-c1
	-- https://stackoverflow.com/questions/9605556/how-to-project-a-point-onto-a-plane-in-3d
	local p1a1b1c1_dst = vec3_dot(a1b1c1nn_x, a1b1c1nn_y, a1b1c1nn_z, p1x-a1x, p1y-a1y, p1z-a1z)
	local p1pr_x, p1pr_y, p1pr_z = vec3_sub(p1x,p1y,p1z, a1b1c1nn_x*p1a1b1c1_dst,a1b1c1nn_y*p1a1b1c1_dst,a1b1c1nn_z*p1a1b1c1_dst)

	-- u, v, w for projected point in 1st triangle
	-- 3d: https://gamedev.stackexchange.com/questions/23743/whats-the-most-efficient-way-to-find-barycentric-coordinates
	local p1pra1x, p1pra1y, p1pra1z = vec3_sub(p1pr_x,p1pr_y,p1pr_z, a1x,a1y,a1z);

	local d00 = vec3_dot(b1a1x,b1a1y,b1a1z, b1a1x,b1a1y,b1a1z);
	local d01 = vec3_dot(b1a1x,b1a1y,b1a1z, c1a1x,c1a1y,c1a1z);
	local d11 = vec3_dot(c1a1x,c1a1y,c1a1z, c1a1x,c1a1y,c1a1z);
	local d20 = vec3_dot(p1pra1x,p1pra1y,p1pra1z, b1a1x,b1a1y,b1a1z);
	local d21 = vec3_dot(p1pra1x,p1pra1y,p1pra1z, c1a1x,c1a1y,c1a1z);
	local denom = d00 * d11 - d01 * d01
	if math.abs(denom) < 0.001 then
		-- error("-- bailing: 3")
		return p1x,p1y,p1z
	end
	-- vwu
	-- wvu wuv vuw uwv uvw
	local v = (d11 * d20 - d01 * d21) / denom
	local w = (d00 * d21 - d01 * d20) / denom
	local u = 1.0 - w - v

	-- backward calc: P=uA+vB+wC, z-offset (second triangle normal) is same, but scaled by triangle area change
	local p2pr_x = u*a2x + v*b2x + w*c2x
	local p2pr_y = u*a2y + v*b2y + w*c2y
	local p2pr_z = u*a2z + v*b2z + w*c2z
	local zscale = a2b2c2n_area/a1b1c1n_area
	local p2x = p2pr_x + a2b2c2nn_x*p1a1b1c1_dst*zscale
	local p2y = p2pr_y + a2b2c2nn_y*p1a1b1c1_dst*zscale
	local p2z = p2pr_z + a2b2c2nn_z*p1a1b1c1_dst*zscale

	local dbg_tridiff = vec3_length(a1x-a2x,a1y-a2y,a1z-a2z)+vec3_length(b1x-b2x,b1y-b2y,b1z-b2z)+vec3_length(c1x-c2x,c1y-c2y,c1z-c2z)+vec3_length(a2b2c2n_x-a1b1c1n_x,a2b2c2n_y-a1b1c1n_y,a2b2c2n_z-a1b1c1n_z)
	--error("mode:" .. tostring(vMODE) .. " diffs:" .. tostring(dbg_tridiff).."->"..tostring( vec3_length(p2x-p1x, p2y-p1y, p2z-p1z) ).." prdiff:"..tostring( vec3_length(p2pr_x-p1pr_x, p2pr_y-p1pr_y, p2pr_z-p1pr_z) ))
	--error("mode:" .. tostring(vMODE) .. " diffs:" .. tostring(dbg_tridiff).."->".." ww:"..tostring(u+w+v).."/"..tostring(u).."/"..tostring(v).."/"..tostring(w).." zscale:"..tostring(zscale) ))
	return p2x, p2y, p2z
end