diff --git a/src/Triangle.Examples/Examples/Example1.cs b/src/Triangle.Examples/Examples/Example1.cs index 5419d48..82d09d5 100644 --- a/src/Triangle.Examples/Examples/Example1.cs +++ b/src/Triangle.Examples/Examples/Example1.cs @@ -27,7 +27,7 @@ namespace TriangleNet.Examples if (print) SvgImage.Save(mesh, "example-1.svg", 500); - return true; + return mesh.Triangles.Count > 0; } } } diff --git a/src/Triangle.Examples/Examples/Example10.cs b/src/Triangle.Examples/Examples/Example10.cs new file mode 100644 index 0000000..ff13379 --- /dev/null +++ b/src/Triangle.Examples/Examples/Example10.cs @@ -0,0 +1,125 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using TriangleNet.Geometry; +using TriangleNet.Meshing; +using TriangleNet.Meshing.Algorithm; +using TriangleNet.Rendering.Text; +using TriangleNet.Tools; + +namespace TriangleNet.Examples +{ + /// + /// Scattered data interpolation without USE_Z or USE_ATTRIBS. + /// + internal class Example10 + { + // The function we are sampling. + private static readonly Func F = p => Math.Sin(p.X) * Math.Cos(p.Y); + + public static bool Run(bool print = false) + { + // The input domain. + var r = new Rectangle(0d, 0d, 10d, 10d); + + var mesh = GetScatteredDataMesh(r, out double[] data); + //var mesh = GetStructuredDataMesh(r, out double[] data); + + if (print) SvgImage.Save(mesh, "example-10.svg", 500); + + // The points to interpolate. + var xy = Generate.RandomPoints(50, r); + + var xyData = InterpolateData((Mesh)mesh, data, xy); + + double error = xy.Max(p => Math.Abs(xyData[p.ID] - F(p))); + + // L2 error + // double error = Math.Sqrt(xy.Sum(p => Math.Pow(xyData[p.ID] - F(p), 2))); + + return error < 0.5; + } + + private static IMesh GetStructuredDataMesh(Rectangle domain, out double[] data) + { + var mesh = GenericMesher.StructuredMesh(domain, 20, 20); + + mesh.Renumber(); + + // Generate function values for mesh points. + data = new double[mesh.Vertices.Count]; + + foreach (var item in mesh.Vertices) + { + data[item.ID] = F(item); + } + + return mesh; + } + + private static IMesh GetScatteredDataMesh(Rectangle domain, out double[] data) + { + var r = new Rectangle(domain); + + double h = domain.Width / 20; + + // Generate a rectangle boundary point set (20 points on each side). + var input = Generate.Rectangle(r, 0.5); + + // Making sure we add some margin to the boundary. + h = -h / 2; + r.Resize(h, h); + + // Add more input points (more sampling points, better interpolation). + input.Points.AddRange(Generate.RandomPoints(350, r)); + + var mesher = new GenericMesher(new Dwyer()); + + // Generate mesh. + var mesh = mesher.Triangulate(input.Points); + + mesh.Renumber(); + + // Generate function values for mesh points. + data = new double[mesh.Vertices.Count]; + + foreach (var item in mesh.Vertices) + { + data[item.ID] = F(item); + } + + return mesh; + } + + private static double[] InterpolateData(Mesh mesh, double[] data, IEnumerable xy) + { + // The interpolated values. + var values = new double[xy.Count()]; + + var qtree = new TriangleQuadTree(mesh); + + int i = 0; + + foreach (var p in xy) + { + var tri = qtree.Query(p.X, p.Y); + + // For easy access of the interpolated values. + p.ID = i; + + if (tri == null) + { + values[i] = float.NaN; + } + else + { + values[i] = Interpolation.InterpolatePoint(tri, p, data); + } + + i++; + } + + return values; + } + } +} \ No newline at end of file diff --git a/src/Triangle.Examples/Examples/Example2.cs b/src/Triangle.Examples/Examples/Example2.cs index e3dd6f1..c497ccf 100644 --- a/src/Triangle.Examples/Examples/Example2.cs +++ b/src/Triangle.Examples/Examples/Example2.cs @@ -22,7 +22,7 @@ if (print) SvgImage.Save(mesh, "example-2.svg", 500); - return true; + return mesh.Triangles.Count > 0; } public static IPolygon CreatePolygon(double h = 0.2) diff --git a/src/Triangle.Examples/Examples/Example3.cs b/src/Triangle.Examples/Examples/Example3.cs index 74fd070..7e0182d 100644 --- a/src/Triangle.Examples/Examples/Example3.cs +++ b/src/Triangle.Examples/Examples/Example3.cs @@ -19,7 +19,7 @@ namespace TriangleNet.Examples if (print) SvgImage.Save(mesh, "example-3.svg", 500); - return true; + return mesh.Triangles.Count > 0; } public static IMesh CreateMesh() diff --git a/src/Triangle.Examples/Examples/Example4.cs b/src/Triangle.Examples/Examples/Example4.cs index cc1fff6..d3c122a 100644 --- a/src/Triangle.Examples/Examples/Example4.cs +++ b/src/Triangle.Examples/Examples/Example4.cs @@ -43,7 +43,7 @@ namespace TriangleNet.Examples if (print) SvgImage.Save(mesh, "example-4.svg", 500); - return true; + return mesh.Triangles.Count > 0; } public static IPolygon CreatePolygon() diff --git a/src/Triangle.Examples/Examples/Example5.cs b/src/Triangle.Examples/Examples/Example5.cs index ac9cbd2..17540af 100644 --- a/src/Triangle.Examples/Examples/Example5.cs +++ b/src/Triangle.Examples/Examples/Example5.cs @@ -25,7 +25,7 @@ namespace TriangleNet.Examples if (print) SvgImage.Save(mesh, "example-5-2.svg", 500, true, false); - return true; + return mesh.Triangles.Count > 0; } /// diff --git a/src/Triangle.Examples/Examples/Example6.cs b/src/Triangle.Examples/Examples/Example6.cs index 739fafa..0ed1b2c 100644 --- a/src/Triangle.Examples/Examples/Example6.cs +++ b/src/Triangle.Examples/Examples/Example6.cs @@ -4,6 +4,7 @@ using TriangleNet; using TriangleNet.Geometry; using TriangleNet.Meshing.Iterators; + using TriangleNet.Rendering.Text; using TriangleNet.Tools; using TriangleNet.Topology; @@ -12,14 +13,14 @@ /// public static class Example6 { - public static bool Run() + public static bool Run(bool print = true) { // Generate the input geometry. var polygon = new Polygon(8, true); // Two intersecting rectangles. - var A = Generate.Rectangle(0.0, 0.0, 4.0, 4.0, 1); - var B = Generate.Rectangle(1.0, 1.0, 4.0, 4.0, 2); + var A = Generate.Rectangle(0d, 0d, 4d, 4d, label: 1); + var B = Generate.Rectangle(1d, 1d, 4d, 4d, label: 2); polygon.Add(A); polygon.Add(B); @@ -27,6 +28,8 @@ // Generate mesh. var mesh = (Mesh)polygon.Triangulate(); + if (print) SvgImage.Save(mesh, "example-6.svg", 500); + // Find a seeding triangle (in this case, the point (2, 2) lies in // both rectangles). var seed = (new TriangleQuadTree(mesh)).Query(2.0, 2.0) as Triangle; @@ -47,7 +50,7 @@ // The xor of A and B. var xor = mesh.Triangles.Where(t => t.Label == 1 || t.Label == 2); - return true; + return intersection.Any() && difference.Any() && xor.Any(); } } } \ No newline at end of file diff --git a/src/Triangle.Examples/Examples/Example7.cs b/src/Triangle.Examples/Examples/Example7.cs index 06bedc6..cdc2f8f 100644 --- a/src/Triangle.Examples/Examples/Example7.cs +++ b/src/Triangle.Examples/Examples/Example7.cs @@ -37,7 +37,7 @@ namespace TriangleNet.Examples if (length > MAX_EDGE_LENGTH) { - Console.WriteLine("Something's wrong in here ..."); + return false; } } diff --git a/src/Triangle.Examples/Examples/Example8.cs b/src/Triangle.Examples/Examples/Example8.cs index 8d131ba..f0a16a5 100644 --- a/src/Triangle.Examples/Examples/Example8.cs +++ b/src/Triangle.Examples/Examples/Example8.cs @@ -16,12 +16,10 @@ namespace TriangleNet.Examples { var mesh = (Mesh)Example3.CreateMesh(); - FindAdjacencyMatrix(mesh); - - return true; + return FindAdjacencyMatrix(mesh); } - private static void FindAdjacencyMatrix(Mesh mesh) + private static bool FindAdjacencyMatrix(Mesh mesh) { mesh.Renumber(); @@ -57,8 +55,10 @@ namespace TriangleNet.Examples // Column pointers should be exactly the same. if (!CompareArray(matrix1.ColumnPointers, matrix2.ColumnPointers)) { - Console.WriteLine("Something's wrong in here ..."); + return false; } + + return true; } private static bool CompareArray(int[] a, int[] b) diff --git a/src/Triangle.Examples/Program.cs b/src/Triangle.Examples/Program.cs index e0ffef4..05269a1 100644 --- a/src/Triangle.Examples/Program.cs +++ b/src/Triangle.Examples/Program.cs @@ -13,10 +13,11 @@ namespace TriangleNet Check("Example 3", Example3.Run()); Check("Example 4", Example4.Run()); Check("Example 5", Example5.Run()); - Check("Example 6", Example8.Run()); - Check("Example 7", Example6.Run()); - Check("Example 8", Example7.Run()); + Check("Example 6", Example6.Run()); + Check("Example 7", Example7.Run()); + Check("Example 8", Example8.Run()); Check("Example 9", Example9.Run()); + Check("Example 10", Example10.Run()); } static void Check(string item, bool success) diff --git a/src/Triangle/Tools/Interpolation.cs b/src/Triangle/Tools/Interpolation.cs index 857275c..227e020 100644 --- a/src/Triangle/Tools/Interpolation.cs +++ b/src/Triangle/Tools/Interpolation.cs @@ -1,10 +1,54 @@ - +// ----------------------------------------------------------------------- +// +// Triangle Copyright (c) 1993, 1995, 1997, 1998, 2002, 2005 Jonathan Richard Shewchuk +// Triangle.NET code by Christian Woltering +// +// ----------------------------------------------------------------------- + namespace TriangleNet.Tools { using TriangleNet.Geometry; public static class Interpolation { + /// + /// Linear interpolation of a point. + /// + /// The triangle containing the point + /// The point to interpolate. + /// The vertex data (z values). + /// The linear interpolation value. + /// + /// IMPORTANT: this method assumes the mesh vertex ids correspond to the data array indices. + /// + public static double InterpolatePoint(ITriangle tri, Point p, double[] data) + { + var org = tri.GetVertex(0); + var dest = tri.GetVertex(1); + var apex = tri.GetVertex(2); + + double xdo = dest.x - org.x; + double ydo = dest.y - org.y; + double xao = apex.x - org.x; + double yao = apex.y - org.y; + + double denominator = 0.5 / (xdo * yao - xao * ydo); + + double dx = p.x - org.x; + double dy = p.y - org.y; + + // To interpolate z value for the given point inserted, define a + // coordinate system with a xi-axis, directed from the triangle's + // origin to its destination, and an eta-axis, directed from its + // origin to its apex. + double xi = (yao * dx - xao * dy) * (2.0 * denominator); + double eta = (xdo * dy - ydo * dx) * (2.0 * denominator); + + double orgz = data[org.id]; + + return orgz + xi * (data[dest.id] - orgz) + eta * (data[apex.id] - orgz); + } + #if USE_ATTRIBS /// /// Linear interpolation of vertex attributes. @@ -17,35 +61,26 @@ namespace TriangleNet.Tools /// public static void InterpolateAttributes(Vertex vertex, ITriangle triangle, int n) { - Vertex org = triangle.GetVertex(0); - Vertex dest = triangle.GetVertex(1); - Vertex apex = triangle.GetVertex(2); + var org = triangle.GetVertex(0); + var dest = triangle.GetVertex(1); + var apex = triangle.GetVertex(2); - double xdo, ydo, xao, yao; - double denominator; - double dx, dy; - double xi, eta; + double xdo = dest.x - org.x; + double ydo = dest.y - org.y; + double xao = apex.x - org.x; + double yao = apex.y - org.y; - // Compute the circumcenter of the triangle. - xdo = dest.x - org.x; - ydo = dest.y - org.y; - xao = apex.x - org.x; - yao = apex.y - org.y; + double denominator = 0.5 / (xdo * yao - xao * ydo); - denominator = 0.5 / (xdo * yao - xao * ydo); - - //dx = (yao * dodist - ydo * aodist) * denominator; - //dy = (xdo * aodist - xao * dodist) * denominator; - - dx = vertex.x - org.x; - dy = vertex.y - org.y; + double dx = vertex.x - org.x; + double dy = vertex.y - org.y; // To interpolate vertex attributes for the new vertex, define a // coordinate system with a xi-axis directed from the triangle's // origin to its destination, and an eta-axis, directed from its // origin to its apex. - xi = (yao * dx - xao * dy) * (2.0 * denominator); - eta = (xdo * dy - ydo * dx) * (2.0 * denominator); + double xi = (yao * dx - xao * dy) * (2.0 * denominator); + double eta = (xdo * dy - ydo * dx) * (2.0 * denominator); for (int i = 0; i < n; i++) { @@ -68,35 +103,27 @@ namespace TriangleNet.Tools /// public static void InterpolateZ(Point p, ITriangle triangle) { - Vertex org = triangle.GetVertex(0); - Vertex dest = triangle.GetVertex(1); - Vertex apex = triangle.GetVertex(2); - - double xdo, ydo, xao, yao; - double denominator; - double dx, dy; - double xi, eta; + var org = triangle.GetVertex(0); + var dest = triangle.GetVertex(1); + var apex = triangle.GetVertex(2); // Compute the circumcenter of the triangle. - xdo = dest.x - org.x; - ydo = dest.y - org.y; - xao = apex.x - org.x; - yao = apex.y - org.y; + double xdo = dest.x - org.x; + double ydo = dest.y - org.y; + double xao = apex.x - org.x; + double yao = apex.y - org.y; - denominator = 0.5 / (xdo * yao - xao * ydo); + double denominator = 0.5 / (xdo * yao - xao * ydo); - //dx = (yao * dodist - ydo * aodist) * denominator; - //dy = (xdo * aodist - xao * dodist) * denominator; - - dx = p.x - org.x; - dy = p.y - org.y; + double dx = p.x - org.x; + double dy = p.y - org.y; // To interpolate z value for the given point inserted, define a // coordinate system with a xi-axis, directed from the triangle's // origin to its destination, and an eta-axis, directed from its // origin to its apex. - xi = (yao * dx - xao * dy) * (2.0 * denominator); - eta = (xdo * dy - ydo * dx) * (2.0 * denominator); + double xi = (yao * dx - xao * dy) * (2.0 * denominator); + double eta = (xdo * dy - ydo * dx) * (2.0 * denominator); p.z = org.z + xi * (dest.z - org.z) + eta * (apex.z - org.z); }