source: OniSplit/Motoko/Quadify.cs

using System;
using System.Collections.Generic;

namespace Oni.Motoko
{
    internal class Quadify
    {
        private readonly Geometry mesh;
        private readonly List<Face> faces;

        public Quadify(Geometry mesh)
        {
            this.mesh = mesh;

            faces = new List<Face>();

            for (int i = 0; i < mesh.Triangles.Length; i += 3)
            {
                var plane = new Plane(
                    mesh.Points[mesh.Triangles[i + 0]],
                    mesh.Points[mesh.Triangles[i + 1]],
                    mesh.Points[mesh.Triangles[i + 2]]);

                faces.Add(new Face(
                    mesh,
                    new[] { mesh.Triangles[i + 0], mesh.Triangles[i + 1], mesh.Triangles[i + 2] },
                    plane.Normal));
            }
        }

        public static List<int[]> Do(Geometry mesh)
        {
            var quadrangulate = new Quadify(mesh);
            return quadrangulate.Execute();
        }

        public List<int[]> Execute()
        {
            GenerateAdjacency();

            var candidates = new List<QuadCandidate>();
            var newFaces = new int[faces.Count][];
            var quadified = new bool[faces.Count];
            int quadCount = 0;

            for (int i = 0; i < faces.Count; i++)
            {
                Face f1 = faces[i];

                if (quadified[i])
                    continue;

                candidates.Clear();

                foreach (Edge e1 in f1.edges)
                {
                    foreach (Edge e2 in e1.adjacency)
                    {
                        if (quadified[faces.IndexOf(e2.face)])
                            continue;

                        candidates.Add(new QuadCandidate(e1, e2));
                    }
                }

                if (candidates.Count > 0)
                {
                    candidates.Sort(new QuadCandidateComparer());
                    newFaces[i] = candidates[0].CreateQuad();

                    int k = faces.IndexOf(candidates[0].e2.face);

                    quadified[i] = true;
                    quadified[k] = true;

                    quadCount++;
                }
            }

            var quadList = new List<int[]>(faces.Count - quadCount);

            for (int i = 0; i < faces.Count; i++)
            {
                if (newFaces[i] != null)
                    quadList.Add(newFaces[i]);
                else if (!quadified[i])
                    quadList.Add(faces[i].indices);
            }

            return quadList;
        }

        private void GenerateAdjacency()
        {
            var points = mesh.Points;
            var pointUseCount = new int[points.Length];
            var pointUsage = new int[points.Length][];

            foreach (Face face in faces)
            {
                foreach (int i in face.indices)
                    pointUseCount[i]++;
            }

            for (int faceIndex = 0; faceIndex < faces.Count; faceIndex++)
            {
                foreach (int pointIndex in faces[faceIndex].indices)
                {
                    int useCount = pointUseCount[pointIndex];
                    int[] usage = pointUsage[pointIndex];

                    if (usage == null)
                    {
                        usage = new int[useCount];
                        pointUsage[pointIndex] = usage;
                    }

                    usage[usage.Length - useCount] = faceIndex;
                    pointUseCount[pointIndex] = useCount - 1;
                }
            }

            var adjacencyBuffer = new List<Edge>();

            foreach (Face f1 in faces)
            {
                foreach (Edge e1 in f1.edges)
                {
                    int[] usage0 = pointUsage[e1.Point0Index];
                    int[] usage1 = pointUsage[e1.Point1Index];

                    if (usage0 == null || usage1 == null)
                        continue;

                    adjacencyBuffer.Clear();

                    foreach (int adjFaceIndex in MatchSortedArrays(usage0, usage1))
                    {
                        Face f2 = faces[adjFaceIndex];

                        if (f2 == f1 || (f2.normal - f1.normal).Length() > 0.01f)
                            continue;

                        foreach (Edge e2 in f2.edges)
                        {
                            if (e1.IsShared(e2))
                                adjacencyBuffer.Add(e2);
                        }
                    }

                    e1.adjacency = adjacencyBuffer.ToArray();
                }
            }
        }

        private static IEnumerable<int> MatchSortedArrays(int[] a1, int[] a2)
        {
            int l1 = a1.Length;
            int l2 = a2.Length;
            int i1 = 0;
            int i2 = 0;

            while (i1 < l1 && i2 < l2)
            {
                int v1 = a1[i1];
                int v2 = a2[i2];

                if (v1 < v2)
                {
                    i1++;
                }
                else if (v1 > v2)
                {
                    i2++;
                }
                else
                {
                    i1++;
                    i2++;

                    yield return v1;
                }
            }
        }

        private class Face
        {
            public readonly Geometry mesh;
            public readonly int[] indices;
            public readonly Edge[] edges;
            public readonly Vector3 normal;

            public Face(Geometry mesh, int[] pointIndices, Vector3 normal)
            {
                this.mesh = mesh;
                this.indices = pointIndices;
                this.normal = normal;

                edges = new Edge[indices.Length];

                for (int i = 0; i < edges.Length; i++)
                    edges[i] = new Edge(this, i);
            }
        }

        private class Edge
        {
            private static readonly Edge[] emptyEdges = new Edge[0];

            public readonly Face face;
            public readonly int i0;
            public readonly int i1;
            public Edge[] adjacency;

            public Edge(Face polygon, int index)
            {
                this.face = polygon;
                this.i0 = index;
                this.i1 = (index + 1) % face.edges.Length;
                this.adjacency = emptyEdges;
            }

            public int Point0Index
            {
                get { return face.indices[i0]; }
            }

            public int Point1Index
            {
                get { return face.indices[i1]; }
            }

            public bool IsShared(Edge e)
            {
                return (Point0Index == e.Point1Index && Point1Index == e.Point0Index);
            }
        }

        private class QuadCandidateComparer : IComparer<QuadCandidate>
        {
            public int Compare(QuadCandidate x, QuadCandidate y)
            {
                return x.length.CompareTo(y.length);
            }
        }

        private class QuadCandidate
        {
            public readonly Edge e1;
            public readonly Edge e2;
            public readonly float length;

            public QuadCandidate(Edge e1, Edge e2)
            {
                this.e1 = e1;
                this.e2 = e2;

                Vector3[] points = e1.face.mesh.Points;
                this.length = (points[e1.Point0Index] - points[e1.Point1Index]).LengthSquared();
            }

            public int[] CreateQuad()
            {
                int[] newPoints = new int[4];
                int l = 0;

                newPoints[l] = e1.face.indices[e1.i1];
                l++;

                for (int k = 0; k < 3; k++)
                {
                    if (k != e1.i0 && k != e1.i1)
                    {
                        newPoints[l] = e1.face.indices[k];
                        l++;

                        break;
                    }
                }

                newPoints[l] = e1.face.indices[e1.i0];
                l++;

                for (int k = 0; k < 3; k++)
                {
                    if (k != e2.i0 && k != e2.i1)
                    {
                        newPoints[l] = e2.face.indices[k];
                        l++;

                        break;
                    }
                }

                return newPoints;
            }
        }
    }
}