huanld/Triangle.NET
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Update examples (add interpolation example 10).

Browse Source
This commit is contained in:
wo80
2022-03-05 21:13:52 +01:00
parent ffbe041f7b
commit fd205fbb9c
11 changed files with 216 additions and 60 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/Triangle.Examples/Examples/Example1.cs
+1 -1
View File
@@ -27,7 +27,7 @@ namespace TriangleNet.Examples
if (print) SvgImage.Save(mesh, "example-1.svg", 500);
return true;
return mesh.Triangles.Count > 0;
}
}
}
src/Triangle.Examples/Examples/Example10.cs
+125
View File
@@ -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
{
/// <summary>
/// Scattered data interpolation without USE_Z or USE_ATTRIBS.
/// </summary>
internal class Example10
{
// The function we are sampling.
private static readonly Func<Point, double> 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<Point> 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;
}
}
}
src/Triangle.Examples/Examples/Example2.cs
+1 -1
View File
@@ -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)
src/Triangle.Examples/Examples/Example3.cs
+1 -1
View File
@@ -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()
src/Triangle.Examples/Examples/Example4.cs
+1 -1
View File
@@ -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()
src/Triangle.Examples/Examples/Example5.cs
+1 -1
View File
@@ -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;
}
/// <summary>
src/Triangle.Examples/Examples/Example6.cs
+7 -4
View File
@@ -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 @@
/// </summary>
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();
}
}
}
src/Triangle.Examples/Examples/Example7.cs
+1 -1
View File
@@ -37,7 +37,7 @@ namespace TriangleNet.Examples
if (length > MAX_EDGE_LENGTH)
{
Console.WriteLine("Something's wrong in here ...");
return false;
}
}
src/Triangle.Examples/Examples/Example8.cs
+5 -5
View File
@@ -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)
src/Triangle.Examples/Program.cs
+4 -3
View File
@@ -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)
src/Triangle/Tools/Interpolation.cs
+69 -42
View File
@@ -1,10 +1,54 @@
// -----------------------------------------------------------------------
// <copyright file="Interpolation.cs">
// Triangle Copyright (c) 1993, 1995, 1997, 1998, 2002, 2005 Jonathan Richard Shewchuk
// Triangle.NET code by Christian Woltering
// </copyright>
// -----------------------------------------------------------------------
namespace TriangleNet.Tools
{
using TriangleNet.Geometry;
public static class Interpolation
{
/// <summary>
/// Linear interpolation of a point.
/// </summary>
/// <param name="tri">The triangle containing the point <paramref name="p"/></param>
/// <param name="p">The point to interpolate.</param>
/// <param name="data">The vertex data (z values).</param>
/// <returns>The linear interpolation value.</returns>
/// <remarks>
/// IMPORTANT: this method assumes the mesh vertex ids correspond to the data array indices.
/// </remarks>
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
/// <summary>
/// Linear interpolation of vertex attributes.
@@ -17,35 +61,26 @@ namespace TriangleNet.Tools
/// </remarks>
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
/// </remarks>
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);
}