Various mesh generation techniques for spheres.

SphereMesh.cs

/*
 * MIT License
 *
 * Copyright (c) 2020 Eric Freed
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 *
 * SphereMesh.cs is a port of https://github.com/caosdoar/spheres.
 */
using System;
using System.Collections.Generic;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;

namespace MyNamespace
{
	/// <summary>
	/// Base class for a standard mesh, containing a set of vertices and elements.
	/// </summary>
	public abstract class Mesh<TVertex, TElement> where TVertex : struct where TElement : struct
	{
		/// <summary>
		/// Gets the number of bytes of the vertex data.
		/// </summary>
		public int VerticesSize => Marshal.SizeOf<TVertex>() * Vertices.Count;

		/// <summary>
		/// Gets the number of bytes of the element data/
		/// </summary>
		public int ElementSize => Marshal.SizeOf<TElement>() * Elements.Count;
		
		/// <summary>
		/// A collection of vertices describing the mesh.
		/// </summary>
		public List<TVertex> Vertices { get; }
		
		/// <summary>
		/// A collection of elements (indices) for rendering the mesh using <c>GL_TRIANGLES</c>.
		/// </summary>
		public List<TElement> Elements { get; }

		/// <summary>
		/// Gets the number of triangles in the mesh.
		/// </summary>
		public int TriangleCount => Elements.Count / 3;

		/// <summary>
		/// Clears the mesh, removing all of tis vertices and elements.
		/// </summary>
		public void Clear()
		{
			Vertices.Clear();
			Elements.Clear();
		}

		/// <summary>
		/// Base class constructor.
		/// </summary>
		protected Mesh()
		{
			Vertices = new List<TVertex>();
			Elements = new List<TElement>();
		}
	}
	
	/// <summary>
	/// Mesh for creating a sphere.
	/// <para>All meshes are built assuming it will rendered using a <c>GL_TRIANGLES</c> primitive type.</para>
	/// </summary>
    public class SphereMesh : Mesh<Vector3, uint>
	{
		/// <summary>
		/// The value of π as a single-precision floating point number.
		/// </summary>
		private const float PI = (float) Math.PI;
	    
	    /// <summary>
	    /// This method divides the sphere using meridians (lines from pole to pole) and parallels (lines parallel to
	    /// the equator). It produces faces with bigger area near the equator and smaller ones close to the poles.
	    /// <para>
	    /// This is the most common implementation of a sphere mesh, and is sometimes known as a "UV Sphere" in other
	    /// 3D tools such as Blender.
	    /// </para>
	    /// </summary>
	    /// <param name="meridians">The number lines of longitude to subdivide the sphere.</param>
	    /// <param name="parallels">The number lines of latitude to subdivide the sphere.</param>
	    /// <returns>The generated mesh.</returns>
	    public static SphereMesh Standard(uint meridians, uint parallels)
		{
			var mesh = new SphereMesh();
			
			mesh.Vertices.Add(new Vector3(0.0f, 1.0f, 0.0f));
			for (var j = 0; j < parallels - 1; ++j)
			{

				var polar = PI * (j + 1) / parallels;
				var sp = (float) Math.Sin(polar);
				var cp = (float) Math.Cos(polar);
			    
				for (var i = 0; i < meridians; ++i)
				{
					var azimuth = 2.0f * PI * i / meridians;
					var sa = (float) Math.Sin(azimuth);
					var ca = (float) Math.Cos(azimuth);
					mesh.Vertices.Add(new Vector3(sp * ca, cp, sp * sa));
				}
			}
			mesh.Vertices.Add(new Vector3(0.0f, -1.0f, 0.0f));
		    
			for (var i = 0u; i < meridians; ++i)
			{
				var a = i + 1u;
				var b = (i + 1u) % meridians + 1u;
				mesh.AddTriangle(0, b, a);
			}
		    
			for (var j = 0u; j < parallels - 2u; ++j)
			{
				var aStart = j * meridians + 1u;
				var bStart = (j + 1u) * meridians + 1u;
				for (var i = 0u; i < meridians; ++i)
				{
					var a = aStart + i;
					var a1 = aStart + (i + 1u) % meridians;
					var b = bStart + i;
					var b1 = bStart + (i + 1u) % meridians;
					mesh.AddQuad(a, a1, b1, b);
				}
			}
			
			for (var i = 0u; i < meridians; ++i)
			{
				var a = i + meridians * (parallels - 2u) + 1u;
				var b = (i + 1) % meridians + meridians * (parallels - 2u) + 1u;
				mesh.AddTriangle((uint) mesh.Vertices.Count - 1u, a, b);
			}
			
			return mesh;
		}
		
	    /// <summary>
	    /// This method uses a uniformly subdivided cube where each vertex position is normalized and multiplied by the
	    /// sphere radius. This creates a non uniformly subdivided sphere where the triangles closer to the center of a
	    /// cube face are bigger than the ones closer to the edges of the cube.
	    /// </summary>
	    /// <param name="divisions">The number of subdivisions to create within the cube, where the higher values will
	    /// result in a "smoother" cube</param>
	    /// <returns>The generated mesh.</returns>
		public static SphereMesh NormalizedCube(uint divisions)
		{
			var mesh = new SphereMesh();
			
			var step = 1.0f / divisions;
			var step3 = new Vector3(step, step, step);

			for (var face = 0u; face < 6u; ++face)
			{
				var origin = CUBE_ORIGIN[face];
				var right = CUBE_RIGHT[face];
				var up = CUBE_UP[face];
				for (var j = 0u; j < divisions + 1u; ++j)
				{
					var j3 = new Vector3(j, j, j);
					for (var i = 0u; i < divisions + 1u; ++i)
					{
						var i3 = new Vector3(i, i, i);
						var p = origin + step3 * (i3 * right + j3 * up);
						mesh.Vertices.Add(Vector3.Normalize(p));
					}
				}
			}

			var k = divisions + 1u;
			for (var face = 0u; face < 6u; ++face)
			{
				for (var j = 0u; j < divisions; ++j)
				{
					var bottom = j < (divisions / 2u);
					for (var i = 0u; i < divisions; ++i)
					{
						var left = i < (divisions / 2u);
						var a = (face * k + j) * k + i;
						var b = (face * k + j) * k + i + 1u;
						var c = (face * k + j + 1u) * k + i;
						var d = (face * k + j + 1u) * k + i + 1u;
						if (bottom ^ left) 
							mesh.AddQuadAlt(a, c, d, b);
						else 
							mesh.AddQuad(a, c, d, b);
					}
				}
			}
			
			return mesh;
		}
		
	    /// <summary>
	    /// This method is based on a subdivided cube as well but it tries to create more uniform divisions in the
	    /// sphere. The area of the faces and the length of the edges suffer less variation, but the sphere still has
	    /// some obvious deformation as points get closer to the corners of the original cube.
	    /// </summary>
	    /// <param name="divisions">The number of subdivisions to create within the cube, where the higher values will
	    /// result in a "smoother" cube</param>
	    /// <returns>The generated mesh.</returns>
		public static SphereMesh SpherifiedCube(uint divisions)
		{
			var mesh = new SphereMesh();
			var step = 1.0f / divisions;
			var step3 = new Vector3(step, step, step);

			for (var face = 0u; face < 6u; ++face)
			{
				var origin = CUBE_ORIGIN[face];
				var right = CUBE_RIGHT[face];
				var up = CUBE_UP[face];
			    
				for (var j = 0u; j < divisions + 1u; ++j)
				{
					var j3 = new Vector3(j, j, j);
					for (uint i = 0; i < divisions + 1; ++i)
					{
						var i3 = new Vector3(i, i, i);
						var p = origin + step3 * (i3 * right + j3 * up);
						var p2 = p * p;
						var n = new Vector3
						(
							p.X * (float) Math.Sqrt(1.0f - 0.5f * (p2.Y + p2.Z) + p2.Y * p2.Z / 3.0f),
							p.Y * (float) Math.Sqrt(1.0f - 0.5f * (p2.Z + p2.X) + p2.Z * p2.X / 3.0f),
							p.Z * (float) Math.Sqrt(1.0f - 0.5f * (p2.X + p2.Y) + p2.X * p2.Y / 3.0f)
						);
						mesh.Vertices.Add(n);
					}
				}
			}

			var k = divisions + 1u;
			for (var face = 0u; face < 6u; ++face)
			{
				for (var j = 0u; j < divisions; ++j)
				{
					var bottom = j < (divisions / 2u);
					for (var i = 0u; i < divisions; ++i)
					{
						var left = i < (divisions / 2u);
						var a = (face * k + j) * k + i;
						var b = (face * k + j) * k + i + 1u;
						var c = (face * k + j + 1u) * k + i;
						var d = (face * k + j + 1u) * k + i + 1u;
						if (bottom ^ left) 
							mesh.AddQuadAlt(a, c, d, b);
						else 
							mesh.AddQuad(a, c, d, b);
					}
				}
			}
			
			return mesh;
		}

	    /// <summary>
	    /// Initially creates a icosahedron, which is a polyhedron 20 identical equilateral triangles, each with the
	    /// same area and each vertex identical distance from its neighbors.
	    /// <para>
	    /// To get a higher number of triangles we need to subdivide each triangle into four triangles by creating a
	    /// new vertex at the middle point of each edge which is then normalized, to make it lie in the sphere surface.
	    /// Sadly this breaks the initial properties of the icosahedron, the triangles are not equilateral anymore and
	    /// neither the area nor the distance between adjacent vertices is the same across the mesh.
	    /// </para>
	    /// Although it usually has a higher triangle count than other types, this is still typically the best
	    /// approximation and most accurate of a sphere.
	    /// </summary>
	    /// <param name="subdivisions">The number of times to subdivide the mesh, where each pass will generate 4 new
	    /// triangles for each existing one.</param>
	    /// <returns></returns>
		public static SphereMesh Icosahedron(int subdivisions = 0)
		{
			var mesh = new SphereMesh();
			mesh.Vertices.AddRange(ICOSAHEDRON_VERTICES);
			mesh.Elements.AddRange(ICOSAHEDRON_ELEMENTS);
			
			while (subdivisions-- > 0)
				mesh.Subdivide();
			return mesh;
		}

	    /// <summary>
	    /// Scales the mesh multiplying each vertex by the specified <paramref name="scalar"/> value.
	    /// </summary>
	    /// <param name="scalar">The </param>
		public void Scale(float scalar)
		{
			for (var i = 0; i < Vertices.Count; i++)
				Vertices[i] *= scalar;
		}

	    /// <summary>
	    /// Subdivides each triangle in a mesh into 4 new sub-triangles to smooth its edges.
	    /// </summary>
	    [MethodImpl(MethodImplOptions.AggressiveInlining)]
	    private void Subdivide()
		{
			var divisions = new Dictionary<Edge, uint>();
			var elementCount = Elements.Count;
			
			for (var i = 0; i < elementCount; ++i)
			{
				var f0 = Elements[i * 3];
				var f1 = Elements[i * 3 + 1];
				var f2 = Elements[i * 3 + 2];

				var v0 = Vertices[(int) f0];
				var v1 = Vertices[(int) f1];
				var v2 = Vertices[(int) f2];

				var f3 = SubdivideEdge(f0, f1, v0, v1, divisions);
				var f4 = SubdivideEdge(f1, f2, v1, v2, divisions);
				var f5 = SubdivideEdge(f2, f0, v2, v0, divisions);

				AddTriangle(f0, f3, f5);
				AddTriangle(f3, f1, f4);
				AddTriangle(f4, f2, f5);
				AddTriangle(f3, f4, f5);
			}
		}
		
		[MethodImpl(MethodImplOptions.AggressiveInlining)]
		private uint SubdivideEdge(uint f0, uint f1, Vector3 v0, Vector3 v1, IDictionary<Edge, uint> divisions)
		{
			var edge = new Edge(f0, f1);
			if (divisions.TryGetValue(edge, out var value))
				return value;

			var v = Vector3.Normalize(new Vector3(0.5f) * (v0 + v1));
			var f = (uint) Vertices.Count;
			Vertices.Add(v);
			divisions.Add(edge, f);
			return f;
		}

		[MethodImpl(MethodImplOptions.AggressiveInlining)]
        private void AddTriangle(uint a, uint b, uint c)
        {
	        Elements.Add(a);
	        Elements.Add(b);
	        Elements.Add(c);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private void AddQuad(uint a, uint b, uint c, uint d)
        {
	        Elements.Add(a);
	        Elements.Add(b);
	        Elements.Add(c);
	        Elements.Add(a);
	        Elements.Add(c);
	        Elements.Add(d);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private void AddQuadAlt(uint a, uint b, uint c, uint d)
        {
	        Elements.Add(a);
	        Elements.Add(b);
	        Elements.Add(d);
	        Elements.Add(b);
	        Elements.Add(c);
	        Elements.Add(d);
        }
        
        private static readonly Vector3[] CUBE_ORIGIN;
        private static readonly Vector3[] CUBE_RIGHT;
        private static readonly Vector3[] CUBE_UP;
        private static readonly Vector3[] ICOSAHEDRON_VERTICES;
        private static readonly uint[] ICOSAHEDRON_ELEMENTS;
        
	    static SphereMesh()
	    {
		    CUBE_ORIGIN = new[]
		    {
			    new Vector3(-1.0f, -1.0f, -1.0f),
			    new Vector3(1.0f, -1.0f, -1.0f),
			    new Vector3(1.0f, -1.0f, 1.0f),
			    new Vector3(-1.0f, -1.0f, 1.0f),
			    new Vector3(-1.0f, 1.0f, -1.0f),
			    new Vector3(-1.0f, -1.0f, 1.0f)
		    };

		    CUBE_RIGHT = new[]
		    {
			    new Vector3(2.0f, 0.0f, 0.0f),
			    new Vector3(0.0f, 0.0f, 2.0f),
			    new Vector3(-2.0f, 0.0f, 0.0f),
			    new Vector3(0.0f, 0.0f, -2.0f),
			    new Vector3(2.0f, 0.0f, 0.0f),
			    new Vector3(2.0f, 0.0f, 0.0f)
		    };

		    CUBE_UP = new[]
		    {
			    new Vector3(0.0f, 2.0f, 0.0f),
			    new Vector3(0.0f, 2.0f, 0.0f),
			    new Vector3(0.0f, 2.0f, 0.0f),
			    new Vector3(0.0f, 2.0f, 0.0f),
			    new Vector3(0.0f, 0.0f, 2.0f),
			    new Vector3(0.0f, 0.0f, -2.0f)
		    };
		    
		    var t = (1.0f + (float) Math.Sqrt(5.0f)) / 2.0f;
		    
		    ICOSAHEDRON_VERTICES = new []
		    {
			    Vector3.Normalize(new Vector3(-1.0f,  t, 0.0f)),
			    Vector3.Normalize(new Vector3( 1.0f,  t, 0.0f)),
			    Vector3.Normalize(new Vector3(-1.0f, -t, 0.0f)),
			    Vector3.Normalize(new Vector3( 1.0f, -t, 0.0f)),
			    Vector3.Normalize(new Vector3(0.0f, -1.0f,  t)),
			    Vector3.Normalize(new Vector3(0.0f,  1.0f,  t)),
			    Vector3.Normalize(new Vector3(0.0f, -1.0f, -t)),
			    Vector3.Normalize(new Vector3(0.0f,  1.0f, -t)),
			    Vector3.Normalize(new Vector3( t, 0.0f, -1.0f)),
			    Vector3.Normalize(new Vector3( t, 0.0f,  1.0f)),
			    Vector3.Normalize(new Vector3(-t, 0.0f, -1.0f)),
			    Vector3.Normalize(new Vector3(-t, 0.0f,  1.0f))
		    };
		    
		    ICOSAHEDRON_ELEMENTS = new uint[]
		    {
			    0, 11, 5,
			    0, 5, 1,
			    0, 1, 7,
			    0, 7, 10,
			    0, 10, 11,
			    1, 5, 9,
			    5, 11, 4,
			    11, 10, 2,
			    10, 7, 6,
			    7, 1, 8,
			    3, 9, 4,
			    3, 4, 2,
			    3, 2, 6,
			    3, 6, 8,
			    3, 8, 9,
			    4, 9, 5,
			    2, 4, 11,
			    6, 2, 10,
			    8, 6, 7,
			    9, 8, 1
		    };
	    }
	    
	    /// <summary>
	    /// Internal type used for calculating subdivisions.
	    /// </summary>
	    struct Edge : IEquatable<Edge>
	    {
		    private readonly uint v0;
		    private readonly uint v1;

		    public Edge(uint v0, uint v1)
		    {
			    this.v0 = v0 < v1 ? v0 : v1;
			    this.v1 = v0 < v1 ? v1 : v0;
		    }

		    public bool Equals(Edge other) => v0 == other.v0 && v1 == other.v1;

		    public override bool Equals(object obj) => obj is Edge other && Equals(other);

		    public override int GetHashCode()
		    {
			    unchecked
			    {
				    return ((int) v0 * 601) ^ (int) v1;
			    }
		    }
	    }
	    
    }
}