1 files changed, 323 insertions, 0 deletions
diff --git a/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.cpp b/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.cpp
new file mode 100644
index 0000000000..683ee603cd
--- /dev/null
+++ b/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.cpp
@@ -0,0 +1,323 @@
+// Copyright NVIDIA Corporation 2008 -- Ignacio Castano <icastano@nvidia.com>
+
+#include "nvmesh.h" // pch
+
+#include "ParameterizationQuality.h"
+
+#include "nvmesh/halfedge/Mesh.h"
+#include "nvmesh/halfedge/Face.h"
+#include "nvmesh/halfedge/Vertex.h"
+#include "nvmesh/halfedge/Edge.h"
+
+#include "nvmath/Vector.inl"
+
+#include "nvcore/Debug.h"
+
+#include <float.h>
+
+
+using namespace nv;
+
+#if 0
+/*
+float triangleConformalEnergy(Vector3 q[3], Vector2 p[3])
+{
+const Vector3 v1 = q[0];
+const Vector3 v2 = q[1];
+const Vector3 v3 = q[2];
+
+const Vector2 w1 = p[0];
+const Vector2 w2 = p[1];
+const Vector2 w3 = p[2];
+
+float x1 = v2.x() - v1.x();
+float x2 = v3.x() - v1.x();
+float y1 = v2.y() - v1.y();
+float y2 = v3.y() - v1.y();
+float z1 = v2.z() - v1.z();
+float z2 = v3.z() - v1.z();
+
+float s1 = w2.x() - w1.x();
+float s2 = w3.x() - w1.x();
+float t1 = w2.y() - w1.y();
+float t2 = w3.y() - w1.y();
+
+float r = 1.0f / (s1 * t2 - s2 * t1);
+Vector3 sdir((t2 * x1 - t1 * x2) * r, (t2 * y1 - t1 * y2) * r, (t2 * z1 - t1 * z2) * r);
+Vector3 tdir((s1 * x2 - s2 * x1) * r, (s1 * y2 - s2 * y1) * r, (s1 * z2 - s2 * z1) * r);
+
+Vector3 N = cross(v3-v1, v2-v1);
+
+// Rotate 90 around N.
+}
+*/
+
+static float triangleConformalEnergy(Vector3 q[3], Vector2 p[3])
+{
+    // Using Denis formulas:
+    Vector3 c0 = q[1] - q[2];
+    Vector3 c1 = q[2] - q[0];
+    Vector3 c2 = q[0] - q[1];
+
+    Vector3 N = cross(-c0, c1);
+    float T = length(N);	// 2T
+    N = normalize(N, 0);
+
+    float cot_alpha0 = dot(-c1, c2) / length(cross(-c1, c2));
+    float cot_alpha1 = dot(-c2, c0) / length(cross(-c2, c0));
+    float cot_alpha2 = dot(-c0, c1) / length(cross(-c0, c1));
+
+    Vector3 t0 = -cot_alpha1 * c1 + cot_alpha2 * c2;
+    Vector3 t1 = -cot_alpha2 * c2 + cot_alpha0 * c0;
+    Vector3 t2 = -cot_alpha0 * c0 + cot_alpha1 * c1;
+
+    nvCheck(equal(length(t0), length(c0)));
+    nvCheck(equal(length(t1), length(c1)));
+    nvCheck(equal(length(t2), length(c2)));
+    nvCheck(equal(dot(t0, c0), 0));
+    nvCheck(equal(dot(t1, c1), 0));
+    nvCheck(equal(dot(t2, c2), 0));
+
+    // Gradients
+    Vector3 grad_u = 1.0f / T * (p[0].x * t0 + p[1].x * t1 + p[2].x * t2);
+    Vector3 grad_v = 1.0f / T * (p[0].y * t0 + p[1].y * t1 + p[2].y * t2);
+
+    // Rotated gradients
+    Vector3 Jgrad_u = 1.0f / T * (p[0].x * c0 + p[1].x * c1 + p[2].x * c2);
+    Vector3 Jgrad_v = 1.0f / T * (p[0].y * c0 + p[1].y * c1 + p[2].y * c2);
+
+    // Using Lengyel's formulas:
+    { 
+        const Vector3 v1 = q[0];
+        const Vector3 v2 = q[1];
+        const Vector3 v3 = q[2];
+
+        const Vector2 w1 = p[0];
+        const Vector2 w2 = p[1];
+        const Vector2 w3 = p[2];
+
+        float x1 = v2.x - v1.x;
+        float x2 = v3.x - v1.x;
+        float y1 = v2.y - v1.y;
+        float y2 = v3.y - v1.y;
+        float z1 = v2.z - v1.z;
+        float z2 = v3.z - v1.z;
+
+        float s1 = w2.x - w1.x;
+        float s2 = w3.x - w1.x;
+        float t1 = w2.y - w1.y;
+        float t2 = w3.y - w1.y;
+
+        float r = 1.0f / (s1 * t2 - s2 * t1);
+        Vector3 sdir((t2 * x1 - t1 * x2) * r, (t2 * y1 - t1 * y2) * r, (t2 * z1 - t1 * z2) * r);
+        Vector3 tdir((s1 * x2 - s2 * x1) * r, (s1 * y2 - s2 * y1) * r, (s1 * z2 - s2 * z1) * r);
+
+        Vector3 Jsdir = cross(N, sdir);
+        Vector3 Jtdir = cross(N, tdir);
+
+        float x = 3;
+    }
+
+    // check: sdir == grad_u
+    // check: tdir == grad_v
+
+    return length(grad_u - Jgrad_v);
+}
+#endif // 0
+
+
+ParameterizationQuality::ParameterizationQuality()
+{
+    m_totalTriangleCount = 0;
+    m_flippedTriangleCount = 0;
+    m_zeroAreaTriangleCount = 0;
+
+    m_parametricArea = 0.0f;
+    m_geometricArea = 0.0f;
+
+    m_stretchMetric = 0.0f;
+    m_maxStretchMetric = 0.0f;
+
+    m_conformalMetric = 0.0f;
+    m_authalicMetric = 0.0f;
+}
+
+ParameterizationQuality::ParameterizationQuality(const HalfEdge::Mesh * mesh)
+{
+    nvDebugCheck(mesh != NULL);
+
+    m_totalTriangleCount = 0;
+    m_flippedTriangleCount = 0;
+    m_zeroAreaTriangleCount = 0;
+
+    m_parametricArea = 0.0f;
+    m_geometricArea = 0.0f;
+
+    m_stretchMetric = 0.0f;
+    m_maxStretchMetric = 0.0f;
+
+    m_conformalMetric = 0.0f;
+    m_authalicMetric = 0.0f;
+
+    const uint faceCount = mesh->faceCount();
+    for (uint f = 0; f < faceCount; f++)
+    {
+        const HalfEdge::Face * face = mesh->faceAt(f);
+        const HalfEdge::Vertex * vertex0 = NULL;
+
+        Vector3 p[3];
+        Vector2 t[3];
+
+        for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance())
+        {
+            const HalfEdge::Edge * edge = it.current();
+
+            if (vertex0 == NULL)
+            {
+                vertex0 = edge->vertex;
+
+                p[0] = vertex0->pos;
+                t[0] = vertex0->tex;
+            }
+            else if (edge->to() != vertex0)
+            {
+                p[1] = edge->from()->pos;
+                p[2] = edge->to()->pos;
+                t[1] = edge->from()->tex;
+                t[2] = edge->to()->tex;
+
+                processTriangle(p, t);
+            }
+        }
+    }
+
+    if (m_flippedTriangleCount + m_zeroAreaTriangleCount == faceCount)
+    {
+        // If all triangles are flipped, then none is.
+        m_flippedTriangleCount = 0;
+    }
+
+    nvDebugCheck(isFinite(m_parametricArea) && m_parametricArea >= 0);
+    nvDebugCheck(isFinite(m_geometricArea) && m_geometricArea >= 0);
+    nvDebugCheck(isFinite(m_stretchMetric));
+    nvDebugCheck(isFinite(m_maxStretchMetric));
+    nvDebugCheck(isFinite(m_conformalMetric));
+    nvDebugCheck(isFinite(m_authalicMetric));
+}
+
+bool ParameterizationQuality::isValid() const
+{
+    return m_flippedTriangleCount == 0; // @@ Does not test for self-overlaps.
+}
+
+float ParameterizationQuality::rmsStretchMetric() const
+{
+    if (m_geometricArea == 0) return 0.0f;
+    float normFactor = sqrtf(m_parametricArea / m_geometricArea);
+    return sqrtf(m_stretchMetric / m_geometricArea) * normFactor;
+}
+
+float ParameterizationQuality::maxStretchMetric() const
+{
+    if (m_geometricArea == 0) return 0.0f;
+    float normFactor = sqrtf(m_parametricArea / m_geometricArea);
+    return m_maxStretchMetric * normFactor;
+}
+
+float ParameterizationQuality::rmsConformalMetric() const
+{
+    if (m_geometricArea == 0) return 0.0f;
+    return sqrtf(m_conformalMetric / m_geometricArea);
+}
+
+float ParameterizationQuality::maxAuthalicMetric() const
+{
+    if (m_geometricArea == 0) return 0.0f;
+    return sqrtf(m_authalicMetric / m_geometricArea);
+}
+
+void ParameterizationQuality::operator += (const ParameterizationQuality & pq)
+{
+    m_totalTriangleCount += pq.m_totalTriangleCount;
+    m_flippedTriangleCount += pq.m_flippedTriangleCount;
+    m_zeroAreaTriangleCount += pq.m_zeroAreaTriangleCount;
+
+    m_parametricArea += pq.m_parametricArea;
+    m_geometricArea += pq.m_geometricArea;
+
+    m_stretchMetric += pq.m_stretchMetric;
+    m_maxStretchMetric = max(m_maxStretchMetric, pq.m_maxStretchMetric);
+
+    m_conformalMetric += pq.m_conformalMetric;
+    m_authalicMetric += pq.m_authalicMetric;
+}
+
+
+void ParameterizationQuality::processTriangle(Vector3 q[3], Vector2 p[3])
+{
+    m_totalTriangleCount++;
+
+    // Evaluate texture stretch metric. See:
+    // - "Texture Mapping Progressive Meshes", Sander, Snyder, Gortler & Hoppe
+    // - "Mesh Parameterization: Theory and Practice", Siggraph'07 Course Notes, Hormann, Levy & Sheffer.
+
+    float t1 = p[0].x;
+    float s1 = p[0].y;
+    float t2 = p[1].x;
+    float s2 = p[1].y;
+    float t3 = p[2].x;
+    float s3 = p[2].y;
+
+    float geometricArea = length(cross(q[1] - q[0], q[2] - q[0])) / 2;
+    float parametricArea = ((s2 - s1)*(t3 - t1) - (s3 - s1)*(t2 - t1)) / 2;
+    
+    if (isZero(parametricArea))
+    {
+        m_zeroAreaTriangleCount++;
+        return;
+    }
+
+    Vector3 Ss = (q[0] * (t2- t3) + q[1] * (t3 - t1) + q[2] * (t1 - t2)) / (2 * parametricArea);
+    Vector3 St = (q[0] * (s3- s2) + q[1] * (s1 - s3) + q[2] * (s2 - s1)) / (2 * parametricArea);
+
+    float a = dot(Ss, Ss); // E
+    float b = dot(Ss, St); // F
+    float c = dot(St, St); // G
+
+    // Compute eigen-values of the first fundamental form:
+    float sigma1 = sqrtf(0.5f * max(0.0f, a + c - sqrtf(square(a - c) + 4 * square(b)))); // gamma uppercase, min eigenvalue.
+    float sigma2 = sqrtf(0.5f * max(0.0f, a + c + sqrtf(square(a - c) + 4 * square(b)))); // gamma lowercase, max eigenvalue.
+    nvCheck(sigma2 >= sigma1);
+
+    // isometric: sigma1 = sigma2 = 1
+    // conformal: sigma1 / sigma2 = 1
+    // authalic: sigma1 * sigma2 = 1
+
+    float rmsStretch = sqrtf((a + c) * 0.5f);
+    float rmsStretch2 = sqrtf((square(sigma1) + square(sigma2)) * 0.5f);
+    nvDebugCheck(equal(rmsStretch, rmsStretch2, 0.01f));
+
+    if (parametricArea < 0.0f)
+    {
+        // Count flipped triangles.
+        m_flippedTriangleCount++;
+
+        parametricArea = fabsf(parametricArea);
+    }
+
+    m_stretchMetric += square(rmsStretch) * geometricArea;
+    m_maxStretchMetric = max(m_maxStretchMetric, sigma2);
+
+    if (!isZero(sigma1, 0.000001f)) {
+        // sigma1 is zero when geometricArea is zero.
+        m_conformalMetric += (sigma2 / sigma1) * geometricArea;
+    }
+    m_authalicMetric += (sigma1 * sigma2) * geometricArea;
+
+    // Accumulate total areas.
+    m_geometricArea += geometricArea;
+    m_parametricArea += parametricArea;
+
+
+    //triangleConformalEnergy(q, p);
+}