
PolyUtil = {}

PolyUtil.vertSprite = nil
PolyUtil.vertSpriteData = nil

function PolyUtil.init()
	PolyUtil.vertSprite = allocSprite()

	PolyUtil.vertSpriteData = {}
	PolyUtil.vertSpriteData.position = vec3(0.0, 0.0, 0.0)
	PolyUtil.vertSpriteData.billboard = 1.0
	PolyUtil.vertSpriteData.billboardPixelPerfect = 1.0
	PolyUtil.vertSpriteData.size = vec2(4.0, 4.0)
	PolyUtil.vertSpriteData.color = Config.uiColors.polyEditorVertex
end

-- Directional lighting is done on the GPU by a single cached microShader. It never touches
-- vertex colors: meshLight() only stores the light uniforms on a mesh, and the microShader
-- multiplies the mesh's colors by them at draw time.
-- Lighting is on only for draws made between PolyUtil.beginLighting() and endLighting();
-- to draw a mesh unlit, simply draw it outside such a region.
PolyUtil.lightShader = nil

-- Compiles the lighting microShader once and caches it. The shader reproduces the old CPU math:
--   factor = (ambient + lightColor * clamp(dot(normal, lightDir), 0, 1)) / 255
--   litColor = vec4(base.rgb * factor, base.a)
-- Lighting lives in microVertex() because normals/colors are vertex attributes; microPixel() is a no-op.
-- It reads worldNormal rather than the raw vertexNormal attribute, so the light follows how the
-- model is actually oriented: rotations are respected, and a model mirrored by a negative scale
-- axis is lit correctly with nothing extra to set.
function PolyUtil.ensureLightShader()
	if PolyUtil.lightShader ~= nil then
		return
	end

	PolyUtil.lightShader = allocMicroShader()

	local err = compileMicroShader(PolyUtil.lightShader, [[
		uniform vec3 lightDir;
		uniform vec3 lightColor;
		uniform vec3 ambientColor;

		void microVertex(){
			float lightDot = clamp(dot(worldNormal, lightDir), 0.0, 1.0);
			vec3 factor = (ambientColor + lightColor * lightDot) / 255.0;
			outVertexColor = vec4(vertexColor.rgb * factor, vertexColor.a);
		}

		void microPixel(){
		}
	]])

	if err ~= "" then
		print("PolyUtil lighting shader compile error:", err)
	end
end

-- Stores a mesh's directional light uniforms. The mesh is actually lit only when drawn
-- between PolyUtil.beginLighting()/endLighting(); outside such a region the uniforms are inert.
--   lightDir      - vec3 light direction in world space
--   lightColor    - vec3/vec4 directional light color, 0..255 per channel
--   ambientColor  - vec3/vec4 scene ambient color, 0..255 per channel
function PolyUtil.meshLight(mesh, lightDir, lightColor, ambientColor)
	PolyUtil.ensureLightShader()

	setMeshMicroShader(mesh, "lightDir", lightDir, MICROSHADER_UNIFORM)
	setMeshMicroShader(mesh, "lightColor", vec3(lightColor.x, lightColor.y, lightColor.z), MICROSHADER_UNIFORM)
	setMeshMicroShader(mesh, "ambientColor", vec3(ambientColor.x, ambientColor.y, ambientColor.z), MICROSHADER_UNIFORM)
end

-- Activates the lighting microShader. Wrap mesh draws with beginLighting()/endLighting():
-- every mesh drawn in between must have gone through meshLight() first, and no
-- sprite/wireframe/line/UI draw should happen inside the region.
function PolyUtil.beginLighting()
	PolyUtil.ensureLightShader()

	setActiveMicroShader(PolyUtil.lightShader)
end

-- Restores the standard mesh shader.
function PolyUtil.endLighting()
	setActiveMicroShader(nil)
end

function PolyUtil.makeQuad()
	local quads = {
		vec3(0.5, -0.5, 0.0), 
		vec3(0.5, 0.5, 0.0), 
		vec3(-0.5, 0.5, 0.0), 
		vec3(-0.5, -0.5, 0.0), 
	}

	return quads
end

function PolyUtil.box(size, subdivs)
	 -- Assumes a global or accessible function 'vec3(x, y, z)' exists
	-- that creates a new vector object with each call.

	local s = size or 1.0
	local half_s = s / 2.0
	-- Ensure at least 1 subdivision (the original cube)
	local n = math.max(1, math.floor(subdivs or 1))

	-- The step size between vertices along a single dimension
	local step = s / n

	-- The list of unique vec3 vertices for the entire subdivided mesh
	local vertices = {}

	-- Face definitions:
	-- { axis_fix: Fixed Axis (1=X, 2=Y, 3=Z), sign: +1/-1, axis_1: U-axis, axis_2: V-axis }
	local faces = {
		-- Face 1: Front (+Z) - Original mapping is fine for CCW
		{ axis_fix = 3, sign = 1, axis_1 = 1, axis_2 = 2 },
		-- Face 2: Back (-Z) - X-traversal must be reversed for CCW
		{ axis_fix = 3, sign = -1, axis_1 = 1, axis_2 = 2 },
		-- Face 3: Right (+X) - Z-traversal must be reversed for CCW
		{ axis_fix = 1, sign = 1, axis_1 = 3, axis_2 = 2 },
		-- Face 4: Left (-X) - Z-traversal must be standard for CCW
		{ axis_fix = 1, sign = -1, axis_1 = 3, axis_2 = 2 },
		-- Face 5: Top (+Y) - Z-traversal must be reversed for CCW
		{ axis_fix = 2, sign = 1, axis_1 = 1, axis_2 = 3 },
		-- Face 6: Bottom (-Y) - Z-traversal must be standard for CCW
		{ axis_fix = 2, sign = -1, axis_1 = 1, axis_2 = 3 },
	}

	-- Helper function to set the vector components based on the face definition
	local function get_vertex_coords(face, p1_idx, p2_idx)
		local coords = {0, 0, 0}
		
		local is_plus_x_face = (face.axis_fix == 1 and face.sign == 1)
		local is_plus_y_face = (face.axis_fix == 2 and face.sign == 1)

		-- The fixed coordinate
		coords[face.axis_fix] = face.sign * half_s
		
		-- Default traversal is -half_s + idx * step

		-- 1. Calculate P1 (U-axis) coordinate
		local p1_val = -half_s + p1_idx * step
		
		-- FIX: For the Back face (-Z), the U-axis (X) must be traversed in reverse for CCW winding.
		local is_minus_z_face = (face.axis_fix == 3 and face.sign == -1)
		if is_minus_z_face then
			 p1_val = half_s - p1_idx * step
		end

		-- 2. Calculate P2 (V-axis) coordinate
		local p2_val = -half_s + p2_idx * step
		
		-- FIX: For +X (Right) and +Y (Top), the V-axis (Z-axis in this case) must be traversed in reverse 
		-- to achieve CCW winding without flipping the quad order.
		if is_plus_x_face or is_plus_y_face then
			-- Note: For +X, axis_2 is Y. For +Y, axis_2 is Z. 
			-- We need to check if the spanning axis is Z (axis index 3).
			if face.axis_2 == 3 then
				p2_val = half_s - p2_idx * step
			-- The Right face (+X) needs its U-axis (Z) reversed, but since axis_1 is Z (index 3) 
			-- for the Right face, we should reverse axis_1 instead of axis_2.
			elseif is_plus_x_face and face.axis_1 == 3 then
				p1_val = half_s - p1_idx * step
			end
		end


		coords[face.axis_1] = p1_val
		coords[face.axis_2] = p2_val
		
		return coords[1], coords[2], coords[3]
	end

	-- Iterate over the 6 faces
	for _, face in ipairs(faces) do
		-- Iterate over the grid segments along the two spanning axes
		for i = 0, n - 1 do -- P1 axis index (0 to n-1)
			for j = 0, n - 1 do -- P2 axis index (0 to n-1)
				-- Get the 4 corners of the current quad segment (i, j)
				-- Standard CCW order: (i, j), (i+1, j), (i+1, j+1), (i, j+1)

				-- V1: (i, j)
				local x1, y1, z1 = get_vertex_coords(face, i, j)
				
				-- V2: (i+1, j)
				local x2, y2, z2 = get_vertex_coords(face, i + 1, j)

				-- V3: (i+1, j+1)
				local x3, y3, z3 = get_vertex_coords(face, i + 1, j + 1)
				
				-- V4: (i, j+1)
				local x4, y4, z4 = get_vertex_coords(face, i, j + 1)

				-- Insert the 4 unique vec3 objects for this quad
				table.insert(vertices, vec3(x1, y1, z1))
				table.insert(vertices, vec3(x2, y2, z2))
				table.insert(vertices, vec3(x3, y3, z3))
				table.insert(vertices, vec3(x4, y4, z4))
			end
		end
	end

	return vertices
end

function PolyUtil.sphere(radius, u_res, v_res)
	-- Assumes a global or accessible function 'vec3(x, y, z)' exists.

	local r = radius or 1.0
	-- Ensure resolutions are at least 3 for a meaningful sphere, and integers.
	local u = math.max(3, math.floor(u_res or 16)) -- Longitude segments (around the equator)
	local v = math.max(3, math.floor(v_res or 8))  -- Latitude segments (from pole to pole)

	local vertices = {}

	-- Angular step sizes
	-- U: Goes from 0 to 2*pi (full circle)
	local u_step = 2 * math.pi / u
	-- V: Goes from 0 to pi (pole to pole)
	local v_step = math.pi / v

	-- Function to calculate coordinates for a given U and V index
	local function get_sphere_coords(i, j)
		local polesOffset = 0.0
		-- if j == 0 then
		--     polesOffset = (QUANTIZE_STEP * 4.0)
		-- elseif j == v then
		--     polesOffset = - (QUANTIZE_STEP * 4.0)
		-- end

		-- i is the U index (longitude), j is the V index (latitude)
		
		-- V Angle (Phi): Latitude. Ranges from 0 (North Pole) to pi (South Pole).
		local phi = (j + polesOffset) * v_step
		
		-- U Angle (Theta): Longitude. Ranges from 0 to 2*pi.
		local theta = i * u_step
		
		-- Spherical to Cartesian conversion for Y-UP pole alignment:
		-- Y (vertical axis) = r * cos(phi)
		-- X and Z are calculated in the horizontal plane (using sin(phi))
		
		local sin_phi = math.sin(phi)
		local cos_phi = math.cos(phi)
		
		-- X-coordinate: 
		local x = r * sin_phi * math.cos(theta)
		
		-- Y-coordinate (Pole Axis): 
		local y = r * cos_phi
		
		-- Z-coordinate:
		local z = r * sin_phi * math.sin(theta)

		return x, y, z
	end

	-- Iterate over the quad segments (faces)
	-- i: U-segments (longitude)
	for i = 0, u - 1 do
		-- j: V-segments (latitude)
		for j = 0, v - 1 do
			-- Define the four corners of the quad (i, j) in CCW order:
			
			-- V1: (i, j)
			local x1, y1, z1 = get_sphere_coords(i, j)
			
			-- V2: (i+1, j) (Next U segment, wrapped)
			local x2, y2, z2 = get_sphere_coords((i + 1) % u, j)

			-- V3: (i+1, j+1)
			local x3, y3, z3 = get_sphere_coords((i + 1) % u, j + 1)
			
			-- V4: (i, j+1)
			local x4, y4, z4 = get_sphere_coords(i, j + 1)
			
			-- Insert the 4 unique vec3 objects for this quad
			table.insert(vertices, vec3(x1, y1, z1))
			table.insert(vertices, vec3(x2, y2, z2))
			table.insert(vertices, vec3(x3, y3, z3))
			table.insert(vertices, vec3(x4, y4, z4))
		end
	end

	return vertices
end

function PolyUtil.cylinder(radius, height, u_res, v_res)
	-- Assumes a global or accessible function 'vec3(x, y, z)' exists.

	local r = radius or 0.5
	local h = height or 1.0
	local half_h = h / 2.0
	
	-- U_res: Radial segments (around the circle/wall)
	local u = math.max(3, math.floor(u_res or 16)) 
	-- V_res: Height segments (along the wall)
	local v = math.max(1, math.floor(v_res or 1))  
	
	local vertices = {}

	-- Angular step size for the circumference (U-axis)
	local u_step = 2 * math.pi / u
	-- Vertical step size for the wall (V-axis)
	local v_step = h / v

	-- === 1. SIDE WALL MESH (Inverted Winding) ===
	local y_start = -half_h 
	
	for j = 0, v - 1 do -- Iterate over V-segments (height)
		local y1 = y_start + j * v_step
		local y2 = y_start + (j + 1) * v_step
		
		for i = 0, u - 1 do -- Iterate over U-segments (radial)
			local angle1 = i * u_step
			local angle2 = (i + 1) * u_step
			
			local x1 = r * math.cos(angle1)
			local z1 = r * math.sin(angle1)
			local x2 = r * math.cos(angle2)
			local z2 = r * math.sin(angle2)

			-- Original CCW order: (x1, y1, z1), (x2, y1, z2), (x2, y2, z2), (x1, y2, z1)
			-- NEW CW order (Inverted):
			
			-- V1: (i, j) - Bottom-start
			table.insert(vertices, vec3(x1, y1, z1)) 
			-- V2: (i, j+1) - Top-start (Swapped V2 and V4)
			table.insert(vertices, vec3(x1, y2, z1)) 
			-- V3: (i+1, j+1) - Top-end (Swapped V3 and V2)
			table.insert(vertices, vec3(x2, y2, z2)) 
			-- V4: (i+1, j) - Bottom-end
			table.insert(vertices, vec3(x2, y1, z2)) 
		end
	end

	-- === 2. TOP CAP MESH (+Y Cap, Inverted Winding) ===
	local y_cap = half_h
	
	for i = 0, u - 1 do -- Iterate over radial segments
		local angle1 = i * u_step
		local angle2 = (i + 1) * u_step
		
		local x1 = r * math.cos(angle1)
		local z1 = r * math.sin(angle1)
		local x2 = r * math.cos(angle2)
		local z2 = r * math.sin(angle2)

		-- Original CCW order (Center, Edge_i, Edge_{i+1}, Center)
		-- Inverted CW order: (Center, Edge_{i+1}, Edge_i, Center)

		-- V1: Center (unique object)
		table.insert(vertices, vec3(0, y_cap, 0)) 
		-- V2: Corner 2 (Edge i+1) - Swapped V2 and V3
		table.insert(vertices, vec3(x2, y_cap, z2))
		-- V3: Corner 1 (Edge i)
		table.insert(vertices, vec3(x1, y_cap, z1))
		-- V4: Center (unique object)
		table.insert(vertices, vec3(0, y_cap, 0))
	end
	
	-- === 3. BOTTOM CAP MESH (-Y Cap, Inverted Winding) ===
	local y_bottom = -half_h
	
	for i = 0, u - 1 do -- Iterate over radial segments
		local angle1 = i * u_step
		local angle2 = (i + 1) * u_step
		
		local x1 = r * math.cos(angle1)
		local z1 = r * math.sin(angle1)
		local x2 = r * math.cos(angle2)
		local z2 = r * math.sin(angle2)

		-- Original CCW order (Center, Edge_{i+1}, Edge_i, Center)
		-- Inverted CW order: (Center, Edge_i, Edge_{i+1}, Center)

		-- V1: Center (unique object)
		table.insert(vertices, vec3(0, y_bottom, 0)) 
		-- V2: Corner 1 (Edge i) - Swapped V2 and V3
		table.insert(vertices, vec3(x1, y_bottom, z1))
		-- V3: Corner 2 (Edge i+1)
		table.insert(vertices, vec3(x2, y_bottom, z2))
		-- V4: Center (unique object)
		table.insert(vertices, vec3(0, y_bottom, 0))
	end

	return vertices
end

function PolyUtil.cone(radius, height, u_res, v_res)
	-- Assumes a global or accessible function 'vec3(x, y, z)' exists.

	local r = radius or 0.5
	local h = height or 1.0
	local half_h = h / 2.0
	
	-- U_res: Radial segments (around the circle/wall)
	local u = math.max(3, math.floor(u_res or 16)) 
	-- V_res: Height segments (along the wall)
	local v = math.max(1, math.floor(v_res or 1))  
	
	local vertices = {}

	-- Angular step size for the circumference (U-axis)
	local u_step = 2 * math.pi / u
	
	-- === 1. SIDE WALL MESH (Tapering from Base to Apex - Inverted Winding) ===
	
	for j = 0, v - 1 do -- Iterate over V-segments (height/taper)
		
		-- t: Taper factor, ranges from 0 (base) to 1 (apex).
		local t1 = j / v
		local t2 = (j + 1) / v 
		
		-- Radius and Y-level for the current segment's bottom (t1)
		local r1 = r * (1 - t1)
		local y1 = -half_h + h * t1
		
		-- Radius and Y-level for the current segment's top (t2)
		local r2 = r * (1 - t2)
		local y2 = -half_h + h * t2
		
		for i = 0, u - 1 do -- Iterate over U-segments (radial)
			local angle1 = i * u_step
			local angle2 = (i + 1) * u_step
			
			-- Base coordinates (r1, y1)
			local x1_r1 = r1 * math.cos(angle1)
			local z1_r1 = r1 * math.sin(angle1)
			local x2_r1 = r1 * math.cos(angle2)
			local z2_r1 = r1 * math.sin(angle2)
			
			-- Top coordinates (r2, y2)
			local x1_r2 = r2 * math.cos(angle1)
			local z1_r2 = r2 * math.sin(angle1)
			local x2_r2 = r2 * math.cos(angle2)
			local z2_r2 = r2 * math.sin(angle2)

			-- Original CCW Order: V1 (i, j), V2 (i+1, j), V3 (i+1, j+1), V4 (i, j+1)
			-- NEW CW Order (Inverted): V1 (i, j), V2 (i, j+1), V3 (i+1, j+1), V4 (i+1, j)

			-- V1: (i, j) - Bottom-start
			table.insert(vertices, vec3(x1_r1, y1, z1_r1)) 
			-- V2: (i, j+1) - Top-start (Swapped V2 and V4)
			table.insert(vertices, vec3(x1_r2, y2, z1_r2)) 
			-- V3: (i+1, j+1) - Top-end (Swapped V3 and V2)
			table.insert(vertices, vec3(x2_r2, y2, z2_r2)) 
			-- V4: (i+1, j) - Bottom-end
			table.insert(vertices, vec3(x2_r1, y1, z2_r1)) 
		end
	end

	-- === 2. BOTTOM CAP MESH (-Y Cap, Inverted Winding) ===
	local y_bottom = -half_h
	
	for i = 0, u - 1 do -- Iterate over radial segments
		local angle1 = i * u_step
		local angle2 = (i + 1) * u_step
		
		local x1 = r * math.cos(angle1)
		local z1 = r * math.sin(angle1)
		local x2 = r * math.cos(angle2)
		local z2 = r * math.sin(angle2)

		-- Original CCW order: (Center, Edge_{i+1}, Edge_i, Center)
		-- Inverted CW order: (Center, Edge_i, Edge_{i+1}, Center)

		-- V1: Center (unique object)
		table.insert(vertices, vec3(0, y_bottom, 0)) 
		-- V2: Corner 1 (Edge i) - Swapped V2 and V3
		table.insert(vertices, vec3(x1, y_bottom, z1))
		-- V3: Corner 2 (Edge i+1)
		table.insert(vertices, vec3(x2, y_bottom, z2))
		-- V4: Center (unique object, for the quad structure)
		table.insert(vertices, vec3(0, y_bottom, 0))
	end

	return vertices
end

function PolyUtil.drawNormals(quads, normals, cameraMat)
	setDepthEnabled(false)
	local div = 1.0 / 4
	for i = 1, #quads / 4 do
		local n = normals[((i - 1) * 4) + 1]
		if vec3Dot(cameraMat.forwardAxis, n) > 0.0 then
			local p = ((quads[((i - 1) * 4) + 1] + quads[((i - 1) * 4) + 2] + quads[((i - 1) * 4) + 3] + quads[((i - 1) * 4) + 4]) * div)

			drawLine(p, p + (n * 0.4))
		end
	end
	setDepthEnabled(true)
end

function PolyUtil.drawVert(vertPos, size, color)
	setDepthEnabled(false)

	PolyUtil.vertSpriteData.size.x = size
	PolyUtil.vertSpriteData.size.y = size

	PolyUtil.vertSpriteData.color = color

	PolyUtil.vertSpriteData.position = vertPos
	setSprite(PolyUtil.vertSprite, PolyUtil.vertSpriteData)
	drawSprite(PolyUtil.vertSprite)

	setDepthEnabled(true)
end

function PolyUtil.drawWireframeQuad(quadId, quads, color)
	setDepthEnabled(false)

	local p1 = quads[((quadId - 1) * 4) + 1]
	local p2 = quads[((quadId - 1) * 4) + 2]
	drawLine(p1, p2, color)

	p1 = quads[((quadId - 1) * 4) + 2]
	p2 = quads[((quadId - 1) * 4) + 3]
	drawLine(p1, p2, color)

	p1 = quads[((quadId - 1) * 4) + 3]
	p2 = quads[((quadId - 1) * 4) + 4]
	drawLine(p1, p2, color)

	p1 = quads[((quadId - 1) * 4) + 4]
	p2 = quads[((quadId - 1) * 4) + 1]
	drawLine(p1, p2, color)

	setDepthEnabled(true)
end

function PolyUtil.getAllOverlappingVerts(mesh, verts)
	outOverlappingVerts = {}
	vertsDict = {}

	for i = 1, #verts do
		local vertId = verts[i]

		local overlappingVerts = getOverlappingVertices(mesh, vertId)
		
		for vIt = 1, #overlappingVerts do
			local vId = overlappingVerts[vIt]
			
			vertsDict[vId] = 1
		end
	end

	for vId, v in pairs(vertsDict) do
		outOverlappingVerts[#outOverlappingVerts + 1] = vId
	end

	return outOverlappingVerts
end

function PolyUtil.selectedVertsToQuads(mesh, selectedVerts)
	-- fill holes?
	-- local vertsDict = {}
	-- for i = 1, #selectedVerts do
	-- 	local vertId = selectedVerts[i]
		
	-- 	local overlapVerts = getOverlappingVertices(mesh, vertId)
	-- 	for j = 1, #overlapVerts do
	-- 		local overlapVertId = overlapVerts[j]
	-- 		vertsDict[overlapVertId] = 1
	-- 	end
	-- end

	-- local verts = {}
	-- for vId, v in pairs(vertsDict) do
	-- 	verts[#verts + 1] = vId
	-- end

	local quadIds = {}
	for i = 1, #selectedVerts do
		local quadId = ceil(selectedVerts[i] / 4)

		local vId1 = ((quadId - 1) * 4) + 1
		local vId2 = ((quadId - 1) * 4) + 2
		local vId3 = ((quadId - 1) * 4) + 3
		local vId4 = ((quadId - 1) * 4) + 4

		local selectedVertsInQuad = 0

		for j = 1, #selectedVerts do
			local vId = selectedVerts[j]

			if vId == vId1 or vId == vId2 or vId == vId3 or vId == vId4 then
				selectedVertsInQuad = selectedVertsInQuad + 1
			end
		end
		
		if selectedVertsInQuad == 4 then
			quadIds[quadId] = 1
		end
	end

	return quadIds
end

function PolyUtil.quadFromPoints(p1, p2, p3, p4)
	local EPSILON = 0.00001

	local original_points = {p1, p2, p3, p4}
	local sortable_list = {}

	-- 1. Calculate the Centroid (C)
	local cx, cy, cz = 0, 0, 0
	for _, p in ipairs(original_points) do
		cx = cx + p.x
		cy = cy + p.y
		cz = cz + p.z
	end
	local C = vec3(cx / 4, cy / 4, cz / 4)

	-- 2. Determine the Normalized Plane Normal (N)
	-- Vectors: V_12 = P2 - P1, V_13 = P3 - P1
	local V_12 = vec3(p2.x - p1.x, p2.y - p1.y, p2.z - p1.z)
	local V_13 = vec3(p3.x - p1.x, p3.y - p1.y, p3.z - p1.z)
	local N = vec3Normalize(vec3Cross(V_12, V_13))

	-- 3. Determine a Normalized Reference Vector (R)
	-- Vector from Centroid (C) to the first point (p1).
	local R = vec3Normalize(vec3(p1.x - C.x, p1.y - C.y, p1.z - C.z))

	-- 4. Calculate the angle for each point and store in the sortable list
	for i, p in ipairs(original_points) do
		-- Vector from Centroid to current point (P_i)
		local V = vec3(p.x - C.x, p.y - C.y, p.z - C.z)
		local V_norm = vec3Normalize(V)
		
		-- Angle magnitude (theta) via Dot Product
		local cos_theta = vec3Dot(R, V_norm)
		-- Clamp to ensure math.acos doesn't error due to float inaccuracies
		cos_theta = min(1.0, max(-1.0, cos_theta))
		local angle_magnitude = acos(cos_theta)

		-- Determine the sign (direction) via Cross Product and Normal Vector
		-- Sign = (R x V_norm) . N
		local sign = vec3Dot(vec3Cross(R, V_norm), N)

		local angle = angle_magnitude
		if sign < 0 then
			angle = -angle_magnitude
		end
		
		-- Store the original point and its angle in a new temporary table
		table.insert(sortable_list, {
			point = p,
			angle = angle
		})
	end

	-- 5. Sort the list based on the calculated angle
	table.sort(sortable_list, function(a, b)
		return a.angle < b.angle
	end)

	-- 6. Extract the original points in the sorted order
	local sorted_points = {}
	for _, item in ipairs(sortable_list) do
		table.insert(sorted_points, item.point)
	end

	return sorted_points
end

function PolyUtil.splitQuads(quadPoints, splitType, cameraMat, t)
	local newQuads = {}
	
	-- - Vertical Split separates L/R, so we sort by U_coord (Width).
	-- - Horizontal Split separates B/T, so we sort by V_coord (Height).
	local primarySortKey = (splitType == "vertical") and "v_coord" or "u_coord"
	local secondarySortKey = (splitType == "vertical") and "u_coord" or "v_coord"
	
	-- Iterate through the points, stepping by 4 for each quad
	for i = 1, #quadPoints, 4 do
		-- 1. Extract the quad vertices
		local P1, P2, P3, P4 = quadPoints[i], quadPoints[i + 1], quadPoints[i + 2], quadPoints[i + 3]

		-- orient quad to respect vertical or horizontal
		local M12 = vec3Lerp(P1, P2, 0.5)
		local M43 = vec3Lerp(P4, P3, 0.5)
		local U_axis = M43 - M12 
		if abs(vec3Dot(cameraMat.upAxis, U_axis)) < 0.2 then
			P1 = quadPoints[i + 3]
			P2 = quadPoints[i]
			P3 = quadPoints[i + 1]
			P4 = quadPoints[i + 2]
		end

		-- 2. Define the Intrinsic Axes (U: Horizontal, V: Vertical)
		local edge_a = P2 - P1
		local edge_b = P4 - P1
		local N = vec3Cross(edge_a, edge_b)

		local M12 = vec3Lerp(P1, P2, 0.5)
		local M43 = vec3Lerp(P4, P3, 0.5)
		local U_axis = M43 - M12 
		local V_axis = vec3Cross(U_axis, N) 
		
		-- 3. Create and populate temporary data structure for sorting
		local vertexData = { {p = P1}, {p = P2}, {p = P3}, {p = P4} }

		-- 4. Calculate projection coordinates
		for _, data in ipairs(vertexData) do
			data.v_coord = vec3Dot(data.p, V_axis) -- Height projection
			data.u_coord = vec3Dot(data.p, U_axis) -- Width projection
		end

		-- 5. Primary Sort: Sort by the axis that defines the two halves (U for Vertical, V for Horizontal)
		table.sort(vertexData, function(a, b)
			return a[primarySortKey] < b[primarySortKey]
		end)

		-- H1/H2 are the points of the first half (e.g., Left side or Bottom side)
		-- H3/H4 are the points of the second half (e.g., Right side or Top side)
		local H1, H2, H3, H4 = vertexData[1], vertexData[2], vertexData[3], vertexData[4]
		
		-- 6. Secondary Re-ordering for Edge Definition
		
		-- Sort the H1/H2 and H3/H4 pairs by the secondary axis to correctly identify the edge points.
		local H1_val, H2_val = H1[secondarySortKey], H2[secondarySortKey]
		local H3_val, H4_val = H3[secondarySortKey], H4[secondarySortKey]
		
		-- E1 and E2 define the FIRST EDGE/SIDE (e.g., Left edge or Bottom side)
		local E1 = H1_val < H2_val and H1.p or H2.p -- E1 is the start of the edge (e.g., BL, or TL)
		local E2 = H1_val < H2_val and H2.p or H1.p -- E2 is the end of the edge (e.g., TR, or BR)

		-- E3 and E4 define the SECOND EDGE/SIDE (e.g., Right edge or Top side)
		local E3 = H3_val < H4_val and H3.p or H4.p
		local E4 = H3_val < H4_val and H4.p or H3.p

		-- 7. Calculate Midpoints and Define New Quads

		if splitType == "vertical" then
			-- Vertical Split (Left/Right)
			-- E1, E2 are the two points on the LEFT SIDE (e.g., BL and TL). E3, E4 are on the RIGHT SIDE (e.g., BR and TR).
			-- The new shared edge connects the Bottom Midpoint (MB) and Top Midpoint (MT).
			
			-- Midpoint MB: between the two bottom points (E1 and E3)
			local MB = vec3Lerp(E1, E3, t) 
			-- Midpoint MT: between the two top points (E2 and E4)
			local MT = vec3Lerp(E2, E4, t) 

			local vA1 = E1
			local vB1 = MB
			local vC1 = MT
			local vD1 = E2

			local vA2 = MB
			local vB2 = E3
			local vC2 = E4
			local vD2 = MT

			local edge_a = E2 - E1
			local edge_b = E4 - E1
			local N = vec3Cross(edge_a, edge_b)

			-- Quad A (Bottom Half): E1, MB, MT, E2 (Winding: Left-Bottom -> Right-Bottom -> Right-Top -> Left-Top)
			-- Quad A (Left Half): E1, E3, MT, E2  (Wrong)
			
			-- Winding for Left Half (using BL, TL, TR, BR convention for the original quad): 
			-- BL(E1), BR(E3), MT(M2), MB(M1)
			
			-- Quad A (Left Half): E1 (BL), MB, MT, E2 (TL)
			table.insert(newQuads, vA1) -- BL
			table.insert(newQuads, vB1) -- New Bottom-Midpoint
			table.insert(newQuads, vC1) -- New Top-Midpoint
			table.insert(newQuads, vD1) -- TL
			
			-- Quad B (Right Half): MB, E3, E4, MT
			table.insert(newQuads, vA2) -- MB
			table.insert(newQuads, vB2) -- BR
			table.insert(newQuads, vC2) -- TR
			table.insert(newQuads, vD2) -- MT

		else -- splitType == "horizontal"
			-- Horizontal Split (Bottom/Top)
			-- E1, E2 are the two points on the BOTTOM SIDE (e.g., BL and BR). E3, E4 are on the TOP SIDE (e.g., TL and TR).
			-- The new shared edge connects the Left Midpoint (ML) and Right Midpoint (MR).
			
			-- Midpoint ML: between the two left points (E1 and E3)
			local ML = vec3Lerp(E1, E3, t)
			-- Midpoint MR: between the two right points (E2 and E4)
			local MR = vec3Lerp(E2, E4, t)

			local vA1 = E1
			local vB1 = E2
			local vC1 = MR
			local vD1 = ML

			local vA2 = ML
			local vB2 = MR
			local vC2 = E4
			local vD2 = E3

			local edge_a = E2 - E1
			local edge_b = E4 - E1
			local N = vec3Cross(edge_a, edge_b)

			-- Quad A (Bottom Half): E1, E2, MR, ML
			table.insert(newQuads, vA1) -- BL
			table.insert(newQuads, vB1) -- BR
			table.insert(newQuads, vC1) -- New Right-Midpoint
			table.insert(newQuads, vD1) -- New Left-Midpoint
			
			-- Quad B (Top Half): ML, MR, E4, E3
			table.insert(newQuads, vA2) -- ML
			table.insert(newQuads, vB2) -- MR
			table.insert(newQuads, vC2) -- TR
			table.insert(newQuads, vD2) -- TL
		end
	end

	return newQuads
end

function PolyUtil.fixNormalsFromNeighbors(selectedVerts, mesh)
	local currentMeshQuads = getMeshVertices(mesh)
	local currentMeshNormals = getMeshVertices(mesh, VERT_NORMAL)

	local quadIds = PolyUtil.selectedVertsToQuads(mesh, selectedVerts)

	for quadId, v in pairs(quadIds) do
		local numFlippedNormals = 0

		for vIt = 1, 4 do
			local vId = ((quadId - 1) * 4) + vIt
			
			local overlappingVerts = getOverlappingVertices(mesh, vId)

			local neighborsNormal = vec3(0.0, 0.0, 0.0)

			for overlapVertIt = 1, #overlappingVerts do
				local overlapVertId = overlappingVerts[overlapVertIt]
				if overlapVertId ~= vId then
					neighborsNormal = neighborsNormal + currentMeshNormals[overlapVertId]
				end
			end

			neighborsNormal = vec3Normalize(neighborsNormal)

			local vertNormal = currentMeshNormals[vId]
			if vec3Dot(vertNormal, neighborsNormal) < 0.0 then
				numFlippedNormals = numFlippedNormals + 1
			end
		end

		if numFlippedNormals > 1 then
			flipNormal(mesh, quadId)
		end
	end
end

function PolyUtil.extrude(selectedVerts, currentMesh)
	local currentMeshQuads = getMeshVertices(currentMesh)
	local currentMeshNormals = getMeshVertices(currentMesh, VERT_NORMAL)

	local quadIds = PolyUtil.selectedVertsToQuads(currentMesh, selectedVerts)

	local extrudeNormal = vec3()

	local numOverlapsPerVert = {}

	for quadId, v in pairs(quadIds) do
		for vIt = 1, 4 do
			local vId = ((quadId - 1) * 4) + vIt

			extrudeNormal = extrudeNormal + currentMeshNormals[vId]

			local overlapVerts = getOverlappingVertices(currentMesh, vId)
			numOverlapsPerVert[vId] = #overlapVerts
		end
	end

	extrudeNormal = vec3Normalize(extrudeNormal)

	local extrudedPointToBasePosMap = {}

	local extrudeOffset = QUANTIZE_STEP * 20.0

	for quadId, v in pairs(quadIds) do
		for vIt = 1, 4 do
			local vId = ((quadId - 1) * 4) + vIt

			local vertPos = currentMeshQuads[vId]
			local newPos = vertPos + (extrudeNormal * extrudeOffset)

			currentMeshQuads[vId] = newPos

			extrudedPointToBasePosMap[vId] = vertPos
		end
	end
	
	-- find boundary edges
	local boundaryEdges = {}

	function isBoundary(vId)
		local overlapVerts = getOverlappingVertices(currentMesh, vId)

		if #overlapVerts <= numOverlapsPerVert[vId] then
			return true
		end

		return false
	end

	function edgeIsShared(quadId, edgeVId1, edgeVId2, quadIds)
	 local hash1 = getVertexHash(currentMesh, edgeVId1)
	 local hash2 = getVertexHash(currentMesh, edgeVId2)

	 for otherQuadId, v in pairs(quadIds) do
		 if otherQuadId ~= quadId then
			 local quadVId1 = ((otherQuadId - 1) * 4) + 1
			 local quadVId2 = ((otherQuadId - 1) * 4) + 2
			 local quadVId3 = ((otherQuadId - 1) * 4) + 3
			 local quadVId4 = ((otherQuadId - 1) * 4) + 4

			 local quadHash1 = getVertexHash(currentMesh, quadVId1)
			 local quadHash2 = getVertexHash(currentMesh, quadVId2)
			 local quadHash3 = getVertexHash(currentMesh, quadVId3)
			 local quadHash4 = getVertexHash(currentMesh, quadVId4)

			 local hash1Shared = false
			 if hash1 == quadHash1 or 
				 hash1 == quadHash2 or
				 hash1 == quadHash3 or
				 hash1 == quadHash4 then

				 hash1Shared = true
			 end

			 local hash2Shared = false
			 if hash2 == quadHash1 or 
				 hash2 == quadHash2 or
				 hash2 == quadHash3 or
				 hash2 == quadHash4 then

				 hash2Shared = true
			 end

			 if hash1Shared == true and hash2Shared == true then
				 return true
			 end
		 end
	 end

	 return false    
	end

	for quadId, v in pairs(quadIds) do
		local vId1 = ((quadId - 1) * 4) + 1
		local vId2 = ((quadId - 1) * 4) + 2
		local vId3 = ((quadId - 1) * 4) + 3
		local vId4 = ((quadId - 1) * 4) + 4

		local v1IsBoundary = isBoundary(vId1);
		local v2IsBoundary = isBoundary(vId2);
		local v3IsBoundary = isBoundary(vId3);
		local v4IsBoundary = isBoundary(vId4);

		if v1IsBoundary == true and v2IsBoundary == true and edgeIsShared(quadId, vId1, vId2, quadIds) == false then
			local edge = {}
			edge.p1 = currentMeshQuads[vId1]
			edge.p2 = currentMeshQuads[vId2]

			edge.vId1 = vId1
			edge.vId2 = vId2

			boundaryEdges[#boundaryEdges + 1] = edge
		end

		if v2IsBoundary == true and v3IsBoundary == true and edgeIsShared(quadId, vId2, vId3, quadIds) == false then
			local edge = {}
			edge.p1 = currentMeshQuads[vId2]
			edge.p2 = currentMeshQuads[vId3]

			edge.vId1 = vId2
			edge.vId2 = vId3

			boundaryEdges[#boundaryEdges + 1] = edge
		end

		if v3IsBoundary == true and v4IsBoundary == true and edgeIsShared(quadId, vId3, vId4, quadIds) == false then
			local edge = {}
			edge.p1 = currentMeshQuads[vId3]
			edge.p2 = currentMeshQuads[vId4]

			edge.vId1 = vId3
			edge.vId2 = vId4

			boundaryEdges[#boundaryEdges + 1] = edge
		end

		if v4IsBoundary == true and v1IsBoundary == true and edgeIsShared(quadId, vId4, vId1, quadIds) == false then
			local edge = {}
			edge.p1 = currentMeshQuads[vId4]
			edge.p2 = currentMeshQuads[vId1]

			edge.vId1 = vId4
			edge.vId2 = vId1

			boundaryEdges[#boundaryEdges + 1] = edge
		end
	end

	local newQuads = {}
	for i = 1, #boundaryEdges do
		local edge = boundaryEdges[i]

		local basePoint1 = extrudedPointToBasePosMap[edge.vId1]
		local basePoint2 = extrudedPointToBasePosMap[edge.vId2]

		local newQuad = PolyUtil.quadFromPoints(edge.p1, edge.p2, basePoint1, basePoint2)

		currentMeshQuads[#currentMeshQuads + 1] = newQuad[4]
		currentMeshQuads[#currentMeshQuads + 1] = newQuad[3]
		currentMeshQuads[#currentMeshQuads + 1] = newQuad[2]
		currentMeshQuads[#currentMeshQuads + 1] = newQuad[1]
	end

	return currentMeshQuads
end
