diff --git a/src/Triangle.Examples/Examples/Example4.cs b/src/Triangle.Examples/Examples/Example4.cs index 270054d..8f43295 100644 --- a/src/Triangle.Examples/Examples/Example4.cs +++ b/src/Triangle.Examples/Examples/Example4.cs @@ -21,19 +21,21 @@ namespace TriangleNet.Examples // The ideal area if triangles were equilateral. var area = Math.Sqrt(3) / 4 * h * h; - var quality = new QualityMeasure(mesh); + var quality = new QualityMeasure(); + + quality.Update(mesh); if (print) { Console.WriteLine($" Ideal area: {area}"); - Console.WriteLine($" Min. area: {quality.AreaMinimum}"); - Console.WriteLine($" Max. area: {quality.AreaMaximum}"); + Console.WriteLine($" Min. area: {quality.Area.Minimum}"); + Console.WriteLine($" Max. area: {quality.Area.Maximum}"); Console.WriteLine($" Avg. area: {quality.AreaTotal / mesh.Triangles.Count}"); SvgImage.Save(mesh, "example-4.svg", 500); } - return quality.AreaMinimum < area && quality.AreaMaximum > area; + return quality.Area.Minimum < area && quality.Area.Maximum > area; } // The boundary segment size of the input geometry. diff --git a/src/Triangle.Tests/Tools/QualityMeasureTest.cs b/src/Triangle.Tests/Tools/QualityMeasureTest.cs index b43e6a4..46f1ed0 100644 --- a/src/Triangle.Tests/Tools/QualityMeasureTest.cs +++ b/src/Triangle.Tests/Tools/QualityMeasureTest.cs @@ -24,13 +24,17 @@ namespace TriangleNet.Tests.Tools var mesh = new Dwyer().Triangulate(vertices, new Configuration()); - var quality = new QualityMeasure(mesh); + var quality = new QualityMeasure(); - Assert.AreEqual(quality.AreaMinimum, quality.AreaMaximum); - Assert.AreEqual(1.0, quality.AlphaMinimum); - Assert.AreEqual(1.0, quality.AlphaMaximum); - Assert.AreEqual(1.0, quality.Q_Minimum); - Assert.AreEqual(1.0, quality.Q_Maximum); + quality.Update(mesh); + + Assert.AreEqual(quality.Area.Minimum, quality.Area.Maximum); + Assert.AreEqual(1.0, quality.Alpha.Minimum); + Assert.AreEqual(1.0, quality.Alpha.Maximum); + Assert.AreEqual(1.0, quality.Eta.Minimum); + Assert.AreEqual(1.0, quality.Eta.Maximum); + Assert.AreEqual(1.0, quality.Q.Minimum); + Assert.AreEqual(1.0, quality.Q.Maximum); } } } diff --git a/src/Triangle.Viewer/Views/StatisticView.cs b/src/Triangle.Viewer/Views/StatisticView.cs index ac6fdaf..4f6730b 100644 --- a/src/Triangle.Viewer/Views/StatisticView.cs +++ b/src/Triangle.Viewer/Views/StatisticView.cs @@ -89,13 +89,15 @@ namespace MeshExplorer.Views lbAngleMax.Text = Util.AngleToString(statistic.LargestAngle); // Update quality - var quality = new QualityMeasure(mesh); + var quality = new QualityMeasure(); - lbQualAlphaMin.Text = Util.DoubleToString(quality.AlphaMinimum); - lbQualAlphaAve.Text = Util.DoubleToString(quality.AlphaAverage); + quality.Update(mesh); - lbQualAspectMin.Text = Util.DoubleToString(quality.Q_Minimum); - lbQualAspectAve.Text = Util.DoubleToString(quality.Q_Average); + lbQualAlphaMin.Text = Util.DoubleToString(quality.Alpha.Minimum); + lbQualAlphaAve.Text = Util.DoubleToString(quality.Alpha.Average); + + lbQualAspectMin.Text = Util.DoubleToString(quality.Q.Minimum); + lbQualAspectAve.Text = Util.DoubleToString(quality.Q.Average); } #endregion diff --git a/src/Triangle/Geometry/Vertex.cs b/src/Triangle/Geometry/Vertex.cs index 8b3a1c3..e84b9e4 100644 --- a/src/Triangle/Geometry/Vertex.cs +++ b/src/Triangle/Geometry/Vertex.cs @@ -53,6 +53,19 @@ namespace TriangleNet.Geometry type = VertexType.InputVertex; } + /// + /// Initializes a new instance of the class. + /// + /// The x coordinate of the vertex. + /// The y coordinate of the vertex. + /// The boundary mark. + /// The vertex type. + public Vertex(double x, double y, int mark, VertexType type) + : base(x, y, mark) + { + this.type = type; + } + #if USE_ATTRIBS /// /// Initializes a new instance of the class. diff --git a/src/Triangle/Tools/QualityMeasure.cs b/src/Triangle/Tools/QualityMeasure.cs index 08780fb..dcf5144 100644 --- a/src/Triangle/Tools/QualityMeasure.cs +++ b/src/Triangle/Tools/QualityMeasure.cs @@ -1,6 +1,5 @@ // ----------------------------------------------------------------------- // -// Original Matlab code by John Burkardt, Florida State University // Triangle.NET code by Christian Woltering // // ----------------------------------------------------------------------- @@ -8,9 +7,70 @@ namespace TriangleNet.Tools { using System; + using System.Collections.Generic; using TriangleNet.Geometry; using TriangleNet.Meshing; + /// + /// Base class for quality measures. + /// + public abstract class Measure + { + protected double min, max, avg, are; + + /// + /// Gets the minimum value of the measure. + /// + public double Minimum => min; + + /// + /// Gets the maximum value of the measure. + /// + public double Maximum => max; + + /// + /// Gets the average value of the measure. + /// + public double Average => avg; + + /// + /// Gets the value averaged over all triangles and weighted by area. + /// + public double Area => are; + + /// + /// Initialize the measure. + /// + public virtual void Initialize() + { + min = double.MaxValue; + max = -double.MaxValue; + avg = 0d; + are = 0d; + } + + /// + /// Compute measure of given triangle. + /// + /// Side length ab. + /// Side length bc. + /// Side length ca. + /// Triangle area. + public abstract double Update(double ab, double bc, double ca, double area); + + /// + /// Finalize the measure. + /// + /// Total number of triangles measured. + /// Total area of triangles measured. + public virtual void Finalize(int n, double totalArea) + { + avg = n > 0 ? avg / n : avg; + + are = totalArea > 0.0 ? are / totalArea : are; + } + } + /// /// Provides mesh quality information. /// @@ -54,84 +114,95 @@ namespace TriangleNet.Tools /// public class QualityMeasure { - AreaMeasure areaMeasure; - AlphaMeasure alphaMeasure; - Q_Measure qMeasure; - + MeasureArea _area; + MeasureAlpha _alpha; + MeasureEta _eta; + MeasureQ _q; #region Public properties /// - /// Minimum triangle area. + /// Gets the area measure. /// - public double AreaMinimum => areaMeasure.area_min; + public Measure Area => _area; /// - /// Maximum triangle area. + /// Gets the alpha measure. /// - public double AreaMaximum => areaMeasure.area_max; + /// + /// The alpha measure computes the minimum angle among all triangles. + /// The best possible value is 1, and the worst 0. + /// + public Measure Alpha => _alpha; /// - /// Total triangulation area. + /// Gets the eta measure. /// - public double AreaTotal => areaMeasure.area_total; + /// + /// The eta measure relates the area of a triangle a to its edge lengths. + /// The best possible value is 1, and the worst 0. + /// + public Measure Eta => _eta; /// - /// Ratio of maximum and minimum triangle area. + /// Gets the Q measure, also knows as normalized shape ratio (NSR). /// - public double AreaRatio => areaMeasure.area_max / areaMeasure.area_min; + /// + /// The Q measure relates the incircle to the circumcircle radius. + /// In an ideally regular mesh, all triangles would have the same + /// equilateral shape, for which Q = 1. + /// + public Measure Q => _q; /// - /// Smallest angle. + /// Gets the total triangulation area. /// - public double AlphaMinimum => alphaMeasure.alpha_min; - - /// - /// Maximum smallest angle. - /// - public double AlphaMaximum => alphaMeasure.alpha_max; - - /// - /// Average angle. - /// - public double AlphaAverage => alphaMeasure.alpha_ave; - - /// - /// Average angle weighted by area. - /// - public double AlphaArea => alphaMeasure.alpha_area; - - /// - /// Smallest aspect ratio. - /// - public double Q_Minimum => qMeasure.q_min; - - /// - /// Largest aspect ratio. - /// - public double Q_Maximum => qMeasure.q_max; - - /// - /// Average aspect ratio. - /// - public double Q_Average => qMeasure.q_ave; - - /// - /// Average aspect ratio weighted by area. - /// - public double Q_Area => qMeasure.q_area; + public double AreaTotal => _area.Area; #endregion + List measures; + /// /// Initializes a new instance of the class. /// - public QualityMeasure(IMesh mesh) + public QualityMeasure() { - areaMeasure = new AreaMeasure(); - alphaMeasure = new AlphaMeasure(); - qMeasure = new Q_Measure(); + _area = new MeasureArea(); + _alpha = new MeasureAlpha(); + _eta = new MeasureEta(); + _q = new MeasureQ(); + measures = new List() + { + _area, _alpha, _eta, _q + }; + } + + /// + /// Add a custom measure. + /// + /// + public void Add(Measure measure) + { + measures.Add(measure); + } + + /// + /// Update all measures for the given mesh. + /// + /// The mesh. + public void Update(IMesh mesh) + { + Update(mesh.Triangles); + } + + /// + /// Update all measures for the given triangles. + /// + /// The triangles. + public void Update(IEnumerable triangles) + { Point a, b, c; double ab, bc, ca; double lx, ly; @@ -139,13 +210,18 @@ namespace TriangleNet.Tools int n = 0; - foreach (var tri in mesh.Triangles) + foreach (var m in measures) + { + m.Initialize(); + } + + foreach (var tri in triangles) { n++; - a = tri.vertices[0]; - b = tri.vertices[1]; - c = tri.vertices[2]; + a = tri.GetVertex(0); + b = tri.GetVertex(1); + c = tri.GetVertex(2); lx = a.x - b.x; ly = a.y - b.y; @@ -157,14 +233,20 @@ namespace TriangleNet.Tools ly = c.y - a.y; ca = Math.Sqrt(lx * lx + ly * ly); - area = areaMeasure.Measure(a, b, c); - alphaMeasure.Measure(ab, bc, ca, area); - qMeasure.Measure(ab, bc, ca, area); + area = 0.5 * Math.Abs(a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y)); + + foreach (var m in measures) + { + m.Update(ab, bc, ca, area); + } } - // Normalize measures - alphaMeasure.Normalize(n, areaMeasure.area_total); - qMeasure.Normalize(n, areaMeasure.area_total); + var totalArea = _area.Area; + + foreach (var m in measures) + { + m.Finalize(n, totalArea); + } } /// @@ -194,6 +276,8 @@ namespace TriangleNet.Tools /// /// Because the finite element node adjacency relationship is symmetric, /// we are guaranteed that ML = MU. + /// + /// Based on Matlab code by John Burkardt, Florida State University. /// public static int Bandwidth(IMesh mesh) { @@ -223,73 +307,60 @@ namespace TriangleNet.Tools return ml + 1 + mu; } - class AreaMeasure + class MeasureArea : Measure { - // Minimum area - public double area_min = double.MaxValue; - // Maximum area - public double area_max = -double.MaxValue; - // Total area of geometry - public double area_total = 0; + private readonly double EPS = 0.0; + // Number of triangles with zero area - public int area_zero = 0; + public int zero; - /// - /// Compute the area of given triangle. - /// - /// Triangle corner a. - /// Triangle corner b. - /// Triangle corner c. - /// Triangle area. - public double Measure(Point a, Point b, Point c) + /// + public override void Initialize() { - double area = 0.5 * Math.Abs(a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y)); + zero = 0; - area_min = Math.Min(area_min, area); - area_max = Math.Max(area_max, area); - area_total += area; + base.Initialize(); + } - if (area == 0.0) + /// + public override double Update(double ab, double bc, double ca, double area) + { + min = Math.Min(min, area); + max = Math.Max(max, area); + + are += area; + + if (area <= EPS) { - area_zero++; + zero++; } return area; } + + /// + public override void Finalize(int n, double area_total) + { + avg = n > 0 ? are / n : are; + } } /// - /// The alpha measure determines the triangulated pointset quality. + /// The alpha measure determines the triangulated point set quality. /// /// /// The alpha measure evaluates the uniformity of the shapes of the triangles - /// defined by a triangulated pointset. + /// defined by a triangulated point set. /// /// We compute the minimum angle among all the triangles in the triangulated /// dataset and divide by the maximum possible value (which, in degrees, /// is 60). The best possible value is 1, and the worst 0. A good /// triangulation should have an alpha score close to 1. /// - class AlphaMeasure + class MeasureAlpha : Measure { - // Minimum value over all triangles - public double alpha_min = double.MaxValue; - // Maximum value over all triangles - public double alpha_max = -double.MaxValue; - // Value averaged over all triangles - public double alpha_ave = 0; - // Value averaged over all triangles and weighted by area - public double alpha_area = 0; - - /// - /// Compute q value of given triangle. - /// - /// Side length ab. - /// Side length bc. - /// Side length ca. - /// Triangle area. - /// - public double Measure(double ab, double bc, double ca, double area) + /// + public override double Update(double ab, double bc, double ca, double area) { double alpha = double.MaxValue; @@ -345,48 +416,60 @@ namespace TriangleNet.Tools // Normalize angle from [0,pi/3] radians into qualities in [0,1]. alpha = alpha * 3.0 / Math.PI; - alpha_ave += alpha; - alpha_area += area * alpha; + avg += alpha; + are += area * alpha; - alpha_min = Math.Min(alpha, alpha_min); - alpha_max = Math.Max(alpha, alpha_max); + min = Math.Min(alpha, min); + max = Math.Max(alpha, max); return alpha; - double acos(double c) => (c <= -1.0) ? Math.PI : ((1.0 <= c) ? 0.0 : Math.Acos(c)); - } - - /// - /// Normalize values. - /// - public void Normalize(int n, double area_total) - { - if (n > 0) - { - alpha_ave /= n; - } - else - { - alpha_ave = 0.0; - } - - if (0.0 < area_total) - { - alpha_area /= area_total; - } - else - { - alpha_area = 0.0; - } + double acos(double x) => (x <= -1.0) ? Math.PI : ((1.0 <= x) ? 0.0 : Math.Acos(x)); } } /// - /// The Q measure determines the triangulated pointset quality. + /// The eta measure determines the triangulated point set quality. + /// + /// + /// The eta measure evaluates the uniformity of the shapes of the triangles + /// defined by a triangulated point set. + /// + /// The measure relates the area of a triangle a to its edge lengths. The best + /// possible value is 1, and the worst 0. A good triangulation should have an + /// eta score close to 1. + /// + class MeasureEta : Measure + { + // Normalization factor to ensure that a perfect triangle + // (equal edges) has a quality of 1. + private readonly static double Factor = 4d * Math.Sqrt(3); + + /// + public override double Update(double ab, double bc, double ca, double area) + { + double ab2 = ab * ab; + double bc2 = bc * bc; + double ca2 = ca * ca; + + double eta = Factor * area / (ab2 + bc2 + ca2); + + avg += eta; + are += area * eta; + + min = Math.Min(eta, min); + max = Math.Max(eta, max); + + return eta; + } + } + + /// + /// The Q measure determines the triangulated point set quality. /// /// /// The Q measure evaluates the uniformity of the shapes of the triangles - /// defined by a triangulated pointset. It uses the aspect ratio + /// defined by a triangulated point set. It uses the aspect ratio /// /// 2 * (incircle radius) / (circumcircle radius) /// @@ -394,61 +477,21 @@ namespace TriangleNet.Tools /// equilateral shape, for which Q = 1. A good mesh would have /// 0.5 < Q. /// - class Q_Measure + class MeasureQ : Measure { - // Minimum value over all triangles - public double q_min = double.MaxValue; - // Maximum value over all triangles - public double q_max = -double.MaxValue; - // Average value - public double q_ave = 0; - // Average value weighted by the area of each triangle - public double q_area = 0; - - /// - /// Compute q value of given triangle. - /// - /// Side length ab. - /// Side length bc. - /// Side length ca. - /// Triangle area. - /// - public double Measure(double ab, double bc, double ca, double area) + /// + public override double Update(double ab, double bc, double ca, double area) { double q = (bc + ca - ab) * (ca + ab - bc) * (ab + bc - ca) / (ab * bc * ca); - q_min = Math.Min(q_min, q); - q_max = Math.Max(q_max, q); + min = Math.Min(min, q); + max = Math.Max(max, q); - q_ave += q; - q_area += q * area; + avg += q; + are += q * area; return q; } - - /// - /// Normalize values. - /// - public void Normalize(int n, double area_total) - { - if (n > 0) - { - q_ave /= n; - } - else - { - q_ave = 0.0; - } - - if (area_total > 0.0) - { - q_area /= area_total; - } - else - { - q_area = 0.0; - } - } } } } diff --git a/src/Triangle/Topology/DCEL/Face.cs b/src/Triangle/Topology/DCEL/Face.cs index c0c2e7e..db4a10b 100644 --- a/src/Triangle/Topology/DCEL/Face.cs +++ b/src/Triangle/Topology/DCEL/Face.cs @@ -38,6 +38,12 @@ namespace TriangleNet.Topology.DCEL internal HalfEdge edge; internal bool bounded; + /// + /// If part of a Voronoi diagram, returns the generator vertex + /// of the face. Otherwise null. + /// + public Point Generator => generator; + /// /// Gets or sets the face id. /// @@ -83,11 +89,12 @@ namespace TriangleNet.Topology.DCEL { this.generator = generator; this.edge = edge; - this.bounded = true; + + bounded = true; if (generator != null) { - this.id = generator.ID; + id = generator.ID; } } @@ -97,7 +104,7 @@ namespace TriangleNet.Topology.DCEL /// public IEnumerable EnumerateEdges() { - var edge = this.Edge; + var edge = Edge; int first = edge.ID; do