//	GGVRMeshes.swift
//
//	© 2024 by Jeff Weeks
//	See TermsOfUse.txt

import simd


struct MeshData {
	let positions: [SIMD3<Double>]
	let facets: [UInt32]
}


// MARK: -
// MARK: Unit sphere

//	What polyhedron should makeUnitSphere() start with?
enum SphereBase {

	//	An icosahedron works best for most purposes
	//	because it distributes the curvature more evenly.
	case icosahedron
	
	//	For the 3-Torus Maze app we start with an octahedron
	//	instead, because its intersections with the coordinate planes
	//	are regular 2ⁿ-gons, which mate perfectly with
	//	cylindrical tubes running parallel to the x, y and z axes.
	//
	//		Note:  The octahedron-based mesh also
	//		offers the possibility of drawing only
	//		a half, a quarter or an eighth of the sphere.
	//
	case octahedron
}

func makeUnitSphere(
	baseShape: SphereBase,
	refinementLevel: Int
) -> MeshData {

	//	Starting with an icosahedron or octahedron, we call
	//	subdivideUnitRadiusMesh() several times to subdivide the mesh,
	//	each time replacing each old facet with four new ones:
	//
	//		  /\		//
	//		 /__\		//
	//		/_\/_\		//
	//
	//	So if we start with an icosahedron, the refinementLevel
	//	determines the number of faces as follows:
	//
	//		0	  20 faces
	//		1	  80 faces
	//		2	 320 faces
	//		3	1280 faces
	//		etc.
	//
	//	Aesthetically, setting refinementLevel = 2 is just barely acceptable.
	//	Setting refinementLevel = 3 looks much better.
	//	The GPU performance penalty for the latter is small:
	//	some quick profiling on my iPhone (under poorly controlled conditions)
	//	found that increasing refinementLevel from 2 to 3
	//
	//	  - significantly increased the time in the vertex program,
	//		but it nevertheless remained small ( < 1 ms in all cases),
	//
	//	  - increased the time in the fragment program
	//			from ~ 7 ms to ~ 9 ms.
	//
	//	so even with refinementLevel = 3 the app had no problem
	//	maintaining its 60 Hz frame rate.
	//
	//
	//	Similar comments apply when starting with an octahedron,
	//	with the number of faces being
	//
	//		0	   8 faces
	//		1	  32 faces
	//		2	 128 faces
	//		3	 512 faces
	//		etc.
	//

	var theMesh: MeshData
	switch baseShape {
	
	case .icosahedron:
		theMesh = makeUnitIcosahedron()
		
	case .octahedron:
		theMesh = makeUnitOctahedron()
	}
	
	for _ in 0 ..< refinementLevel {
	
		theMesh = subdivideUnitRadiusMesh(theMesh)
	}
	
	return theMesh
}

func makeUnitIcosahedron(
) -> MeshData {

	//	The icosahedron's (unnormalized) vertices sit at
	//
	//			( 0, ±1, ±φ)
	//			(±1, ±φ,  0)
	//			(±φ,  0, ±1)
	//
	//	where φ is the golden ratio.  The golden ratio is a root
	//	of the irreducible polynomial φ² - φ - 1,
	//	with numerical value φ = (1 + √5)/2 ≈ 1.6180339887…

	let φ = 0.5 * (1.0 + sqrt(5.0))
		//	golden ratio
		//	≈ 1.61803398874989484820
	
	let theNormalizationFactor = 1.0 / sqrt(φ*φ + 1.0)
		//	inverse length of unnormalized vertex
		//	 = 1/√(φ² + 1) ≈ 0.52573111211913360603
	
	let a = 1.0 * theNormalizationFactor
	let b =  φ  * theNormalizationFactor

	//	The icosahedron's 12 vertices, normalized to unit length
	let thePositions = [
	
		SIMD3<Double>( 0.0, -a,  -b  ),
		SIMD3<Double>( 0.0, +a,  -b  ),
		SIMD3<Double>( 0.0, -a,  +b  ),
		SIMD3<Double>( 0.0, +a,  +b  ),

		SIMD3<Double>( -a,  -b,  0.0 ),
		SIMD3<Double>( +a,  -b,  0.0 ),
		SIMD3<Double>( -a,  +b,  0.0 ),
		SIMD3<Double>( +a,  +b,  0.0 ),

		SIMD3<Double>( -b,  0.0, -a  ),
		SIMD3<Double>( -b,  0.0, +a  ),
		SIMD3<Double>( +b,  0.0, -a  ),
		SIMD3<Double>( +b,  0.0, +a  )
	]

	//	The icosahedron's 20 faces
	let theFacets: [UInt32] = [

		//	Winding order is counterclockwise when viewed
		//	from outside the icosahedron in a right-handed coordinate system.

		//	side-based faces

		 0,  8,  1,
		 1, 10,  0,
		 2, 11,  3,
		 3,  9,  2,

		 4,  0,  5,
		 5,  2,  4,
		 6,  3,  7,
		 7,  1,  6,

		 8,  4,  9,
		 9,  6,  8,
		10,  7, 11,
		11,  5, 10,

		//	corner-based faces
		
		 0,  4,  8,
		 2,  9,  4,
		 1,  8,  6,
		 3,  6,  9,
		 0, 10,  5,
		 2,  5, 11,
		 1,  7, 10,
		 3, 11,  7
	 ]
	
	return MeshData(
			positions: thePositions,
			facets: theFacets)
}

func makeUnitOctahedron(
) -> MeshData {

	//	The octahedron's vertices sit at
	//
	//			(±1,  0,  0)
	//			( 0, ±1,  0)
	//			( 0,  0, ±1)
	//
	let thePositions = [
	
		SIMD3<Double>(-1.0,  0.0,  0.0),
		SIMD3<Double>( 1.0,  0.0,  0.0),
		SIMD3<Double>( 0.0, -1.0,  0.0),
		SIMD3<Double>( 0.0,  1.0,  0.0),
		SIMD3<Double>( 0.0,  0.0, -1.0),
		SIMD3<Double>( 0.0,  0.0,  1.0)
	]

	//	The octahedron's faces
	let theFacets: [UInt32] = [

		//	Winding order is counterclockwise when viewed
		//	from outside the octahedron in a right-handed coordinate system.
		
		//	In RealityKit apps it's not worth messing with partial spheres.
		//	The "Design note" in 3-Torus Maze's makeNodeMesh() function
		//	explains why.

		//		Legacy version:
		//		For the 1/8 sphere at x ≥ 0, y ≥ 0, z ≥ 0,
		//		use only this face (and its subdivisions):
		5, 1, 3,

		//		Legacy version:
		//		For the 1/4 sphere at x ≥ 0, y ≥ 0,
		//		use the previous face along with this one
		//		(and their subdivisions):
		4, 3, 1,

		//		Legacy version:
		//		For the 1/2 sphere at x ≥ 0,
		//		use the previous faces along with these ones
		//		(and their subdivisions):
		4, 1, 2,
		5, 2, 1,
		
		//		Legacy version:
		//		For the full sphere,
		//		use the previous faces along with these ones
		//		(and their subdivisions):
		4, 0, 3,
		4, 2, 0,
		5, 3, 0,
		5, 0, 2
	 ]
	
	return MeshData(
			positions: thePositions,
			facets: theFacets)
}

func subdivideUnitRadiusMesh(
	_ srcMesh: MeshData
) -> MeshData {

	let srcP = srcMesh.positions
	let srcF = srcMesh.facets
	
	let srcVertexCount = srcP.count
	let srcIndexCount = srcF.count
	precondition(srcIndexCount % 3 == 0, "srcMesh.facets.count must be a multiple of 3")
	let srcFacetCount = srcIndexCount / 3

	//	We'll subdivide the mesh, replacing each old facet with four new ones.
	//
	//		  /\		//
	//		 /__\		//
	//		/_\/_\		//
	//
	let dstIndexCount = 4 * srcIndexCount

	//	Euler's formula says that for any cell decomposition of a sphere
	//
	//		v - e + f = 2
	//
	//	Our meshes are comprised exclusively of triangles, so
	//
	//		e = 3f/2
	//
	//	Combining those two formulas gives
	//
	//		v = f/2 + 2
	//
	//	The diagram above shows that
	//
	//		subdivisionF = 4 * originalF
	//
	//	so
	//
	//		subdivisionV = subdivisionF/2 + 2
	//					 = (4*originalF)/2 + 2
	//					 = 2*originalF + 2
	//
	let dstVertexCount =  2 * srcFacetCount  +  2
	
	//	Pre-allocate the output arrays.
	var dstP: [SIMD3<Double>] = Array(repeating: .zero, count: dstVertexCount)
	var dstF: [UInt32] = Array(repeating: 0, count: dstIndexCount)

	//	Keep a running count of the number of vertices we've written into dstP,
	//	so we always know where to write the next one.
	var theVertexCount = 0

	//	First copy the source mesh vertex positions to the destination mesh...
	for i in 0 ..< srcVertexCount {
		dstP[i] = srcP[i]
	}
	theVertexCount += srcVertexCount

	//	...and then create one new vertex on each edge.
	//
	//	Use theNewVertexTable[][] to index them,
	//	so two facets sharing an edge can share the same vertex.
	//
	//	theNewVertexTable[v0][v1] takes the indices v0 and v1 of two old vertices,
	//	and gives the index of the new vertex that sits at the midpoint of the edge
	//	that connects v0 and v1.
	//
	//	The size of theNewVertexTable grows as the square of the number of source vertices.
	//	For modest meshes this won't be a problem.  For larger meshes a fancier algorithm,
	//	with linear rather than quadratic memory demands, could be used.
	//
	//	Initialize all of theNewVertexTable[][]'s entries to nil
	//	to indicate that no new vertices have yet been created.
	//
	var theNewVertexTable: [ [ UInt32? ] ]
		= Array(
			repeating: Array(
					repeating: nil,
					count: srcVertexCount),
			count: srcVertexCount
		)

	//	For each edge in the source mesh, create a new vertex at its midpoint.
	for i in 0 ..< srcFacetCount {
		for j in 0..<3 {
		
			let v0 = Int(srcF[3*i +    j   ])	//	Cast from UInt32 to Int
			let v1 = Int(srcF[3*i + (j+1)%3])
			
			if theNewVertexTable[v0][v1] == nil {
			
				precondition(
					theVertexCount < dstVertexCount,
					"Internal error #1 in subdivideUnitRadiusMesh()")

				//	The new vertex will sit midway between vertices v0 and v1.
				//
				//		Note:  normalize() should be safe here,
				//		because we'll never pass it a zero vector.
				//
				dstP[theVertexCount] = normalize(0.5 * (dstP[v0] + dstP[v1]))

				//	Record the new vertex at [v1][v0] as well as [v0][v1],
				//	so only one new vertex will get created for each edge.
				theNewVertexTable[v0][v1] = UInt32(theVertexCount)
				theNewVertexTable[v1][v0] = UInt32(theVertexCount)

				theVertexCount += 1
			}
		}
	}
	precondition(
		theVertexCount == dstVertexCount,
		"Internal error #2 in subdivideUnitRadiusMesh()")

	//	For each facet in the source mesh,
	//	create four smaller facets in the subdivision.
	for i in 0 ..< srcFacetCount {
	
		//	Call the old vertices incident to this facet v0, v1 and v2
		let v0 = srcF[3*i + 0]
		let v1 = srcF[3*i + 1]
		let v2 = srcF[3*i + 2]

		//	Call the new vertices -- which sit at the midpoints
		//	of the old edges -- vv0, vv1, and vv2.
		//	Each vv_j sits opposite the corresponding v_j .
		guard let vv0 = theNewVertexTable[Int(v1)][Int(v2)],
			  let vv1 = theNewVertexTable[Int(v2)][Int(v0)],
			  let vv2 = theNewVertexTable[Int(v0)][Int(v1)]
		else {
			preconditionFailure("Unexpected nil values in theNewVertexTable[][] in subdivideUnitRadiusMesh()")
		}

		//	Create the new facets.

		dstF[3*(4*i + 0) + 0] = vv0
		dstF[3*(4*i + 0) + 1] = vv1
		dstF[3*(4*i + 0) + 2] = vv2

		dstF[3*(4*i + 1) + 0] = v0
		dstF[3*(4*i + 1) + 1] = vv2
		dstF[3*(4*i + 1) + 2] = vv1

		dstF[3*(4*i + 2) + 0] = v1
		dstF[3*(4*i + 2) + 1] = vv0
		dstF[3*(4*i + 2) + 2] = vv2

		dstF[3*(4*i + 3) + 0] = v2
		dstF[3*(4*i + 3) + 1] = vv1
		dstF[3*(4*i + 3) + 2] = vv0
	}

	return MeshData(
			positions: dstP,
			facets: dstF)
}

// MARK: -
// MARK: Unit cylinder

enum CylinderStyle {

	//	An antiprism works best for most purposes
	//	because it provides twice as many potential
	//	horizon lines as the user walks around it,
	//	compared to a cylinder with the same number
	//	of facets.
	case antiprism
	
	//	For the 3-Torus Maze app we use an ordinary prism
	//	instead, so that it may mate perfectly
	//	with the octahedron-based spheres that it connects to.
	case prism
}

func makeUnitCylinder(
	style cylinderStyle: CylinderStyle,
	numVerticesPerEnd: Int
) -> MeshData {

	let theVertexCount = 2 * numVerticesPerEnd
	let theFacetCount = theVertexCount
	let theIndexCount = 3 * theFacetCount

	let Δθ = 2.0 * Double.pi /  Double(theVertexCount)

	let thePositions = (0 ..< theVertexCount).map() { i in
	
		let theNumHalfSteps: Int
		switch cylinderStyle {
		case .antiprism:	theNumHalfSteps = i
		case .prism:		theNumHalfSteps = (i - i%2)	//	(0,1,2,3,4,…) -> (0,0,2,2,4,…)
		}
		let θ = Double(theNumHalfSteps) * Δθ
		
		return SIMD3<Double>(
			cos(θ),
			sin(θ),
			i%2 == 0 ? +1.0 : -1.0)
	}
	
	var theFacets = Array<UInt32>(repeating: 0, count: theIndexCount)
	for i in 0 ..< theFacetCount {
	
		if i%2 == 0 {
		
			theFacets[3*i + 0] = UInt32(  i  )
			theFacets[3*i + 1] = UInt32(i + 1)
			theFacets[3*i + 2] = UInt32(i + 2)
		
		} else {
		
			theFacets[3*i + 0] = UInt32(i + 1)
			theFacets[3*i + 1] = UInt32(  i  )
			theFacets[3*i + 2] = UInt32(i + 2)
		}
	}
	for i in theIndexCount - 6 ... theIndexCount - 1 {	//	the last few indices "wrap around the cylinder"
	
		if  theFacets[i] >= theVertexCount {
			theFacets[i] -= UInt32(theVertexCount)
		}
	}
	
	return MeshData(
			positions: thePositions,
			facets: theFacets)
}
