Added (bounded) Voronoi diagram class
git-svn-id: https://triangle.svn.codeplex.com/svn@67930 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5
This commit is contained in:
SND\wo80_cp
@@ -12,6 +12,7 @@ namespace MeshExplorer.Controls
|
||||
using System.Text;
|
||||
using System.Windows.Forms;
|
||||
using System.Drawing;
|
||||
using System.Drawing.Drawing2D;
|
||||
|
||||
/// <summary>
|
||||
/// TODO: Update summary.
|
||||
@@ -25,6 +26,24 @@ namespace MeshExplorer.Controls
|
||||
base.OnRenderItemText(e);
|
||||
}
|
||||
|
||||
protected override void OnRenderItemCheck(ToolStripItemImageRenderEventArgs e)
|
||||
{
|
||||
int h = e.Item.Height;
|
||||
int size = h / 2 - 2; // Box size
|
||||
Rectangle rect = new Rectangle(10, (h - size) / 2, size - 2, size);
|
||||
|
||||
var mode = e.Graphics.SmoothingMode;
|
||||
|
||||
using (Pen pen = new Pen(Color.White, 1.6f))
|
||||
{
|
||||
e.Graphics.SmoothingMode = SmoothingMode.AntiAlias;
|
||||
e.Graphics.DrawLine(pen, rect.Left, rect.Bottom - size / 2, rect.Left + size / 2.5f, rect.Bottom - 2);
|
||||
e.Graphics.DrawLine(pen, rect.Left + size / 2.6f, rect.Bottom - 2, rect.Right, rect.Top);
|
||||
}
|
||||
|
||||
e.Graphics.SmoothingMode = mode;
|
||||
}
|
||||
|
||||
protected override void OnRenderArrow(ToolStripArrowRenderEventArgs e)
|
||||
{
|
||||
if (!e.Item.Selected && !e.Item.Pressed)
|
||||
|
||||
@@ -32,6 +32,8 @@ namespace MeshExplorer.Controls
|
||||
Zoom zoom;
|
||||
RenderData data;
|
||||
bool initialized = false;
|
||||
VoronoiRenderer voronoi;
|
||||
bool showVoronoi = false;
|
||||
|
||||
string coordinate = String.Empty;
|
||||
|
||||
@@ -39,6 +41,21 @@ namespace MeshExplorer.Controls
|
||||
|
||||
public long RenderTime { get; private set; }
|
||||
public RenderData Data { get { return data; } }
|
||||
public bool ShowVoronoi
|
||||
{
|
||||
get { return showVoronoi; }
|
||||
set
|
||||
{
|
||||
showVoronoi = value;
|
||||
|
||||
if (voronoi != null && showVoronoi)
|
||||
{
|
||||
voronoi.Update();
|
||||
}
|
||||
|
||||
this.Render();
|
||||
}
|
||||
}
|
||||
|
||||
public MeshRenderer()
|
||||
{
|
||||
@@ -83,6 +100,9 @@ namespace MeshExplorer.Controls
|
||||
|
||||
public void SetData(Mesh mesh)
|
||||
{
|
||||
voronoi = new VoronoiRenderer(mesh);
|
||||
voronoi.Update();
|
||||
|
||||
data.SetData(mesh);
|
||||
|
||||
initialized = true;
|
||||
@@ -259,6 +279,11 @@ namespace MeshExplorer.Controls
|
||||
this.RenderTriangles(g);
|
||||
}
|
||||
|
||||
if (voronoi != null && this.showVoronoi)
|
||||
{
|
||||
voronoi.Render(g, zoom);
|
||||
}
|
||||
|
||||
if (data.Segments != null)
|
||||
{
|
||||
this.RenderSegments(g);
|
||||
|
||||
@@ -13,6 +13,7 @@ namespace MeshExplorer
|
||||
using TriangleNet;
|
||||
using TriangleNet.IO;
|
||||
using TriangleNet.Geometry;
|
||||
using MeshExplorer.IO;
|
||||
|
||||
/// <summary>
|
||||
/// Code of the online examples.
|
||||
|
||||
@@ -4,7 +4,7 @@
|
||||
// </copyright>
|
||||
// -----------------------------------------------------------------------
|
||||
|
||||
namespace MeshExplorer
|
||||
namespace MeshExplorer.IO
|
||||
{
|
||||
using System;
|
||||
using System.Drawing;
|
||||
@@ -14,6 +14,7 @@ namespace MeshExplorer
|
||||
using TriangleNet;
|
||||
using TriangleNet.Data;
|
||||
using TriangleNet.IO;
|
||||
using TriangleNet.Tools;
|
||||
|
||||
/// <summary>
|
||||
/// Writes an image of the mesh to disk.
|
||||
@@ -187,7 +188,7 @@ namespace MeshExplorer
|
||||
filename = String.Format("mesh-{0}.png", DateTime.Now.ToString("yyyy-M-d-hh-mm-ss"));
|
||||
}
|
||||
|
||||
VoronoiData data = DataWriter.WriteVoronoi(mesh);
|
||||
Voronoi data = null; // DataWriter.WriteVoronoi(mesh);
|
||||
|
||||
int n = data.Points.Length;
|
||||
|
||||
@@ -332,7 +333,7 @@ namespace MeshExplorer
|
||||
/// <summary>
|
||||
/// Draw mesh to the graphics object.
|
||||
/// </summary>
|
||||
private static void DrawVoronoi(Graphics g, Mesh mesh, VoronoiData voronoi, float scale)
|
||||
private static void DrawVoronoi(Graphics g, Mesh mesh, Voronoi voronoi, float scale)
|
||||
{
|
||||
g.SmoothingMode = SmoothingMode.AntiAlias;
|
||||
// Colors
|
||||
@@ -399,16 +400,16 @@ namespace MeshExplorer
|
||||
}
|
||||
|
||||
// Draw input points
|
||||
n = voronoi.InputPoints.Length;
|
||||
//n = voronoi.InputPoints.Length;
|
||||
|
||||
for (int i = 0; i < n; i++)
|
||||
{
|
||||
x = (float)voronoi.InputPoints[i].X;
|
||||
y = (float)voronoi.InputPoints[i].Y;
|
||||
//for (int i = 0; i < n; i++)
|
||||
//{
|
||||
// x = (float)voronoi.InputPoints[i].X;
|
||||
// y = (float)voronoi.InputPoints[i].Y;
|
||||
|
||||
g.FillEllipse(spBrush, x - radius, y - radius, 2 * radius, 2 * radius);
|
||||
//g.DrawEllipse(spBrush, x - radius, y - radius, 2 * radius, 2 * radius);
|
||||
}
|
||||
// g.FillEllipse(spBrush, x - radius, y - radius, 2 * radius, 2 * radius);
|
||||
// //g.DrawEllipse(spBrush, x - radius, y - radius, 2 * radius, 2 * radius);
|
||||
//}
|
||||
|
||||
bgBrush.Dispose();
|
||||
ptBrush.Dispose();
|
||||
@@ -418,12 +419,13 @@ namespace MeshExplorer
|
||||
trBrush.Dispose();
|
||||
}
|
||||
|
||||
private static PointF VoronoiBoxIntersection(BBox bounds, TriangleNet.Geometry.Point pt, double[] direction)
|
||||
private static PointF VoronoiBoxIntersection(BBox bounds, TriangleNet.Geometry.Point pt,
|
||||
TriangleNet.Geometry.Point direction)
|
||||
{
|
||||
double x = pt.X;
|
||||
double y = pt.Y;
|
||||
double dx = direction[0];
|
||||
double dy = direction[1];
|
||||
double dx = direction.X;
|
||||
double dy = direction.Y;
|
||||
|
||||
double t1, x1, y1, t2, x2, y2;
|
||||
|
||||
@@ -94,7 +94,7 @@
|
||||
<Compile Include="FormQuality.Designer.cs">
|
||||
<DependentUpon>FormQuality.cs</DependentUpon>
|
||||
</Compile>
|
||||
<Compile Include="ImageWriter.cs" />
|
||||
<Compile Include="IO\ImageWriter.cs" />
|
||||
<Compile Include="IO\FileProcessor.cs" />
|
||||
<Compile Include="IO\Formats\DatFile.cs" />
|
||||
<Compile Include="IO\Formats\JsonFile.cs" />
|
||||
@@ -105,6 +105,7 @@
|
||||
<Compile Include="Program.cs" />
|
||||
<Compile Include="Properties\AssemblyInfo.cs" />
|
||||
<Compile Include="Rendering\RenderData.cs" />
|
||||
<Compile Include="Rendering\VoronoiRenderer.cs" />
|
||||
<Compile Include="Rendering\Zoom.cs" />
|
||||
<Compile Include="Settings.cs" />
|
||||
<Compile Include="Util.cs" />
|
||||
|
||||
@@ -5,11 +5,11 @@ using System.Runtime.InteropServices;
|
||||
// General Information about an assembly is controlled through the following
|
||||
// set of attributes. Change these attribute values to modify the information
|
||||
// associated with an assembly.
|
||||
[assembly: AssemblyTitle("Example1")]
|
||||
[assembly: AssemblyTitle("Mesh Explorer")]
|
||||
[assembly: AssemblyDescription("")]
|
||||
[assembly: AssemblyConfiguration("")]
|
||||
[assembly: AssemblyCompany("")]
|
||||
[assembly: AssemblyProduct("Example1")]
|
||||
[assembly: AssemblyProduct("Mesh Explorer")]
|
||||
[assembly: AssemblyCopyright("Copyright © 2012")]
|
||||
[assembly: AssemblyTrademark("")]
|
||||
[assembly: AssemblyCulture("")]
|
||||
|
||||
@@ -0,0 +1,273 @@
|
||||
// -----------------------------------------------------------------------
|
||||
// <copyright file="VoronoiRenderer.cs" company="">
|
||||
// Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/triangle/
|
||||
// </copyright>
|
||||
// -----------------------------------------------------------------------
|
||||
|
||||
namespace MeshExplorer.Rendering
|
||||
{
|
||||
using System;
|
||||
using System.Collections.Generic;
|
||||
using System.Drawing;
|
||||
using System.Linq;
|
||||
using System.Text;
|
||||
using TriangleNet;
|
||||
using TriangleNet.Tools;
|
||||
|
||||
/// <summary>
|
||||
/// Renders a (bounded) Voronoi diagram.
|
||||
/// </summary>
|
||||
public class VoronoiRenderer
|
||||
{
|
||||
Mesh mesh;
|
||||
Voronoi simpleVoro;
|
||||
BoundedVoronoi boundedVoro;
|
||||
|
||||
Pen lines = new Pen(Color.FromArgb(40, 50, 60));
|
||||
|
||||
public VoronoiRenderer(Mesh mesh)
|
||||
{
|
||||
this.mesh = mesh;
|
||||
|
||||
//if (mesh.NumberOfSegments > 0)
|
||||
if (mesh.IsPolygon)
|
||||
{
|
||||
boundedVoro = new BoundedVoronoi(mesh);
|
||||
}
|
||||
else
|
||||
{
|
||||
simpleVoro = new Voronoi(mesh);
|
||||
}
|
||||
}
|
||||
|
||||
public void Update()
|
||||
{
|
||||
if (simpleVoro != null)
|
||||
{
|
||||
simpleVoro.Generate();
|
||||
}
|
||||
|
||||
if (boundedVoro != null)
|
||||
{
|
||||
boundedVoro.Generate();
|
||||
}
|
||||
}
|
||||
|
||||
public void Reset()
|
||||
{
|
||||
simpleVoro = null;
|
||||
boundedVoro = null;
|
||||
}
|
||||
|
||||
internal void Render(Graphics g, Zoom zoom)
|
||||
{
|
||||
if (simpleVoro != null)
|
||||
{
|
||||
RenderSimple(g, zoom);
|
||||
return;
|
||||
}
|
||||
|
||||
if (boundedVoro != null)
|
||||
{
|
||||
RenderBounded(g, zoom);
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
private void RenderSimple(Graphics g, Zoom zoom)
|
||||
{
|
||||
PointF p0, p1;
|
||||
|
||||
BBox bounds = new BBox(mesh.Bounds);
|
||||
|
||||
TriangleNet.Geometry.Point[] points = simpleVoro.Points;
|
||||
|
||||
// Enlarge 50%
|
||||
bounds.Extend(0.5f);
|
||||
|
||||
// Draw edges
|
||||
int n = simpleVoro.Edges == null ? 0 : simpleVoro.Edges.Length;
|
||||
|
||||
for (int i = 0; i < n; i++)
|
||||
{
|
||||
var seg = simpleVoro.Edges[i];
|
||||
|
||||
if (seg.P1 == -1)
|
||||
{
|
||||
// Infinite voronoi edge
|
||||
p0 = new PointF((float)points[seg.P0].X, (float)points[seg.P0].Y);
|
||||
|
||||
if (!BoxRayIntersection(bounds, points[seg.P0], simpleVoro.Directions[i], out p1))
|
||||
{
|
||||
continue;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
p0 = new PointF((float)points[seg.P0].X, (float)points[seg.P0].Y);
|
||||
p1 = new PointF((float)points[seg.P1].X, (float)points[seg.P1].Y);
|
||||
}
|
||||
|
||||
// Draw line
|
||||
RenderEdge(g, zoom, p0, p1);
|
||||
}
|
||||
|
||||
// Shrink 50%
|
||||
bounds.Extend(-0.5f);
|
||||
|
||||
// Scale the points radius to 2 pixel.
|
||||
//float radius = 1.5f / scale, x, y;
|
||||
|
||||
// Draw points
|
||||
//n = voro.Points.Length;
|
||||
|
||||
//for (int i = 0; i < n; i++)
|
||||
//{
|
||||
// x = (float)voro.Points[i].X;
|
||||
// y = (float)voro.Points[i].Y;
|
||||
|
||||
// g.FillEllipse(Brushes.Blue, x - radius, y - radius, 2 * radius, 2 * radius);
|
||||
//}
|
||||
}
|
||||
|
||||
private void RenderBounded(Graphics g, Zoom zoom)
|
||||
{
|
||||
PointF p0, p1;
|
||||
int n;
|
||||
|
||||
foreach (var cell in boundedVoro.Cells)
|
||||
{
|
||||
n = cell.Length;
|
||||
|
||||
for (int i = 1; i < n; i++)
|
||||
{
|
||||
p0 = new PointF((float)cell[i - 1].X, (float)cell[i - 1].Y);
|
||||
p1 = new PointF((float)cell[i].X, (float)cell[i].Y);
|
||||
|
||||
RenderEdge(g, zoom, p0, p1);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
private void RenderEdge(Graphics g, Zoom zoom, PointF p0, PointF p1)
|
||||
{
|
||||
if (zoom.ViewportContains(p0) ||
|
||||
zoom.ViewportContains(p1))
|
||||
{
|
||||
p0 = zoom.WorldToScreen(p0);
|
||||
p1 = zoom.WorldToScreen(p1);
|
||||
|
||||
g.DrawLine(lines, p0, p1);
|
||||
}
|
||||
}
|
||||
|
||||
private bool BoxRayIntersection(BBox bounds, TriangleNet.Geometry.Point pt,
|
||||
TriangleNet.Geometry.Point direction, out PointF intersect)
|
||||
{
|
||||
double x = pt.X;
|
||||
double y = pt.Y;
|
||||
double dx = direction.X;
|
||||
double dy = direction.Y;
|
||||
|
||||
double t1, x1, y1, t2, x2, y2;
|
||||
|
||||
intersect = new PointF();
|
||||
|
||||
// Check if point is inside the bounds
|
||||
if (x < bounds.MinX || x > bounds.MaxX || y < bounds.MinY || y > bounds.MaxY)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
// Calculate the cut through the vertical boundaries
|
||||
if (dx < 0)
|
||||
{
|
||||
// Line going to the left: intersect with x = bounds.MinX
|
||||
t1 = (bounds.MinX - x) / dx;
|
||||
x1 = bounds.MinX;
|
||||
y1 = y + t1 * dy;
|
||||
}
|
||||
else if (dx > 0)
|
||||
{
|
||||
// Line going to the right: intersect with x = bounds.MaxX
|
||||
t1 = (bounds.MaxX - x) / dx;
|
||||
x1 = bounds.MaxX;
|
||||
y1 = y + t1 * dy;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Line going straight up or down: no intersection possible
|
||||
t1 = double.MaxValue;
|
||||
x1 = y1 = 0;
|
||||
}
|
||||
|
||||
// Calculate the cut through upper and lower boundaries
|
||||
if (dy < 0)
|
||||
{
|
||||
// Line going downwards: intersect with y = bounds.MinY
|
||||
t2 = (bounds.MinY - y) / dy;
|
||||
x2 = x + t2 * dx;
|
||||
y2 = bounds.MinY;
|
||||
}
|
||||
else if (dx > 0)
|
||||
{
|
||||
// Line going upwards: intersect with y = bounds.MaxY
|
||||
t2 = (bounds.MaxY - y) / dy;
|
||||
x2 = x + t2 * dx;
|
||||
y2 = bounds.MaxY;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Horizontal line: no intersection possible
|
||||
t2 = double.MaxValue;
|
||||
x2 = y2 = 0;
|
||||
}
|
||||
|
||||
if (t1 < t2)
|
||||
{
|
||||
intersect.X = (float)x1;
|
||||
intersect.Y = (float)y1;
|
||||
}
|
||||
else
|
||||
{
|
||||
intersect.X = (float)x2;
|
||||
intersect.Y = (float)y2;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Bounding box.
|
||||
/// </summary>
|
||||
struct BBox
|
||||
{
|
||||
public float MinX;
|
||||
public float MaxX;
|
||||
public float MinY;
|
||||
public float MaxY;
|
||||
|
||||
public float Width { get { return MaxX - MinX; } }
|
||||
public float Height { get { return MaxY - MinY; } }
|
||||
|
||||
public BBox(TriangleNet.Geometry.BoundingBox box)
|
||||
{
|
||||
MinX = (float)box.Xmin;
|
||||
MaxX = (float)box.Xmax;
|
||||
MinY = (float)box.Ymin;
|
||||
MaxY = (float)box.Ymax;
|
||||
}
|
||||
|
||||
public void Extend(float amount)
|
||||
{
|
||||
float dx = amount * this.Width;
|
||||
float dy = amount * this.Height;
|
||||
|
||||
MinX -= dx;
|
||||
MaxX += dx;
|
||||
MinY -= dy;
|
||||
MaxY += dy;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -15,7 +15,7 @@ namespace MeshExplorer.Rendering
|
||||
/// <summary>
|
||||
/// Manages the current world to screen transformation
|
||||
/// </summary>
|
||||
class Zoom
|
||||
internal class Zoom
|
||||
{
|
||||
// The complete mesh
|
||||
Rectangle Screen { get; set; }
|
||||
|
||||
@@ -1,24 +0,0 @@
|
||||
// -----------------------------------------------------------------------
|
||||
// <copyright file="TriangulateIO.cs">
|
||||
// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html
|
||||
// Triangle.NET code by Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/index.html
|
||||
// </copyright>
|
||||
// -----------------------------------------------------------------------
|
||||
|
||||
namespace TriangleNet.IO
|
||||
{
|
||||
using TriangleNet.Geometry;
|
||||
|
||||
/// <summary>
|
||||
/// Stores the voronoi data (output only).
|
||||
/// </summary>
|
||||
public class VoronoiData
|
||||
{
|
||||
public Point[] InputPoints;
|
||||
public Point[] Points;
|
||||
public Edge[] Edges;
|
||||
|
||||
// Stores the direction for infinite voronoi edges
|
||||
public double[][] Directions;
|
||||
}
|
||||
}
|
||||
@@ -152,6 +152,11 @@ namespace TriangleNet
|
||||
/// </summary>
|
||||
public int NumberOfEdges { get { return this.edges; } }
|
||||
|
||||
/// <summary>
|
||||
/// Indicates whether the input is a PSLG or a point set.
|
||||
/// </summary>
|
||||
public bool IsPolygon { get { return this.insegments > 0; } }
|
||||
|
||||
#endregion
|
||||
|
||||
/// <summary>
|
||||
@@ -2040,7 +2045,7 @@ namespace TriangleNet
|
||||
/// triangles, but in the end every vertex will point to some triangle
|
||||
/// that contains it.
|
||||
/// </remarks>
|
||||
private void MakeVertexMap()
|
||||
internal void MakeVertexMap()
|
||||
{
|
||||
Otri tri = default(Otri);
|
||||
Vertex triorg;
|
||||
|
||||
@@ -142,7 +142,7 @@ namespace TriangleNet
|
||||
/// </summary>
|
||||
public bool CheckDelaunay()
|
||||
{
|
||||
Otri triangleloop = default(Otri);
|
||||
Otri loop = default(Otri);
|
||||
Otri oppotri = default(Otri);
|
||||
Osub opposubseg = default(Osub);
|
||||
Vertex triorg, tridest, triapex;
|
||||
@@ -157,18 +157,18 @@ namespace TriangleNet
|
||||
horrors = 0;
|
||||
|
||||
// Run through the list of triangles, checking each one.
|
||||
foreach (var t in mesh.triangles.Values)
|
||||
foreach (var tri in mesh.triangles.Values)
|
||||
{
|
||||
triangleloop.triangle = t;
|
||||
loop.triangle = tri;
|
||||
|
||||
// Check all three edges of the triangle.
|
||||
for (triangleloop.orient = 0; triangleloop.orient < 3;
|
||||
triangleloop.orient++)
|
||||
for (loop.orient = 0; loop.orient < 3;
|
||||
loop.orient++)
|
||||
{
|
||||
triorg = triangleloop.Org();
|
||||
tridest = triangleloop.Dest();
|
||||
triapex = triangleloop.Apex();
|
||||
triangleloop.Sym(ref oppotri);
|
||||
triorg = loop.Org();
|
||||
tridest = loop.Dest();
|
||||
triapex = loop.Apex();
|
||||
loop.Sym(ref oppotri);
|
||||
oppoapex = oppotri.Apex();
|
||||
// Only test that the edge is locally Delaunay if there is an
|
||||
// adjoining triangle whose pointer is larger (to ensure that
|
||||
@@ -187,7 +187,7 @@ namespace TriangleNet
|
||||
{
|
||||
// If a subsegment separates the triangles, then the edge is
|
||||
// constrained, so no local Delaunay test should be done.
|
||||
triangleloop.SegPivot(ref opposubseg);
|
||||
loop.SegPivot(ref opposubseg);
|
||||
if (opposubseg.seg != Mesh.dummysub)
|
||||
{
|
||||
shouldbedelaunay = false;
|
||||
|
||||
@@ -0,0 +1,551 @@
|
||||
// -----------------------------------------------------------------------
|
||||
// <copyright file="BoundedVoronoi.cs" company="">
|
||||
// Triangle.NET code by Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/index.html
|
||||
// </copyright>
|
||||
// -----------------------------------------------------------------------
|
||||
|
||||
namespace TriangleNet.Tools
|
||||
{
|
||||
using System;
|
||||
using System.Collections.Generic;
|
||||
using System.Linq;
|
||||
using System.Text;
|
||||
using TriangleNet.Data;
|
||||
using TriangleNet.Geometry;
|
||||
|
||||
/// <summary>
|
||||
/// The Bounded Voronoi Diagram is the dual of a PSLG triangulation.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// 2D Centroidal Voronoi Tessellations with Constraints, 2010,
|
||||
/// Jane Tournois, Pierre Alliez and Olivier Devillers
|
||||
/// </remarks>
|
||||
public class BoundedVoronoi
|
||||
{
|
||||
Mesh mesh;
|
||||
Dictionary<int, Segment> subsegMap;
|
||||
List<Point[]> voronoi;
|
||||
|
||||
/// <summary>
|
||||
/// Initializes a new instance of the <see cref="BoundedVoronoi" /> class.
|
||||
/// </summary>
|
||||
/// <param name="mesh">Mesh instance.</param>
|
||||
public BoundedVoronoi(Mesh mesh)
|
||||
{
|
||||
this.mesh = mesh;
|
||||
this.voronoi = new List<Point[]>();
|
||||
|
||||
// TODO: introduce isVertexMapValid in mesh and don't call
|
||||
// MakeVertexMap() if not necessary.
|
||||
this.mesh.MakeVertexMap();
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Gets the list of Voronoi cells.
|
||||
/// </summary>
|
||||
public List<Point[]> Cells
|
||||
{
|
||||
get { return voronoi; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Computes the bounded voronoi diagram.
|
||||
/// </summary>
|
||||
public void Generate()
|
||||
{
|
||||
voronoi.Clear();
|
||||
|
||||
TagBlindTriangles();
|
||||
|
||||
foreach (var v in mesh.vertices.Values)
|
||||
{
|
||||
if (v.type == VertexType.FreeVertex)
|
||||
{
|
||||
voronoi.Add(ConstructBvdCell(v));
|
||||
}
|
||||
else
|
||||
{
|
||||
voronoi.Add(ConstructBoundaryBvdCell(v));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
private void TagBlindTriangles()
|
||||
{
|
||||
int blinded = 0;
|
||||
|
||||
Stack<Triangle> triangles;
|
||||
subsegMap = new Dictionary<int, Segment>();
|
||||
|
||||
Otri f = default(Otri);
|
||||
Otri f0 = default(Otri);
|
||||
Osub e = default(Osub);
|
||||
Osub sub1 = default(Osub);
|
||||
|
||||
// Tag all triangles non-blind
|
||||
foreach (var t in mesh.triangles.Values)
|
||||
{
|
||||
// Use the infected flag for 'blinded' attribute.
|
||||
t.infected = false;
|
||||
}
|
||||
|
||||
// for each constrained edge e of cdt do
|
||||
foreach (var ss in mesh.subsegs.Values)
|
||||
{
|
||||
// Create a stack: triangles
|
||||
triangles = new Stack<Triangle>();
|
||||
|
||||
// for both adjacent triangles fe to e tagged non-blind do
|
||||
// Push fe into triangles
|
||||
e.seg = ss;
|
||||
e.orient = 0;
|
||||
e.TriPivot(ref f);
|
||||
|
||||
if (f.triangle != Mesh.dummytri && !f.triangle.infected)
|
||||
{
|
||||
triangles.Push(f.triangle);
|
||||
}
|
||||
|
||||
e.SymSelf();
|
||||
e.TriPivot(ref f);
|
||||
|
||||
if (f.triangle != Mesh.dummytri && !f.triangle.infected)
|
||||
{
|
||||
triangles.Push(f.triangle);
|
||||
}
|
||||
|
||||
// while triangles is non-empty
|
||||
while (triangles.Count > 0)
|
||||
{
|
||||
// Pop f from stack triangles
|
||||
f.triangle = triangles.Pop();
|
||||
f.orient = 0;
|
||||
|
||||
// if f is blinded by e (use P) then
|
||||
if (TriangleIsBlinded(ref f, ref e))
|
||||
{
|
||||
// Tag f as blinded by e
|
||||
f.triangle.infected = true;
|
||||
blinded++;
|
||||
|
||||
// Store association triangle -> subseg
|
||||
subsegMap.Add(f.triangle.hash, e.seg);
|
||||
|
||||
// for each adjacent triangle f0 to f do
|
||||
for (f.orient = 0; f.orient < 3; f.orient++)
|
||||
{
|
||||
f.Sym(ref f0);
|
||||
|
||||
f0.SegPivot(ref sub1);
|
||||
|
||||
// if f0 is finite and tagged non-blind & the common edge
|
||||
// between f and f0 is unconstrained then
|
||||
if (f0.triangle != Mesh.dummytri && !f0.triangle.infected && sub1.seg == Mesh.dummysub)
|
||||
{
|
||||
// Push f0 into triangles.
|
||||
triangles.Push(f0.triangle);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
blinded = 0;
|
||||
}
|
||||
|
||||
private bool TriangleIsBlinded(ref Otri tri, ref Osub seg)
|
||||
{
|
||||
Point cc, pt = new Point();
|
||||
double xi = 0, eta = 0;
|
||||
|
||||
Vertex torg = tri.Org();
|
||||
Vertex tdest = tri.Dest();
|
||||
Vertex tapex = tri.Apex();
|
||||
|
||||
cc = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false);
|
||||
|
||||
if (SegmentsIntersect(ref seg, cc, torg, ref pt, true))
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
if (SegmentsIntersect(ref seg, cc, tdest, ref pt, true))
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
if (SegmentsIntersect(ref seg, cc, tapex, ref pt, true))
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
private Point[] ConstructBvdCell(Vertex x)
|
||||
{
|
||||
Otri f = default(Otri);
|
||||
Otri f_init = default(Otri);
|
||||
Otri f_next = default(Otri);
|
||||
Osub sf = default(Osub);
|
||||
Osub sfn = default(Osub);
|
||||
|
||||
Vertex torg, tdest, tapex;
|
||||
Point cc_f, cc_f_next, p;
|
||||
double xi = 0, eta = 0;
|
||||
|
||||
// Call P the polygon (cell) in construction
|
||||
List<Point> P = new List<Point>();
|
||||
|
||||
// Call f_init a triangle incident to x
|
||||
x.tri.Copy(ref f_init);
|
||||
|
||||
if (f_init.Org() != x)
|
||||
{
|
||||
throw new Exception("ConstructBvdCell: inconsistent topology.");
|
||||
}
|
||||
|
||||
// Let f be initialized to f_init
|
||||
f_init.Copy(ref f);
|
||||
// Call f_next the next triangle counterclockwise around x
|
||||
f_init.Onext(ref f_next);
|
||||
|
||||
// repeat ... until f = f_init
|
||||
do
|
||||
{
|
||||
// Call Lffnext the line going through the circumcenters of f and f_next
|
||||
torg = f.Org();
|
||||
tdest = f.Dest();
|
||||
tapex = f.Apex();
|
||||
cc_f = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false);
|
||||
|
||||
torg = f_next.Org();
|
||||
tdest = f_next.Dest();
|
||||
tapex = f_next.Apex();
|
||||
cc_f_next = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false);
|
||||
|
||||
// if f is tagged non-blind then
|
||||
if (!f.triangle.infected)
|
||||
{
|
||||
// Insert the circumcenter of f into P
|
||||
P.Add(cc_f);
|
||||
|
||||
if (f_next.triangle.infected)
|
||||
{
|
||||
// Call S_fnext the constrained edge blinding f_next
|
||||
sfn.seg = subsegMap[f_next.triangle.hash];
|
||||
|
||||
p = new Point();
|
||||
// Insert point Lf,f_next /\ Sf_next into P
|
||||
if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Call Sf the constrained edge blinding f
|
||||
sf.seg = subsegMap[f.triangle.hash];
|
||||
|
||||
// if f_next is tagged non-blind then
|
||||
if (!f_next.triangle.infected)
|
||||
{
|
||||
p = new Point();
|
||||
// Insert point Lf,f_next /\ Sf into P
|
||||
if (SegmentsIntersect(ref sf, cc_f, cc_f_next, ref p, true))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Call Sf_next the constrained edge blinding f_next
|
||||
sfn.seg = subsegMap[f_next.triangle.hash];
|
||||
|
||||
// if Sf != Sf_next then
|
||||
if (!sf.Equal(sfn))
|
||||
{
|
||||
p = new Point();
|
||||
// Insert Lf,fnext /\ Sf and Lf,fnext /\ Sfnext into P
|
||||
if (SegmentsIntersect(ref sf, cc_f, cc_f_next, ref p, true))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
|
||||
p = new Point();
|
||||
if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// f <- f_next
|
||||
f_next.Copy(ref f);
|
||||
|
||||
// Call f_next the next triangle counterclockwise around x
|
||||
f_next.OnextSelf();
|
||||
}
|
||||
while (!f.Equal(f_init));
|
||||
|
||||
// Output: Bounded Voronoi cell of x in counterclockwise order.
|
||||
return P.ToArray();
|
||||
}
|
||||
|
||||
private Point[] ConstructBoundaryBvdCell(Vertex x)
|
||||
{
|
||||
Otri f = default(Otri);
|
||||
Otri f_init = default(Otri);
|
||||
Otri f_next = default(Otri);
|
||||
Otri f_prev = default(Otri);
|
||||
Osub sf = default(Osub);
|
||||
Osub sfn = default(Osub);
|
||||
|
||||
Vertex torg, tdest, tapex;
|
||||
Point cc_f, cc_f_next, p;
|
||||
double xi = 0, eta = 0;
|
||||
|
||||
// Call P the polygon (cell) in construction
|
||||
List<Point> P = new List<Point>();
|
||||
|
||||
// Call f_init a triangle incident to x
|
||||
x.tri.Copy(ref f_init);
|
||||
|
||||
if (f_init.Org() != x)
|
||||
{
|
||||
throw new Exception("ConstructBoundaryBvdCell: inconsistent topology.");
|
||||
}
|
||||
// Let f be initialized to f_init
|
||||
f_init.Copy(ref f);
|
||||
// Call f_next the next triangle counterclockwise around x
|
||||
f_init.Onext(ref f_next);
|
||||
|
||||
f_init.Oprev(ref f_prev);
|
||||
|
||||
// Is the border to the left?
|
||||
if (f_prev.triangle != Mesh.dummytri)
|
||||
{
|
||||
// Go clockwise until we reach the border (or the initial triangle)
|
||||
while (f_prev.triangle != Mesh.dummytri && !f_prev.Equal(f_init))
|
||||
{
|
||||
f_prev.Copy(ref f);
|
||||
f_prev.OprevSelf();
|
||||
}
|
||||
|
||||
f.Copy(ref f_init);
|
||||
f.Onext(ref f_next);
|
||||
}
|
||||
|
||||
// Add vertex on border
|
||||
P.Add(x);
|
||||
|
||||
// Add midpoint of start triangles' edge.
|
||||
torg = f.Org();
|
||||
tdest = f.Dest();
|
||||
P.Add(new Point((torg.X + tdest.X) / 2, (torg.Y + tdest.Y) / 2));
|
||||
|
||||
// repeat ... until f = f_init
|
||||
do
|
||||
{
|
||||
// Call Lffnext the line going through the circumcenters of f and f_next
|
||||
torg = f.Org();
|
||||
tdest = f.Dest();
|
||||
tapex = f.Apex();
|
||||
cc_f = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false);
|
||||
|
||||
if (f_next.triangle == Mesh.dummytri)
|
||||
{
|
||||
if (!f.triangle.infected)
|
||||
{
|
||||
// Add last circumcenter
|
||||
P.Add(cc_f);
|
||||
}
|
||||
|
||||
// Add midpoint of last triangles' edge (chances are it has already
|
||||
// been added, so post process cell to remove duplicates???)
|
||||
torg = f.Org();
|
||||
tapex = f.Apex();
|
||||
P.Add(new Point((torg.X + tapex.X) / 2, (torg.Y + tapex.Y) / 2));
|
||||
|
||||
break;
|
||||
}
|
||||
|
||||
torg = f_next.Org();
|
||||
tdest = f_next.Dest();
|
||||
tapex = f_next.Apex();
|
||||
cc_f_next = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false);
|
||||
|
||||
// if f is tagged non-blind then
|
||||
if (!f.triangle.infected)
|
||||
{
|
||||
// Insert the circumcenter of f into P
|
||||
P.Add(cc_f);
|
||||
|
||||
if (f_next.triangle.infected)
|
||||
{
|
||||
// Call S_fnext the constrained edge blinding f_next
|
||||
sfn.seg = subsegMap[f_next.triangle.hash];
|
||||
|
||||
p = new Point();
|
||||
// Insert point Lf,f_next /\ Sf_next into P
|
||||
if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Call Sf the constrained edge blinding f
|
||||
sf.seg = subsegMap[f.triangle.hash];
|
||||
|
||||
// if f_next is tagged non-blind then
|
||||
if (!f_next.triangle.infected)
|
||||
{
|
||||
tdest = f.Dest();
|
||||
tapex = f.Apex();
|
||||
|
||||
// Both circumcenters lie on the blinded side, but we
|
||||
// have to add the intersection with the segment.
|
||||
p = new Point();
|
||||
|
||||
// Center of f edge dest->apex
|
||||
Point bisec = new Point((tdest.X + tapex.X) / 2, (tdest.Y + tapex.Y) / 2);
|
||||
|
||||
// Find intersection of seg with line through f's bisector and circumcenter
|
||||
if (SegmentsIntersect(ref sf, bisec, cc_f, ref p, false))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
|
||||
|
||||
|
||||
p = new Point();
|
||||
// Insert point Lf,f_next /\ Sf into P
|
||||
if (SegmentsIntersect(ref sf, cc_f, cc_f_next, ref p, true))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Call Sf_next the constrained edge blinding f_next
|
||||
sfn.seg = subsegMap[f_next.triangle.hash];
|
||||
|
||||
// if Sf != Sf_next then
|
||||
if (!sf.Equal(sfn))
|
||||
{
|
||||
p = new Point();
|
||||
// Insert Lf,fnext /\ Sf and Lf,fnext /\ Sfnext into P
|
||||
if (SegmentsIntersect(ref sf, cc_f, cc_f_next, ref p, true))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
|
||||
p = new Point();
|
||||
if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Both circumcenters lie on the blinded side, but we
|
||||
// have to add the intersection with the segment.
|
||||
p = new Point();
|
||||
|
||||
// Center of f_next edge org->dest
|
||||
Point bisec = new Point((torg.X + tdest.X) / 2, (torg.Y + tdest.Y) / 2);
|
||||
|
||||
// Find intersection of seg with line through f_next's bisector and circumcenter
|
||||
if (SegmentsIntersect(ref sf, bisec, cc_f_next, ref p, false))
|
||||
{
|
||||
P.Add(p);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// f <- f_next
|
||||
f_next.Copy(ref f);
|
||||
|
||||
// Call f_next the next triangle counterclockwise around x
|
||||
f_next.OnextSelf();
|
||||
}
|
||||
while (!f.Equal(f_init));
|
||||
|
||||
// Output: Bounded Voronoi cell of x in counterclockwise order.
|
||||
return P.ToArray();
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Determines the intersection point of the line segment defined by points A and B with the
|
||||
/// line segment defined by points C and D.
|
||||
/// </summary>
|
||||
/// <param name="seg">The first segment AB.</param>
|
||||
/// <param name="pc">Endpoint C of second segment.</param>
|
||||
/// <param name="pd">Endpoint D of second segment.</param>
|
||||
/// <param name="p">Reference to the intersection point.</param>
|
||||
/// <param name="strictIntersect">If false, pa and pb represent a line.</param>
|
||||
/// <returns>Returns true if the intersection point was found, and stores that point in X,Y.
|
||||
/// Returns false if there is no determinable intersection point, in which case X,Y will
|
||||
/// be unmodified.
|
||||
/// </returns>
|
||||
private bool SegmentsIntersect(ref Osub seg, Point pc, Point pd, ref Point p, bool strictIntersect)
|
||||
{
|
||||
Vertex s1 = seg.Org();
|
||||
Vertex s2 = seg.Dest();
|
||||
|
||||
double Ax = s1.X, Ay = s1.Y;
|
||||
double Bx = s2.X, By = s2.Y;
|
||||
double Cx = pc.X, Cy = pc.Y;
|
||||
double Dx = pd.X, Dy = pd.Y;
|
||||
|
||||
double distAB, theCos, theSin, newX, ABpos;
|
||||
|
||||
// Fail if either line segment is zero-length.
|
||||
if (Ax == Bx && Ay == By || Cx == Dx && Cy == Dy) return false;
|
||||
|
||||
// Fail if the segments share an end-point.
|
||||
if (Ax == Cx && Ay == Cy || Bx == Cx && By == Cy
|
||||
|| Ax == Dx && Ay == Dy || Bx == Dx && By == Dy)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
// (1) Translate the system so that point A is on the origin.
|
||||
Bx -= Ax; By -= Ay;
|
||||
Cx -= Ax; Cy -= Ay;
|
||||
Dx -= Ax; Dy -= Ay;
|
||||
|
||||
// Discover the length of segment A-B.
|
||||
distAB = Math.Sqrt(Bx * Bx + By * By);
|
||||
|
||||
// (2) Rotate the system so that point B is on the positive X axis.
|
||||
theCos = Bx / distAB;
|
||||
theSin = By / distAB;
|
||||
newX = Cx * theCos + Cy * theSin;
|
||||
Cy = Cy * theCos - Cx * theSin; Cx = newX;
|
||||
newX = Dx * theCos + Dy * theSin;
|
||||
Dy = Dy * theCos - Dx * theSin; Dx = newX;
|
||||
|
||||
// Fail if segment C-D doesn't cross line A-B.
|
||||
if (Cy < 0 && Dy < 0 || Cy >= 0 && Dy >= 0 && strictIntersect) return false;
|
||||
|
||||
// (3) Discover the position of the intersection point along line A-B.
|
||||
ABpos = Dx + (Cx - Dx) * Dy / (Dy - Cy);
|
||||
|
||||
// Fail if segment C-D crosses line A-B outside of segment A-B.
|
||||
if (ABpos < 0 || ABpos > distAB && strictIntersect) return false;
|
||||
|
||||
// (4) Apply the discovered position to line A-B in the original coordinate system.
|
||||
p.x = Ax + ABpos * theCos;
|
||||
p.y = Ay + ABpos * theSin;
|
||||
|
||||
// Success.
|
||||
return true;
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,541 @@
|
||||
// -----------------------------------------------------------------------
|
||||
// <copyright file="QualityMeasure.cs" company="">
|
||||
// TODO: Update copyright text.
|
||||
// </copyright>
|
||||
// -----------------------------------------------------------------------
|
||||
|
||||
namespace TriangleNet.Tools
|
||||
{
|
||||
using System;
|
||||
using System.Collections.Generic;
|
||||
using System.Linq;
|
||||
using System.Text;
|
||||
using TriangleNet.Geometry;
|
||||
|
||||
/// <summary>
|
||||
/// Provides mesh quality information.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Given a triangle abc with points A (ax, ay), B (bx, by), C (cx, cy).
|
||||
///
|
||||
/// The side lengths are given as
|
||||
/// a = sqrt((cx - bx)^2 + (cy - by)^2) -- side BC opposite of A
|
||||
/// b = sqrt((cx - ax)^2 + (cy - ay)^2) -- side CA opposite of B
|
||||
/// c = sqrt((ax - bx)^2 + (ay - by)^2) -- side AB opposite of C
|
||||
///
|
||||
/// The angles are given as
|
||||
/// ang_a = acos((b^2 + c^2 - a^2) / (2 * b * c)) -- angle at A
|
||||
/// ang_b = acos((c^2 + a^2 - b^2) / (2 * c * a)) -- angle at B
|
||||
/// ang_c = acos((a^2 + b^2 - c^2) / (2 * a * b)) -- angle at C
|
||||
///
|
||||
/// The semiperimeter is given as
|
||||
/// s = (a + b + c) / 2
|
||||
///
|
||||
/// The area is given as
|
||||
/// D = abs(ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) / 2
|
||||
/// = sqrt(s * (s - a) * (s - b) * (s - c))
|
||||
///
|
||||
/// The inradius is given as
|
||||
/// r = D / s
|
||||
///
|
||||
/// The circumradius is given as
|
||||
/// R = a * b * c / (4 * D)
|
||||
///
|
||||
/// The altitudes are given as
|
||||
/// alt_a = 2 * D / a -- altitude above side a
|
||||
/// alt_b = 2 * D / b -- altitude above side b
|
||||
/// alt_c = 2 * D / c -- altitude above side c
|
||||
///
|
||||
/// The aspect ratio may be given as the ratio of the longest to the
|
||||
/// shortest edge or, more commonly as the ratio of the circumradius
|
||||
/// to twice the inradius
|
||||
/// ar = R / (2 * r)
|
||||
/// = a * b * c / (8 * (s - a) * (s - b) * (s - c))
|
||||
/// = a * b * c / ((b + c - a) * (c + a - b) * (a + b - c))
|
||||
/// </remarks>
|
||||
public class QualityMeasure
|
||||
{
|
||||
AreaMeasure areaMeasure;
|
||||
AlphaMeasure alphaMeasure;
|
||||
Q_Measure qMeasure;
|
||||
|
||||
Mesh mesh;
|
||||
|
||||
public QualityMeasure()
|
||||
{
|
||||
areaMeasure = new AreaMeasure();
|
||||
alphaMeasure = new AlphaMeasure();
|
||||
qMeasure = new Q_Measure();
|
||||
}
|
||||
|
||||
#region Public properties
|
||||
|
||||
/// <summary>
|
||||
/// Minimum triangle area.
|
||||
/// </summary>
|
||||
public double AreaMinimum
|
||||
{
|
||||
get { return areaMeasure.area_min; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Maximum triangle area.
|
||||
/// </summary>
|
||||
public double AreaMaximum
|
||||
{
|
||||
get { return areaMeasure.area_max; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Ratio of maximum and minimum triangle area.
|
||||
/// </summary>
|
||||
public double AreaRatio
|
||||
{
|
||||
get { return areaMeasure.area_max / areaMeasure.area_min; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Smallest angle.
|
||||
/// </summary>
|
||||
public double AlphaMinimum
|
||||
{
|
||||
get { return alphaMeasure.alpha_min; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Maximum smallest angle.
|
||||
/// </summary>
|
||||
public double AlphaMaximum
|
||||
{
|
||||
get { return alphaMeasure.alpha_max; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Average angle.
|
||||
/// </summary>
|
||||
public double AlphaAverage
|
||||
{
|
||||
get { return alphaMeasure.alpha_ave; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Average angle weighted by area.
|
||||
/// </summary>
|
||||
public double AlphaArea
|
||||
{
|
||||
get { return alphaMeasure.alpha_area; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Smallest aspect ratio.
|
||||
/// </summary>
|
||||
public double Q_Minimum
|
||||
{
|
||||
get { return qMeasure.q_min; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Largest aspect ratio.
|
||||
/// </summary>
|
||||
public double Q_Maximum
|
||||
{
|
||||
get { return qMeasure.q_max; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Average aspect ratio.
|
||||
/// </summary>
|
||||
public double Q_Average
|
||||
{
|
||||
get { return qMeasure.q_ave; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Average aspect ratio weighted by area.
|
||||
/// </summary>
|
||||
public double Q_Area
|
||||
{
|
||||
get { return qMeasure.q_area; }
|
||||
}
|
||||
|
||||
#endregion
|
||||
|
||||
public void Update(Mesh mesh)
|
||||
{
|
||||
this.mesh = mesh;
|
||||
|
||||
// Reset all measures.
|
||||
areaMeasure.Reset();
|
||||
alphaMeasure.Reset();
|
||||
qMeasure.Reset();
|
||||
|
||||
Compute();
|
||||
}
|
||||
|
||||
private void Compute()
|
||||
{
|
||||
Point a, b, c;
|
||||
double ab, bc, ca;
|
||||
double lx, ly;
|
||||
double area;
|
||||
|
||||
int n = 0;
|
||||
|
||||
foreach (var tri in mesh.triangles.Values)
|
||||
{
|
||||
n++;
|
||||
|
||||
a = tri.vertices[0];
|
||||
b = tri.vertices[1];
|
||||
c = tri.vertices[2];
|
||||
|
||||
lx = a.x - b.x;
|
||||
ly = a.y - b.y;
|
||||
ab = Math.Sqrt(lx * lx + ly * ly);
|
||||
lx = b.x - c.x;
|
||||
ly = b.y - c.y;
|
||||
bc = Math.Sqrt(lx * lx + ly * ly);
|
||||
lx = c.x - a.x;
|
||||
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);
|
||||
}
|
||||
|
||||
// Normalize measures
|
||||
alphaMeasure.Normalize(n, areaMeasure.area_total);
|
||||
qMeasure.Normalize(n, areaMeasure.area_total);
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Determines the bandwidth of the coefficient matrix.
|
||||
/// </summary>
|
||||
/// <returns>Bandwidth of the coefficient matrix.</returns>
|
||||
/// <remarks>
|
||||
/// The quantity computed here is the "geometric" bandwidth determined
|
||||
/// by the finite element mesh alone.
|
||||
///
|
||||
/// If a single finite element variable is associated with each node
|
||||
/// of the mesh, and if the nodes and variables are numbered in the
|
||||
/// same way, then the geometric bandwidth is the same as the bandwidth
|
||||
/// of a typical finite element matrix.
|
||||
///
|
||||
/// The bandwidth M is defined in terms of the lower and upper bandwidths:
|
||||
///
|
||||
/// M = ML + 1 + MU
|
||||
///
|
||||
/// where
|
||||
///
|
||||
/// ML = maximum distance from any diagonal entry to a nonzero
|
||||
/// entry in the same row, but earlier column,
|
||||
///
|
||||
/// MU = maximum distance from any diagonal entry to a nonzero
|
||||
/// entry in the same row, but later column.
|
||||
///
|
||||
/// Because the finite element node adjacency relationship is symmetric,
|
||||
/// we are guaranteed that ML = MU.
|
||||
/// </remarks>
|
||||
public int Bandwidth()
|
||||
{
|
||||
// Lower and upper bandwidth of the matrix
|
||||
int ml = 0, mu = 0;
|
||||
|
||||
int gi, gj;
|
||||
|
||||
foreach (var tri in mesh.triangles.Values)
|
||||
{
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
gi = tri[j];
|
||||
|
||||
for (int k = 0; k < 3; k++)
|
||||
{
|
||||
gj = tri[k];
|
||||
|
||||
mu = Math.Max(mu, gj - gi);
|
||||
ml = Math.Max(ml, gi - gj);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return ml + 1 + mu;
|
||||
}
|
||||
|
||||
class AreaMeasure
|
||||
{
|
||||
// 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;
|
||||
// Nmber of triangles with zero area
|
||||
public int area_zero = 0;
|
||||
|
||||
/// <summary>
|
||||
/// Reset all values.
|
||||
/// </summary>
|
||||
public void Reset()
|
||||
{
|
||||
area_min = double.MaxValue;
|
||||
area_max = -double.MaxValue;
|
||||
area_total = 0;
|
||||
area_zero = 0;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Compute the area of given triangle.
|
||||
/// </summary>
|
||||
/// <param name="a">Triangle corner a.</param>
|
||||
/// <param name="b">Triangle corner b.</param>
|
||||
/// <param name="c">Triangle corner c.</param>
|
||||
/// <returns>Triangle area.</returns>
|
||||
public double Measure(Point a, Point b, Point c)
|
||||
{
|
||||
double area = 0.5 * Math.Abs(a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y));
|
||||
|
||||
area_min = Math.Min(area_min, area);
|
||||
area_max = Math.Max(area_max, area);
|
||||
area_total += area;
|
||||
|
||||
if (area == 0.0)
|
||||
{
|
||||
area_zero = area_zero + 1;
|
||||
}
|
||||
|
||||
return area;
|
||||
}
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// The alpha measure determines the triangulated pointset quality.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// The alpha measure evaluates the uniformity of the shapes of the triangles
|
||||
/// defined by a triangulated pointset.
|
||||
///
|
||||
/// 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.
|
||||
/// </remarks>
|
||||
class AlphaMeasure
|
||||
{
|
||||
// Minimum value over all triangles
|
||||
public double alpha_min;
|
||||
// Maximum value over all triangles
|
||||
public double alpha_max;
|
||||
// Value averaged over all triangles
|
||||
public double alpha_ave;
|
||||
// Value averaged over all triangles and weighted by area
|
||||
public double alpha_area;
|
||||
|
||||
/// <summary>
|
||||
/// Reset all values.
|
||||
/// </summary>
|
||||
public void Reset()
|
||||
{
|
||||
alpha_min = double.MaxValue;
|
||||
alpha_max = -double.MaxValue;
|
||||
alpha_ave = 0;
|
||||
alpha_area = 0;
|
||||
}
|
||||
|
||||
double acos(double c)
|
||||
{
|
||||
if (c <= -1.0)
|
||||
{
|
||||
return Math.PI;
|
||||
}
|
||||
else if (1.0 <= c)
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
else
|
||||
{
|
||||
return Math.Acos(c);
|
||||
}
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Compute q value of given triangle.
|
||||
/// </summary>
|
||||
/// <param name="ab">Side length ab.</param>
|
||||
/// <param name="bc">Side length bc.</param>
|
||||
/// <param name="ca">Side length ca.</param>
|
||||
/// <param name="area">Triangle area.</param>
|
||||
/// <returns></returns>
|
||||
public double Measure(double ab, double bc, double ca, double area)
|
||||
{
|
||||
double alpha = double.MaxValue;
|
||||
|
||||
double ab2 = ab * ab;
|
||||
double bc2 = bc * bc;
|
||||
double ca2 = ca * ca;
|
||||
|
||||
double a_angle;
|
||||
double b_angle;
|
||||
double c_angle;
|
||||
|
||||
// Take care of a ridiculous special case.
|
||||
if (ab == 0.0 && bc == 0.0 && ca == 0.0)
|
||||
{
|
||||
a_angle = 2.0 * Math.PI / 3.0;
|
||||
b_angle = 2.0 * Math.PI / 3.0;
|
||||
c_angle = 2.0 * Math.PI / 3.0;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (ca == 0.0 || ab == 0.0)
|
||||
{
|
||||
a_angle = Math.PI;
|
||||
}
|
||||
else
|
||||
{
|
||||
a_angle = acos((ca2 + ab2 - bc2) / (2.0 * ca * ab));
|
||||
}
|
||||
|
||||
if (ab == 0.0 || bc == 0.0)
|
||||
{
|
||||
b_angle = Math.PI;
|
||||
}
|
||||
else
|
||||
{
|
||||
b_angle = acos((ab2 + bc2 - ca2) / (2.0 * ab * bc));
|
||||
}
|
||||
|
||||
if (bc == 0.0 || ca == 0.0)
|
||||
{
|
||||
c_angle = Math.PI;
|
||||
}
|
||||
else
|
||||
{
|
||||
c_angle = acos((bc2 + ca2 - ab2) / (2.0 * bc * ca));
|
||||
}
|
||||
}
|
||||
|
||||
alpha = Math.Min(alpha, a_angle);
|
||||
alpha = Math.Min(alpha, b_angle);
|
||||
alpha = Math.Min(alpha, c_angle);
|
||||
|
||||
// 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;
|
||||
|
||||
alpha_min = Math.Min(alpha, alpha_min);
|
||||
alpha_max = Math.Max(alpha, alpha_max);
|
||||
|
||||
return alpha;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Normalize values.
|
||||
/// </summary>
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// The Q measure determines the triangulated pointset quality.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// The Q measure evaluates the uniformity of the shapes of the triangles
|
||||
/// defined by a triangulated pointset. It uses the aspect ratio
|
||||
///
|
||||
/// 2 * (incircle radius) / (circumcircle radius)
|
||||
///
|
||||
/// In an ideally regular mesh, all triangles would have the same
|
||||
/// equilateral shape, for which Q = 1. A good mesh would have
|
||||
/// 0.5 < Q.
|
||||
/// </remarks>
|
||||
class Q_Measure
|
||||
{
|
||||
// Minimum value over all triangles
|
||||
public double q_min;
|
||||
// Maximum value over all triangles
|
||||
public double q_max;
|
||||
// Average value
|
||||
public double q_ave;
|
||||
// Average value weighted by the area of each triangle
|
||||
public double q_area;
|
||||
|
||||
/// <summary>
|
||||
/// Reset all values.
|
||||
/// </summary>
|
||||
public void Reset()
|
||||
{
|
||||
q_min = double.MaxValue;
|
||||
q_max = -double.MaxValue;
|
||||
q_ave = 0;
|
||||
q_area = 0;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Compute q value of given triangle.
|
||||
/// </summary>
|
||||
/// <param name="ab">Side length ab.</param>
|
||||
/// <param name="bc">Side length bc.</param>
|
||||
/// <param name="ca">Side length ca.</param>
|
||||
/// <param name="area">Triangle area.</param>
|
||||
/// <returns></returns>
|
||||
public double Measure(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);
|
||||
|
||||
q_ave += q;
|
||||
q_area += q * area;
|
||||
|
||||
return q;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Normalize values.
|
||||
/// </summary>
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -1,25 +1,60 @@
|
||||
// -----------------------------------------------------------------------
|
||||
// <copyright file="io.cs" company="">
|
||||
// <copyright file="Voronoi.cs">
|
||||
// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html
|
||||
// Triangle.NET code by Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/triangle/
|
||||
// Triangle.NET code by Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/index.html
|
||||
// </copyright>
|
||||
// -----------------------------------------------------------------------
|
||||
|
||||
namespace TriangleNet.IO
|
||||
{
|
||||
using System;
|
||||
using System.IO;
|
||||
using System.Collections.Generic;
|
||||
using System.Globalization;
|
||||
using TriangleNet.Data;
|
||||
using TriangleNet.Geometry;
|
||||
using TriangleNet.Data;
|
||||
using TriangleNet.Geometry;
|
||||
|
||||
namespace TriangleNet.Tools
|
||||
{
|
||||
/// <summary>
|
||||
/// Generates a mesh representaion using arrays.
|
||||
/// The Voronoi Diagram is the dual of a pointset triangulation.
|
||||
/// </summary>
|
||||
public static class DataWriter
|
||||
public class Voronoi
|
||||
{
|
||||
#region Library
|
||||
Mesh mesh;
|
||||
Point[] points;
|
||||
Point[] directions; // Stores the direction for infinite voronoi edges
|
||||
Edge[] edges;
|
||||
|
||||
/// <summary>
|
||||
/// Initializes a new instance of the <see cref="BoundedVoronoi" /> class.
|
||||
/// </summary>
|
||||
/// <param name="mesh"></param>
|
||||
/// <remarks>
|
||||
/// Be sure MakeVertexMap has been called (should always be the case).
|
||||
/// </remarks>
|
||||
public Voronoi(Mesh mesh)
|
||||
{
|
||||
this.mesh = mesh;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Gets the list of Voronoi vertices.
|
||||
/// </summary>
|
||||
public Point[] Points
|
||||
{
|
||||
get { return points; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Gets the directions for infinite Voronoi edges.
|
||||
/// </summary>
|
||||
public Point[] Directions
|
||||
{
|
||||
get { return directions; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Gets the list of Voronoi edges.
|
||||
/// </summary>
|
||||
public Edge[] Edges
|
||||
{
|
||||
get { return edges; }
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Gets the Voronoi diagram as raw output data.
|
||||
@@ -31,42 +66,21 @@ namespace TriangleNet.IO
|
||||
/// Hence, the Voronoi vertices are listed by traversing the Delaunay
|
||||
/// triangles, and the Voronoi edges are listed by traversing the Delaunay
|
||||
/// edges.
|
||||
///
|
||||
/// WARNING: In order to assign numbers to the Voronoi vertices, this
|
||||
/// procedure messes up the subsegments or the extra nodes of every
|
||||
/// element. Hence, you should call this procedure last.</remarks>
|
||||
public static VoronoiData WriteVoronoi(Mesh mesh)
|
||||
///</remarks>
|
||||
public void Generate()
|
||||
{
|
||||
VoronoiData data = new VoronoiData();
|
||||
|
||||
Otri tri = default(Otri), trisym = default(Otri);
|
||||
Vertex torg, tdest, tapex;
|
||||
Point circumcenter;
|
||||
double xi = 0, eta = 0;
|
||||
int p1, p2;
|
||||
|
||||
int i = 0;
|
||||
|
||||
// Copy input points (actually not part of the voronoi diagram)
|
||||
data.InputPoints = new Vertex[mesh.vertices.Count];
|
||||
|
||||
foreach (var item in mesh.vertices.Values)
|
||||
{
|
||||
if (item.type != VertexType.UndeadVertex)
|
||||
{
|
||||
data.InputPoints[i] = item;
|
||||
i++;
|
||||
}
|
||||
}
|
||||
|
||||
// Allocate memory for Voronoi vertices.
|
||||
data.Points = new Vertex[mesh.triangles.Count];
|
||||
|
||||
int index = 0;
|
||||
this.points = new Point[mesh.triangles.Count];
|
||||
|
||||
tri.orient = 0;
|
||||
|
||||
i = 0;
|
||||
int i = 0;
|
||||
foreach (var item in mesh.triangles.Values)
|
||||
{
|
||||
tri.triangle = item;
|
||||
@@ -76,19 +90,19 @@ namespace TriangleNet.IO
|
||||
circumcenter = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false);
|
||||
|
||||
// X and y coordinates.
|
||||
data.Points[i] = new Point(circumcenter.x, circumcenter.y);
|
||||
this.points[i] = new Point(circumcenter.x, circumcenter.y);
|
||||
|
||||
// Update element id
|
||||
tri.triangle.id = i++;
|
||||
}
|
||||
|
||||
// Allocate memory for output Voronoi edges.
|
||||
data.Edges = new Edge[mesh.edges];
|
||||
this.edges = new Edge[mesh.edges];
|
||||
|
||||
// Allocate memory for output Voronoi norms.
|
||||
data.Directions = new double[mesh.edges][];
|
||||
this.directions = new Point[mesh.edges];
|
||||
|
||||
index = 0;
|
||||
i = 0;
|
||||
// To loop over the set of edges, loop over all triangles, and look at
|
||||
// the three edges of each triangle. If there isn't another triangle
|
||||
// adjacent to the edge, operate on the edge. If there is another
|
||||
@@ -113,8 +127,8 @@ namespace TriangleNet.IO
|
||||
tdest = tri.Dest();
|
||||
|
||||
// Copy an infinite ray. Index of one endpoint, and -1.
|
||||
data.Edges[index] = new Edge(p1, -1);
|
||||
data.Directions[index] = new double[] { tdest[1] - torg[1], torg[0] - tdest[0] };
|
||||
this.edges[i] = new Edge(p1, -1);
|
||||
this.directions[i] = new Point(tdest.y - torg.y, torg.x - tdest.x);
|
||||
}
|
||||
else
|
||||
{
|
||||
@@ -122,18 +136,14 @@ namespace TriangleNet.IO
|
||||
p2 = trisym.triangle.id;
|
||||
// Finite edge. Write indices of two endpoints.
|
||||
|
||||
data.Edges[index] = new Edge(p1, p2);
|
||||
data.Directions[index] = new double[] { 0, 0 };
|
||||
this.edges[i] = new Edge(p1, p2);
|
||||
this.directions[i] = new Point(0, 0);
|
||||
}
|
||||
|
||||
index++;
|
||||
i++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return data;
|
||||
}
|
||||
|
||||
#endregion
|
||||
}
|
||||
}
|
||||
@@ -64,8 +64,6 @@
|
||||
<Compile Include="IO\DebugWriter.cs" />
|
||||
<Compile Include="IO\FileWriter.cs" />
|
||||
<Compile Include="IO\InputTriangle.cs" />
|
||||
<Compile Include="IO\VoronoiData.cs" />
|
||||
<Compile Include="IO\DataWriter.cs" />
|
||||
<Compile Include="IO\FileReader.cs" />
|
||||
<Compile Include="Log\ILog.cs" />
|
||||
<Compile Include="Log\ILogItem.cs" />
|
||||
@@ -80,8 +78,11 @@
|
||||
<Compile Include="Properties\AssemblyInfo.cs" />
|
||||
<Compile Include="Sampler.cs" />
|
||||
<Compile Include="Smoothing\ISmoother.cs" />
|
||||
<Compile Include="Tools\BoundedVoronoi.cs" />
|
||||
<Compile Include="Tools\QualityMeasure.cs" />
|
||||
<Compile Include="Tools\Statistic.cs" />
|
||||
<Compile Include="Algorithm\SweepLine.cs" />
|
||||
<Compile Include="Tools\Voronoi.cs" />
|
||||
</ItemGroup>
|
||||
<ItemGroup />
|
||||
<Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
|
||||
|
||||
Reference in New Issue
Block a user