Initial commit

This commit is contained in:
AbstractConcept 2022-09-13 00:36:34 -05:00
commit 3c7cc0c973
8391 changed files with 704313 additions and 0 deletions

View file

@ -0,0 +1,286 @@
// -----------------------------------------------------------------------
// <copyright file="AdjacencyMatrix.cs" company="">
// Original Matlab code by John Burkardt, Florida State University
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace UnityEngine.U2D.Animation.TriangleNet
.Tools
{
using System;
/// <summary>
/// The adjacency matrix of the mesh.
/// </summary>
internal class AdjacencyMatrix
{
// Number of adjacency entries.
int nnz;
// Pointers into the actual adjacency structure adj. Information about row k is
// stored in entries pcol(k) through pcol(k+1)-1 of adj. Size: N + 1
int[] pcol;
// The adjacency structure. For each row, it contains the column indices
// of the nonzero entries. Size: nnz
int[] irow;
/// <summary>
/// Gets the number of columns (nodes of the mesh).
/// </summary>
public readonly int N;
/// <summary>
/// Gets the column pointers.
/// </summary>
public int[] ColumnPointers
{
get { return pcol; }
}
/// <summary>
/// Gets the row indices.
/// </summary>
public int[] RowIndices
{
get { return irow; }
}
public AdjacencyMatrix(Mesh mesh)
{
this.N = mesh.vertices.Count;
// Set up the adj_row adjacency pointer array.
this.pcol = AdjacencyCount(mesh);
this.nnz = pcol[N];
// Set up the adj adjacency array.
this.irow = AdjacencySet(mesh, this.pcol);
SortIndices();
}
public AdjacencyMatrix(int[] pcol, int[] irow)
{
this.N = pcol.Length - 1;
this.nnz = pcol[N];
this.pcol = pcol;
this.irow = irow;
if (pcol[0] != 0)
{
throw new ArgumentException("Expected 0-based indexing.", "pcol");
}
if (irow.Length < nnz)
{
throw new ArgumentException();
}
}
/// <summary>
/// Computes the bandwidth of an adjacency matrix.
/// </summary>
/// <returns>Bandwidth of the adjacency matrix.</returns>
public int Bandwidth()
{
int band_hi;
int band_lo;
int col;
int i, j;
band_lo = 0;
band_hi = 0;
for (i = 0; i < N; i++)
{
for (j = pcol[i]; j < pcol[i + 1]; j++)
{
col = irow[j];
band_lo = Math.Max(band_lo, i - col);
band_hi = Math.Max(band_hi, col - i);
}
}
return band_lo + 1 + band_hi;
}
#region Adjacency matrix
/// <summary>
/// Counts adjacencies in a triangulation.
/// </summary>
/// <remarks>
/// This routine is called to count the adjacencies, so that the
/// appropriate amount of memory can be set aside for storage when
/// the adjacency structure is created.
///
/// The triangulation is assumed to involve 3-node triangles.
///
/// Two nodes are "adjacent" if they are both nodes in some triangle.
/// Also, a node is considered to be adjacent to itself.
/// </remarks>
int[] AdjacencyCount(Mesh mesh)
{
int n = N;
int n1, n2, n3;
int tid, nid;
int[] pcol = new int[n + 1];
// Set every node to be adjacent to itself.
for (int i = 0; i < n; i++)
{
pcol[i] = 1;
}
// Examine each triangle.
foreach (var tri in mesh.triangles)
{
tid = tri.id;
n1 = tri.vertices[0].id;
n2 = tri.vertices[1].id;
n3 = tri.vertices[2].id;
// Add edge (1,2) if this is the first occurrence, that is, if
// the edge (1,2) is on a boundary (nid <= 0) or if this triangle
// is the first of the pair in which the edge occurs (tid < nid).
nid = tri.neighbors[2].tri.id;
if (nid < 0 || tid < nid)
{
pcol[n1] += 1;
pcol[n2] += 1;
}
// Add edge (2,3).
nid = tri.neighbors[0].tri.id;
if (nid < 0 || tid < nid)
{
pcol[n2] += 1;
pcol[n3] += 1;
}
// Add edge (3,1).
nid = tri.neighbors[1].tri.id;
if (nid < 0 || tid < nid)
{
pcol[n3] += 1;
pcol[n1] += 1;
}
}
// We used PCOL to count the number of entries in each column.
// Convert it to pointers into the ADJ array.
for (int i = n; i > 0; i--)
{
pcol[i] = pcol[i - 1];
}
pcol[0] = 0;
for (int i = 1; i <= n; i++)
{
pcol[i] = pcol[i - 1] + pcol[i];
}
return pcol;
}
/// <summary>
/// Sets adjacencies in a triangulation.
/// </summary>
/// <remarks>
/// This routine can be used to create the compressed column storage
/// for a linear triangle finite element discretization of Poisson's
/// equation in two dimensions.
/// </remarks>
int[] AdjacencySet(Mesh mesh, int[] pcol)
{
int n = this.N;
int[] col = new int[n];
// Copy of the adjacency rows input.
Array.Copy(pcol, col, n);
int i, nnz = pcol[n];
// Output list, stores the actual adjacency information.
int[] list = new int[nnz];
// Set every node to be adjacent to itself.
for (i = 0; i < n; i++)
{
list[col[i]] = i;
col[i] += 1;
}
int n1, n2, n3; // Vertex numbers.
int tid, nid; // Triangle and neighbor id.
// Examine each triangle.
foreach (var tri in mesh.triangles)
{
tid = tri.id;
n1 = tri.vertices[0].id;
n2 = tri.vertices[1].id;
n3 = tri.vertices[2].id;
// Add edge (1,2) if this is the first occurrence, that is, if
// the edge (1,2) is on a boundary (nid <= 0) or if this triangle
// is the first of the pair in which the edge occurs (tid < nid).
nid = tri.neighbors[2].tri.id;
if (nid < 0 || tid < nid)
{
list[col[n1]++] = n2;
list[col[n2]++] = n1;
}
// Add edge (2,3).
nid = tri.neighbors[0].tri.id;
if (nid < 0 || tid < nid)
{
list[col[n2]++] = n3;
list[col[n3]++] = n2;
}
// Add edge (3,1).
nid = tri.neighbors[1].tri.id;
if (nid < 0 || tid < nid)
{
list[col[n1]++] = n3;
list[col[n3]++] = n1;
}
}
return list;
}
public void SortIndices()
{
int k1, k2, n = N;
int[] list = this.irow;
// Ascending sort the entries for each column.
for (int i = 0; i < n; i++)
{
k1 = pcol[i];
k2 = pcol[i + 1];
Array.Sort(list, k1, k2 - k1);
}
}
#endregion
}
}

View file

@ -0,0 +1,11 @@
fileFormatVersion: 2
guid: 826e5298bd8f8498a9a2f3a04a3b2fa9
MonoImporter:
externalObjects: {}
serializedVersion: 2
defaultReferences: []
executionOrder: 0
icon: {instanceID: 0}
userData:
assetBundleName:
assetBundleVariant:

View file

@ -0,0 +1,676 @@
// -----------------------------------------------------------------------
// <copyright file="CuthillMcKee.cs" company="">
// Original Matlab code by John Burkardt, Florida State University
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace UnityEngine.U2D.Animation.TriangleNet
.Tools
{
using System;
/// <summary>
/// Applies the Cuthill and McKee renumbering algorithm to reduce the bandwidth of
/// the adjacency matrix associated with the mesh.
/// </summary>
internal class CuthillMcKee
{
// The adjacency matrix of the mesh.
AdjacencyMatrix matrix;
/// <summary>
/// Gets the permutation vector for the Reverse Cuthill-McKee numbering.
/// </summary>
/// <param name="mesh">The mesh.</param>
/// <returns>Permutation vector.</returns>
public int[] Renumber(Mesh mesh)
{
// Algorithm needs linear numbering of the nodes.
mesh.Renumber(NodeNumbering.Linear);
return Renumber(new AdjacencyMatrix(mesh));
}
/// <summary>
/// Gets the permutation vector for the Reverse Cuthill-McKee numbering.
/// </summary>
/// <param name="mesh">The mesh.</param>
/// <returns>Permutation vector.</returns>
public int[] Renumber(AdjacencyMatrix matrix)
{
this.matrix = matrix;
int bandwidth1 = matrix.Bandwidth();
var pcol = matrix.ColumnPointers;
// Adjust column pointers (1-based indexing).
Shift(pcol, true);
// TODO: Make RCM work with 0-based matrix.
// Compute the RCM permutation.
int[] perm = GenerateRcm();
int[] perm_inv = PermInverse(perm);
int bandwidth2 = PermBandwidth(perm, perm_inv);
if (Log.Verbose)
{
Log.Instance.Info(String.Format("Reverse Cuthill-McKee (Bandwidth: {0} > {1})",
bandwidth1, bandwidth2));
}
// Adjust column pointers (0-based indexing).
Shift(pcol, false);
return perm_inv;
}
#region RCM
/// <summary>
/// Finds the reverse Cuthill-Mckee ordering for a general graph.
/// </summary>
/// <returns>The RCM ordering.</returns>
/// <remarks>
/// For each connected component in the graph, the routine obtains
/// an ordering by calling RCM.
/// </remarks>
int[] GenerateRcm()
{
// Number of nodes in the mesh.
int n = matrix.N;
int[] perm = new int[n];
int i, num, root;
int iccsze = 0;
int level_num = 0;
/// Index vector for a level structure. The level structure is stored in the
/// currently unused spaces in the permutation vector PERM.
int[] level_row = new int[n + 1];
/// Marks variables that have been numbered.
int[] mask = new int[n];
for (i = 0; i < n; i++)
{
mask[i] = 1;
}
num = 1;
for (i = 0; i < n; i++)
{
// For each masked connected component...
if (mask[i] != 0)
{
root = i;
// Find a pseudo-peripheral node ROOT. The level structure found by
// ROOT_FIND is stored starting at PERM(NUM).
FindRoot(ref root, mask, ref level_num, level_row, perm, num - 1);
// RCM orders the component using ROOT as the starting node.
Rcm(root, mask, perm, num - 1, ref iccsze);
num += iccsze;
// We can stop once every node is in one of the connected components.
if (n < num)
{
return perm;
}
}
}
return perm;
}
/// <summary>
/// RCM renumbers a connected component by the reverse Cuthill McKee algorithm.
/// </summary>
/// <param name="root">the node that defines the connected component. It is used as the starting
/// point for the RCM ordering.</param>
/// <param name="mask">Input/output, int MASK(NODE_NUM), a mask for the nodes. Only those nodes with
/// nonzero input mask values are considered by the routine. The nodes numbered by RCM will have
/// their mask values set to zero.</param>
/// <param name="perm">Output, int PERM(NODE_NUM), the RCM ordering.</param>
/// <param name="iccsze">Output, int ICCSZE, the size of the connected component that has been numbered.</param>
/// <param name="node_num">the number of nodes.</param>
/// <remarks>
/// The connected component is specified by a node ROOT and a mask.
/// The numbering starts at the root node.
///
/// An outline of the algorithm is as follows:
///
/// X(1) = ROOT.
///
/// for ( I = 1 to N-1)
/// Find all unlabeled neighbors of X(I),
/// assign them the next available labels, in order of increasing degree.
///
/// When done, reverse the ordering.
/// </remarks>
void Rcm(int root, int[] mask, int[] perm, int offset, ref int iccsze)
{
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int fnbr;
int i, j, k, l;
int jstop, jstrt;
int lbegin, lnbr, lperm, lvlend;
int nbr, node;
// Number of nodes in the mesh.
int n = matrix.N;
/// Workspace, int DEG[NODE_NUM], a temporary vector used to hold
/// the degree of the nodes in the section graph specified by mask and root.
int[] deg = new int[n];
// Find the degrees of the nodes in the component specified by MASK and ROOT.
Degree(root, mask, deg, ref iccsze, perm, offset);
mask[root] = 0;
if (iccsze <= 1)
{
return;
}
lvlend = 0;
lnbr = 1;
// LBEGIN and LVLEND point to the beginning and
// the end of the current level respectively.
while (lvlend < lnbr)
{
lbegin = lvlend + 1;
lvlend = lnbr;
for (i = lbegin; i <= lvlend; i++)
{
// For each node in the current level...
node = perm[offset + i - 1];
jstrt = pcol[node];
jstop = pcol[node + 1] - 1;
// Find the unnumbered neighbors of NODE.
// FNBR and LNBR point to the first and last neighbors
// of the current node in PERM.
fnbr = lnbr + 1;
for (j = jstrt; j <= jstop; j++)
{
nbr = irow[j - 1];
if (mask[nbr] != 0)
{
lnbr += 1;
mask[nbr] = 0;
perm[offset + lnbr - 1] = nbr;
}
}
// Node has neighbors
if (lnbr > fnbr)
{
// Sort the neighbors of NODE in increasing order by degree.
// Linear insertion is used.
k = fnbr;
while (k < lnbr)
{
l = k;
k = k + 1;
nbr = perm[offset + k - 1];
while (fnbr < l)
{
lperm = perm[offset + l - 1];
if (deg[lperm - 1] <= deg[nbr - 1])
{
break;
}
perm[offset + l] = lperm;
l = l - 1;
}
perm[offset + l] = nbr;
}
}
}
}
// We now have the Cuthill-McKee ordering. Reverse it.
ReverseVector(perm, offset, iccsze);
}
/// <summary>
/// Finds a pseudo-peripheral node.
/// </summary>
/// <param name="root">On input, ROOT is a node in the the component of the graph for
/// which a pseudo-peripheral node is sought. On output, ROOT is the pseudo-peripheral
/// node obtained.</param>
/// <param name="mask">MASK[NODE_NUM], specifies a section subgraph. Nodes for which MASK
/// is zero are ignored by FNROOT.</param>
/// <param name="level_num">Output, int LEVEL_NUM, is the number of levels in the level
/// structure rooted at the node ROOT.</param>
/// <param name="level_row">Output, int LEVEL_ROW(NODE_NUM+1), the level structure array pair
/// containing the level structure found.</param>
/// <param name="level">Output, int LEVEL(NODE_NUM), the level structure array pair
/// containing the level structure found.</param>
/// <param name="node_num">the number of nodes.</param>
/// <remarks>
/// The diameter of a graph is the maximum distance (number of edges)
/// between any two nodes of the graph.
///
/// The eccentricity of a node is the maximum distance between that
/// node and any other node of the graph.
///
/// A peripheral node is a node whose eccentricity equals the
/// diameter of the graph.
///
/// A pseudo-peripheral node is an approximation to a peripheral node;
/// it may be a peripheral node, but all we know is that we tried our
/// best.
///
/// The routine is given a graph, and seeks pseudo-peripheral nodes,
/// using a modified version of the scheme of Gibbs, Poole and
/// Stockmeyer. It determines such a node for the section subgraph
/// specified by MASK and ROOT.
///
/// The routine also determines the level structure associated with
/// the given pseudo-peripheral node; that is, how far each node
/// is from the pseudo-peripheral node. The level structure is
/// returned as a list of nodes LS, and pointers to the beginning
/// of the list of nodes that are at a distance of 0, 1, 2, ...,
/// NODE_NUM-1 from the pseudo-peripheral node.
///
/// Reference:
/// Alan George, Joseph Liu,
/// Computer Solution of Large Sparse Positive Definite Systems,
/// Prentice Hall, 1981.
///
/// Norman Gibbs, William Poole, Paul Stockmeyer,
/// An Algorithm for Reducing the Bandwidth and Profile of a Sparse Matrix,
/// SIAM Journal on Numerical Analysis,
/// Volume 13, pages 236-250, 1976.
///
/// Norman Gibbs,
/// Algorithm 509: A Hybrid Profile Reduction Algorithm,
/// ACM Transactions on Mathematical Software,
/// Volume 2, pages 378-387, 1976.
/// </remarks>
void FindRoot(ref int root, int[] mask, ref int level_num, int[] level_row,
int[] level, int offset)
{
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int iccsze;
int j, jstrt;
int k, kstop, kstrt;
int mindeg;
int nghbor, ndeg;
int node;
int level_num2 = 0;
// Determine the level structure rooted at ROOT.
GetLevelSet(ref root, mask, ref level_num, level_row, level, offset);
// Count the number of nodes in this level structure.
iccsze = level_row[level_num] - 1;
// Extreme cases:
// A complete graph has a level set of only a single level.
// Every node is equally good (or bad).
// or
// A "line graph" 0--0--0--0--0 has every node in its only level.
// By chance, we've stumbled on the ideal root.
if (level_num == 1 || level_num == iccsze)
{
return;
}
// Pick any node from the last level that has minimum degree
// as the starting point to generate a new level set.
for (;;)
{
mindeg = iccsze;
jstrt = level_row[level_num - 1];
root = level[offset + jstrt - 1];
if (jstrt < iccsze)
{
for (j = jstrt; j <= iccsze; j++)
{
node = level[offset + j - 1];
ndeg = 0;
kstrt = pcol[node - 1];
kstop = pcol[node] - 1;
for (k = kstrt; k <= kstop; k++)
{
nghbor = irow[k - 1];
if (mask[nghbor] > 0)
{
ndeg += 1;
}
}
if (ndeg < mindeg)
{
root = node;
mindeg = ndeg;
}
}
}
// Generate the rooted level structure associated with this node.
GetLevelSet(ref root, mask, ref level_num2, level_row, level, offset);
// If the number of levels did not increase, accept the new ROOT.
if (level_num2 <= level_num)
{
break;
}
level_num = level_num2;
// In the unlikely case that ROOT is one endpoint of a line graph,
// we can exit now.
if (iccsze <= level_num)
{
break;
}
}
}
/// <summary>
/// Generates the connected level structure rooted at a given node.
/// </summary>
/// <param name="root">the node at which the level structure is to be rooted.</param>
/// <param name="mask">MASK[NODE_NUM]. On input, only nodes with nonzero MASK are to be processed.
/// On output, those nodes which were included in the level set have MASK set to 1.</param>
/// <param name="level_num">Output, int LEVEL_NUM, the number of levels in the level structure. ROOT is
/// in level 1. The neighbors of ROOT are in level 2, and so on.</param>
/// <param name="level_row">Output, int LEVEL_ROW[NODE_NUM+1], the rooted level structure.</param>
/// <param name="level">Output, int LEVEL[NODE_NUM], the rooted level structure.</param>
/// <param name="node_num">the number of nodes.</param>
/// <remarks>
/// Only nodes for which MASK is nonzero will be considered.
///
/// The root node chosen by the user is assigned level 1, and masked.
/// All (unmasked) nodes reachable from a node in level 1 are
/// assigned level 2 and masked. The process continues until there
/// are no unmasked nodes adjacent to any node in the current level.
/// The number of levels may vary between 2 and NODE_NUM.
///
/// Reference:
/// Alan George, Joseph Liu,
/// Computer Solution of Large Sparse Positive Definite Systems,
/// Prentice Hall, 1981.
/// </remarks>
void GetLevelSet(ref int root, int[] mask, ref int level_num, int[] level_row,
int[] level, int offset)
{
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int i, iccsze;
int j, jstop, jstrt;
int lbegin, lvlend, lvsize;
int nbr;
int node;
mask[root] = 0;
level[offset] = root;
level_num = 0;
lvlend = 0;
iccsze = 1;
// LBEGIN is the pointer to the beginning of the current level, and
// LVLEND points to the end of this level.
for (;;)
{
lbegin = lvlend + 1;
lvlend = iccsze;
level_num += 1;
level_row[level_num - 1] = lbegin;
// Generate the next level by finding all the masked neighbors of nodes
// in the current level.
for (i = lbegin; i <= lvlend; i++)
{
node = level[offset + i - 1];
jstrt = pcol[node];
jstop = pcol[node + 1] - 1;
for (j = jstrt; j <= jstop; j++)
{
nbr = irow[j - 1];
if (mask[nbr] != 0)
{
iccsze += 1;
level[offset + iccsze - 1] = nbr;
mask[nbr] = 0;
}
}
}
// Compute the current level width (the number of nodes encountered.)
// If it is positive, generate the next level.
lvsize = iccsze - lvlend;
if (lvsize <= 0)
{
break;
}
}
level_row[level_num] = lvlend + 1;
// Reset MASK to 1 for the nodes in the level structure.
for (i = 0; i < iccsze; i++)
{
mask[level[offset + i]] = 1;
}
}
/// <summary>
/// Computes the degrees of the nodes in the connected component.
/// </summary>
/// <param name="root">the node that defines the connected component.</param>
/// <param name="mask">MASK[NODE_NUM], is nonzero for those nodes which are to be considered.</param>
/// <param name="deg">Output, int DEG[NODE_NUM], contains, for each node in the connected component, its degree.</param>
/// <param name="iccsze">Output, int ICCSIZE, the number of nodes in the connected component.</param>
/// <param name="ls">Output, int LS[NODE_NUM], stores in entries 1 through ICCSIZE the nodes in the
/// connected component, starting with ROOT, and proceeding by levels.</param>
/// <param name="node_num">the number of nodes.</param>
/// <remarks>
/// The connected component is specified by MASK and ROOT.
/// Nodes for which MASK is zero are ignored.
///
/// Reference:
/// Alan George, Joseph Liu,
/// Computer Solution of Large Sparse Positive Definite Systems,
/// Prentice Hall, 1981.
/// </remarks>
void Degree(int root, int[] mask, int[] deg, ref int iccsze, int[] ls, int offset)
{
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int i, ideg;
int j, jstop, jstrt;
int lbegin, lvlend;
int lvsize = 1;
int nbr, node;
// The sign of ADJ_ROW(I) is used to indicate if node I has been considered.
ls[offset] = root;
pcol[root] = -pcol[root];
lvlend = 0;
iccsze = 1;
// If the current level width is nonzero, generate another level.
while (lvsize > 0)
{
// LBEGIN is the pointer to the beginning of the current level, and
// LVLEND points to the end of this level.
lbegin = lvlend + 1;
lvlend = iccsze;
// Find the degrees of nodes in the current level,
// and at the same time, generate the next level.
for (i = lbegin; i <= lvlend; i++)
{
node = ls[offset + i - 1];
jstrt = -pcol[node];
jstop = Math.Abs(pcol[node + 1]) - 1;
ideg = 0;
for (j = jstrt; j <= jstop; j++)
{
nbr = irow[j - 1];
if (mask[nbr] != 0) // EDIT: [nbr - 1]
{
ideg = ideg + 1;
if (0 <= pcol[nbr]) // EDIT: [nbr - 1]
{
pcol[nbr] = -pcol[nbr]; // EDIT: [nbr - 1]
iccsze = iccsze + 1;
ls[offset + iccsze - 1] = nbr;
}
}
}
deg[node] = ideg;
}
// Compute the current level width.
lvsize = iccsze - lvlend;
}
// Reset ADJ_ROW to its correct sign and return.
for (i = 0; i < iccsze; i++)
{
node = ls[offset + i];
pcol[node] = -pcol[node];
}
}
#endregion
#region Tools
/// <summary>
/// Computes the bandwidth of a permuted adjacency matrix.
/// </summary>
/// <param name="perm">The permutation.</param>
/// <param name="perm_inv">The inverse permutation.</param>
/// <returns>Bandwidth of the permuted adjacency matrix.</returns>
/// <remarks>
/// The matrix is defined by the adjacency information and a permutation.
/// The routine also computes the bandwidth and the size of the envelope.
/// </remarks>
int PermBandwidth(int[] perm, int[] perm_inv)
{
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int col, i, j;
int band_lo = 0;
int band_hi = 0;
int n = matrix.N;
for (i = 0; i < n; i++)
{
for (j = pcol[perm[i]]; j < pcol[perm[i] + 1]; j++)
{
col = perm_inv[irow[j - 1]];
band_lo = Math.Max(band_lo, i - col);
band_hi = Math.Max(band_hi, col - i);
}
}
return band_lo + 1 + band_hi;
}
/// <summary>
/// Produces the inverse of a given permutation.
/// </summary>
/// <param name="n">Number of items permuted.</param>
/// <param name="perm">PERM[N], a permutation.</param>
/// <returns>The inverse permutation.</returns>
int[] PermInverse(int[] perm)
{
int n = matrix.N;
int[] perm_inv = new int[n];
for (int i = 0; i < n; i++)
{
perm_inv[perm[i]] = i;
}
return perm_inv;
}
/// <summary>
/// Reverses the elements of an integer vector.
/// </summary>
/// <param name="size">number of entries in the array.</param>
/// <param name="a">the array to be reversed.</param>
/// <example>
/// Input:
/// N = 5,
/// A = ( 11, 12, 13, 14, 15 ).
///
/// Output:
/// A = ( 15, 14, 13, 12, 11 ).
/// </example>
void ReverseVector(int[] a, int offset, int size)
{
int i;
int j;
for (i = 0; i < size / 2; i++)
{
j = a[offset + i];
a[offset + i] = a[offset + size - 1 - i];
a[offset + size - 1 - i] = j;
}
}
void Shift(int[] a, bool up)
{
int length = a.Length;
if (up)
{
for (int i = 0; i < length; a[i]++, i++) ;
}
else
{
for (int i = 0; i < length; a[i]--, i++) ;
}
}
#endregion
}
}

View file

@ -0,0 +1,11 @@
fileFormatVersion: 2
guid: 0115a1541423b477a94d872384667321
MonoImporter:
externalObjects: {}
serializedVersion: 2
defaultReferences: []
executionOrder: 0
icon: {instanceID: 0}
userData:
assetBundleName:
assetBundleVariant:

View file

@ -0,0 +1,107 @@
namespace UnityEngine.U2D.Animation.TriangleNet
.Tools
{
using Animation.TriangleNet.Geometry;
internal static class Interpolation
{
#if USE_ATTRIBS
/// <summary>
/// Linear interpolation of vertex attributes.
/// </summary>
/// <param name="vertex">The interpolation vertex.</param>
/// <param name="triangle">The triangle containing the vertex.</param>
/// <param name="n">The number of vertex attributes.</param>
/// <remarks>
/// The vertex is expected to lie inside the triangle.
/// </remarks>
internal static void InterpolateAttributes(Vertex vertex, ITriangle triangle, int n)
{
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;
// 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;
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;
// 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);
for (int i = 0; i < n; i++)
{
// Interpolate the vertex attributes.
vertex.attributes[i] = org.attributes[i]
+ xi * (dest.attributes[i] - org.attributes[i])
+ eta * (apex.attributes[i] - org.attributes[i]);
}
}
#endif
#if USE_Z
/// <summary>
/// Linear interpolation of a scalar value.
/// </summary>
/// <param name="p">The interpolation point.</param>
/// <param name="triangle">The triangle containing the point.</param>
/// <remarks>
/// The point is expected to lie inside the triangle.
/// </remarks>
internal 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;
// 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;
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;
// 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);
p.z = org.z + xi * (dest.z - org.z) + eta * (apex.z - org.z);
}
#endif
}
}

View file

@ -0,0 +1,11 @@
fileFormatVersion: 2
guid: 09c24a0e213174105acf5329eadc89bb
MonoImporter:
externalObjects: {}
serializedVersion: 2
defaultReferences: []
executionOrder: 0
icon: {instanceID: 0}
userData:
assetBundleName:
assetBundleVariant:

View file

@ -0,0 +1,226 @@
// -----------------------------------------------------------------------
// <copyright file="IntersectionHelper.cs" company="">
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace UnityEngine.U2D.Animation.TriangleNet
.Tools
{
using Animation.TriangleNet.Geometry;
internal static class IntersectionHelper
{
/// <summary>
/// Compute intersection of two segments.
/// </summary>
/// <param name="p0">Segment 1 start point.</param>
/// <param name="p1">Segment 1 end point.</param>
/// <param name="q0">Segment 2 start point.</param>
/// <param name="q1">Segment 2 end point.</param>
/// <param name="c0">The intersection point.</param>
/// <remarks>
/// This is a special case of segment intersection. Since the calling algorithm assures
/// that a valid intersection exists, there's no need to check for any special cases.
/// </remarks>
internal static void IntersectSegments(Point p0, Point p1, Point q0, Point q1, ref Point c0)
{
double ux = p1.x - p0.x;
double uy = p1.y - p0.y;
double vx = q1.x - q0.x;
double vy = q1.y - q0.y;
double wx = p0.x - q0.x;
double wy = p0.y - q0.y;
double d = (ux * vy - uy * vx);
double s = (vx * wy - vy * wx) / d;
// Intersection point
c0.x = p0.X + s * ux;
c0.y = p0.Y + s * uy;
}
/// <summary>
/// Intersect segment with a bounding box.
/// </summary>
/// <param name="rect">The clip rectangle.</param>
/// <param name="p0">Segment endpoint.</param>
/// <param name="p1">Segment endpoint.</param>
/// <param name="c0">The new location of p0.</param>
/// <param name="c1">The new location of p1.</param>
/// <returns>Returns true, if segment is clipped.</returns>
/// <remarks>
/// Based on Liang-Barsky function by Daniel White:
/// http://www.skytopia.com/project/articles/compsci/clipping.html
/// </remarks>
internal static bool LiangBarsky(Rectangle rect, Point p0, Point p1, ref Point c0, ref Point c1)
{
// Define the x/y clipping values for the border.
double xmin = rect.Left;
double xmax = rect.Right;
double ymin = rect.Bottom;
double ymax = rect.Top;
// Define the start and end points of the line.
double x0 = p0.X;
double y0 = p0.Y;
double x1 = p1.X;
double y1 = p1.Y;
double t0 = 0.0;
double t1 = 1.0;
double dx = x1 - x0;
double dy = y1 - y0;
double p = 0.0, q = 0.0, r;
for (int edge = 0; edge < 4; edge++)
{
// Traverse through left, right, bottom, top edges.
if (edge == 0) { p = -dx; q = -(xmin - x0); }
if (edge == 1) { p = dx; q = (xmax - x0); }
if (edge == 2) { p = -dy; q = -(ymin - y0); }
if (edge == 3) { p = dy; q = (ymax - y0); }
r = q / p;
if (p == 0 && q < 0) return false; // Don't draw line at all. (parallel line outside)
if (p < 0)
{
if (r > t1) return false; // Don't draw line at all.
else if (r > t0) t0 = r; // Line is clipped!
}
else if (p > 0)
{
if (r < t0) return false; // Don't draw line at all.
else if (r < t1) t1 = r; // Line is clipped!
}
}
c0.X = x0 + t0 * dx;
c0.Y = y0 + t0 * dy;
c1.X = x0 + t1 * dx;
c1.Y = y0 + t1 * dy;
return true; // (clipped) line is drawn
}
/// <summary>
/// Intersect a ray with a bounding box.
/// </summary>
/// <param name="rect">The clip rectangle.</param>
/// <param name="p0">The ray startpoint (inside the box).</param>
/// <param name="p1">Any point in ray direction (NOT the direction vector).</param>
/// <param name="c1">The intersection point.</param>
/// <returns>Returns false, if startpoint is outside the box.</returns>
internal static bool BoxRayIntersection(Rectangle rect, Point p0, Point p1, ref Point c1)
{
return BoxRayIntersection(rect, p0, p1.x - p0.x, p1.y - p0.y, ref c1);
}
/// <summary>
/// Intersect a ray with a bounding box.
/// </summary>
/// <param name="rect">The clip rectangle.</param>
/// <param name="p">The ray startpoint (inside the box).</param>
/// <param name="dx">X direction.</param>
/// <param name="dy">Y direction.</param>
/// <returns>Returns false, if startpoint is outside the box.</returns>
internal static Point BoxRayIntersection(Rectangle rect, Point p, double dx, double dy)
{
var intersection = new Point();
if (BoxRayIntersection(rect, p, dx, dy, ref intersection))
{
return intersection;
}
return null;
}
/// <summary>
/// Intersect a ray with a bounding box.
/// </summary>
/// <param name="rect">The clip rectangle.</param>
/// <param name="p">The ray startpoint (inside the box).</param>
/// <param name="dx">X direction.</param>
/// <param name="dy">Y direction.</param>
/// <param name="c">The intersection point.</param>
/// <returns>Returns false, if startpoint is outside the box.</returns>
internal static bool BoxRayIntersection(Rectangle rect, Point p, double dx, double dy, ref Point c)
{
double x = p.X;
double y = p.Y;
double t1, x1, y1, t2, x2, y2;
// Bounding box
double xmin = rect.Left;
double xmax = rect.Right;
double ymin = rect.Bottom;
double ymax = rect.Top;
// Check if point is inside the bounds
if (x < xmin || x > xmax || y < ymin || y > ymax)
{
return false;
}
// Calculate the cut through the vertical boundaries
if (dx < 0)
{
// Line going to the left: intersect with x = minX
t1 = (xmin - x) / dx;
x1 = xmin;
y1 = y + t1 * dy;
}
else if (dx > 0)
{
// Line going to the right: intersect with x = maxX
t1 = (xmax - x) / dx;
x1 = xmax;
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 = minY
t2 = (ymin - y) / dy;
x2 = x + t2 * dx;
y2 = ymin;
}
else if (dy > 0)
{
// Line going upwards: intersect with y = maxY
t2 = (ymax - y) / dy;
x2 = x + t2 * dx;
y2 = ymax;
}
else
{
// Horizontal line: no intersection possible
t2 = double.MaxValue;
x2 = y2 = 0;
}
if (t1 < t2)
{
c.x = x1;
c.y = y1;
}
else
{
c.x = x2;
c.y = y2;
}
return true;
}
}
}

View file

@ -0,0 +1,11 @@
fileFormatVersion: 2
guid: 5427b508fe0244b63ac33dac16710613
MonoImporter:
externalObjects: {}
serializedVersion: 2
defaultReferences: []
executionOrder: 0
icon: {instanceID: 0}
userData:
assetBundleName:
assetBundleVariant:

View file

@ -0,0 +1,246 @@
// -----------------------------------------------------------------------
// <copyright file="PolygonValidator.cs">
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace UnityEngine.U2D.Animation.TriangleNet
.Tools
{
using System;
using System.Collections.Generic;
using Animation.TriangleNet.Geometry;
internal static class PolygonValidator
{
/// <summary>
/// Test the polygon for consistency.
/// </summary>
internal static bool IsConsistent(IPolygon poly)
{
var logger = Log.Instance;
var points = poly.Points;
int horrors = 0;
int i = 0;
int count = points.Count;
if (count < 3)
{
logger.Warning("Polygon must have at least 3 vertices.", "PolygonValidator.IsConsistent()");
return false;
}
foreach (var p in points)
{
if (p == null)
{
horrors++;
logger.Warning(String.Format("Point {0} is null.", i), "PolygonValidator.IsConsistent()");
}
else if (double.IsNaN(p.x) || double.IsNaN(p.y))
{
horrors++;
logger.Warning(String.Format("Point {0} has invalid coordinates.", i), "PolygonValidator.IsConsistent()");
}
else if (double.IsInfinity(p.x) || double.IsInfinity(p.y))
{
horrors++;
logger.Warning(String.Format("Point {0} has invalid coordinates.", i), "PolygonValidator.IsConsistent()");
}
i++;
}
i = 0;
foreach (var seg in poly.Segments)
{
if (seg == null)
{
horrors++;
logger.Warning(String.Format("Segment {0} is null.", i), "PolygonValidator.IsConsistent()");
// Always abort if a NULL-segment is found.
return false;
}
var p = seg.GetVertex(0);
var q = seg.GetVertex(1);
if ((p.x == q.x) && (p.y == q.y))
{
horrors++;
logger.Warning(String.Format("Endpoints of segment {0} are coincident (IDs {1} / {2}).", i, p.id, q.id),
"PolygonValidator.IsConsistent()");
}
i++;
}
if (points[0].id == points[1].id)
{
horrors += CheckVertexIDs(poly, count);
}
else
{
horrors += CheckDuplicateIDs(poly);
}
return horrors == 0;
}
/// <summary>
/// Test the polygon for duplicate vertices.
/// </summary>
internal static bool HasDuplicateVertices(IPolygon poly)
{
var logger = Log.Instance;
int horrors = 0;
var points = poly.Points.ToArray();
VertexSorter.Sort(points);
for (int i = 1; i < points.Length; i++)
{
if (points[i - 1] == points[i])
{
horrors++;
logger.Warning(String.Format("Found duplicate point {0}.", points[i]),
"PolygonValidator.HasDuplicateVertices()");
}
}
return horrors > 0;
}
/// <summary>
/// Test the polygon for 360 degree angles.
/// </summary>
/// <param name="poly">The polygon.</param>
/// <param name="threshold">The angle threshold.</param>
internal static bool HasBadAngles(IPolygon poly, double threshold = 2e-12)
{
var logger = Log.Instance;
int horrors = 0;
int i = 0;
Point p0 = null, p1 = null;
Point q0, q1;
int count = poly.Points.Count;
foreach (var seg in poly.Segments)
{
q0 = p0; // Previous segment start point.
q1 = p1; // Previous segment end point.
p0 = seg.GetVertex(0); // Current segment start point.
p1 = seg.GetVertex(1); // Current segment end point.
if (p0 == p1 || q0 == q1)
{
// Ignore zero-length segments.
continue;
}
if (q0 != null && q1 != null)
{
// The two segments are connected.
if (p0 == q1 && p1 != null)
{
if (IsBadAngle(q0, p0, p1, threshold))
{
horrors++;
logger.Warning(String.Format("Bad segment angle found at index {0}.", i),
"PolygonValidator.HasBadAngles()");
}
}
}
i++;
}
return horrors > 0;
}
private static bool IsBadAngle(Point a, Point b, Point c, double threshold = 0.0)
{
double x = DotProduct(a, b, c);
double y = CrossProductLength(a, b, c);
return Math.Abs(Math.Atan2(y, x)) <= threshold;
}
// Returns the dot product <AB, BC>.
private static double DotProduct(Point a, Point b, Point c)
{
// Calculate the dot product.
return (a.x - b.x) * (c.x - b.x) + (a.y - b.y) * (c.y - b.y);
}
// Returns the length of cross product AB x BC.
private static double CrossProductLength(Point a, Point b, Point c)
{
// Calculate the Z coordinate of the cross product.
return (a.x - b.x) * (c.y - b.y) - (a.y - b.y) * (c.x - b.x);
}
private static int CheckVertexIDs(IPolygon poly, int count)
{
var logger = Log.Instance;
int horrors = 0;
int i = 0;
Vertex p, q;
foreach (var seg in poly.Segments)
{
p = seg.GetVertex(0);
q = seg.GetVertex(1);
if (p.id < 0 || p.id >= count)
{
horrors++;
logger.Warning(String.Format("Segment {0} has invalid startpoint.", i),
"PolygonValidator.IsConsistent()");
}
if (q.id < 0 || q.id >= count)
{
horrors++;
logger.Warning(String.Format("Segment {0} has invalid endpoint.", i),
"PolygonValidator.IsConsistent()");
}
i++;
}
return horrors;
}
private static int CheckDuplicateIDs(IPolygon poly)
{
var ids = new HashSet<int>();
// Check for duplicate ids.
foreach (var p in poly.Points)
{
if (!ids.Add(p.id))
{
Log.Instance.Warning("Found duplicate vertex ids.", "PolygonValidator.IsConsistent()");
return 1;
}
}
return 0;
}
}
}

View file

@ -0,0 +1,11 @@
fileFormatVersion: 2
guid: e7ad5906e3c5748eaa55d2843a8c129e
MonoImporter:
externalObjects: {}
serializedVersion: 2
defaultReferences: []
executionOrder: 0
icon: {instanceID: 0}
userData:
assetBundleName:
assetBundleVariant:

View file

@ -0,0 +1,542 @@
// -----------------------------------------------------------------------
// <copyright file="QualityMeasure.cs" company="">
// Original Matlab code by John Burkardt, Florida State University
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace UnityEngine.U2D.Animation.TriangleNet
.Tools
{
using System;
using Animation.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>
internal 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)
{
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 label to a nonzero
/// label in the same row, but earlier column,
///
/// MU = maximum distance from any diagonal label to a nonzero
/// label 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()
{
if (mesh == null) return 0;
// Lower and upper bandwidth of the matrix
int ml = 0, mu = 0;
int gi, gj;
foreach (var tri in mesh.triangles)
{
for (int j = 0; j < 3; j++)
{
gi = tri.GetVertex(j).id;
for (int k = 0; k < 3; k++)
{
gj = tri.GetVertex(k).id;
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 &lt; 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;
}
}
}
}
}

View file

@ -0,0 +1,11 @@
fileFormatVersion: 2
guid: 187a6ce7f29024eab88bfc87e4b5312c
MonoImporter:
externalObjects: {}
serializedVersion: 2
defaultReferences: []
executionOrder: 0
icon: {instanceID: 0}
userData:
assetBundleName:
assetBundleVariant:

View file

@ -0,0 +1,539 @@
// -----------------------------------------------------------------------
// <copyright file="Statistic.cs">
// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace UnityEngine.U2D.Animation.TriangleNet
.Tools
{
using System;
using Animation.TriangleNet.Topology;
using Animation.TriangleNet.Geometry;
/// <summary>
/// Gather mesh statistics.
/// </summary>
internal class Statistic
{
#region Static members
/// <summary>
/// Number of incircle tests performed.
/// </summary>
internal static long InCircleCount = 0;
internal static long InCircleAdaptCount = 0;
/// <summary>
/// Number of counterclockwise tests performed.
/// </summary>
internal static long CounterClockwiseCount = 0;
internal static long CounterClockwiseAdaptCount = 0;
/// <summary>
/// Number of 3D orientation tests performed.
/// </summary>
internal static long Orient3dCount = 0;
/// <summary>
/// Number of right-of-hyperbola tests performed.
/// </summary>
internal static long HyperbolaCount = 0;
/// <summary>
/// // Number of circumcenter calculations performed.
/// </summary>
internal static long CircumcenterCount = 0;
/// <summary>
/// Number of circle top calculations performed.
/// </summary>
internal static long CircleTopCount = 0;
/// <summary>
/// Number of vertex relocations.
/// </summary>
internal static long RelocationCount = 0;
#endregion
#region Properties
double minEdge = 0;
/// <summary>
/// Gets the shortest edge.
/// </summary>
public double ShortestEdge { get { return minEdge; } }
double maxEdge = 0;
/// <summary>
/// Gets the longest edge.
/// </summary>
public double LongestEdge { get { return maxEdge; } }
//
double minAspect = 0;
/// <summary>
/// Gets the shortest altitude.
/// </summary>
public double ShortestAltitude { get { return minAspect; } }
double maxAspect = 0;
/// <summary>
/// Gets the largest aspect ratio.
/// </summary>
public double LargestAspectRatio { get { return maxAspect; } }
double minArea = 0;
/// <summary>
/// Gets the smallest area.
/// </summary>
public double SmallestArea { get { return minArea; } }
double maxArea = 0;
/// <summary>
/// Gets the largest area.
/// </summary>
public double LargestArea { get { return maxArea; } }
double minAngle = 0;
/// <summary>
/// Gets the smallest angle.
/// </summary>
public double SmallestAngle { get { return minAngle; } }
double maxAngle = 0;
/// <summary>
/// Gets the largest angle.
/// </summary>
public double LargestAngle { get { return maxAngle; } }
int[] angleTable;
/// <summary>
/// Gets the angle histogram.
/// </summary>
public int[] AngleHistogram { get { return angleTable; } }
int[] minAngles;
/// <summary>
/// Gets the min angles histogram.
/// </summary>
public int[] MinAngleHistogram { get { return minAngles; } }
int[] maxAngles;
/// <summary>
/// Gets the max angles histogram.
/// </summary>
public int[] MaxAngleHistogram { get { return maxAngles; } }
double meshArea = 0;
/// <summary>
/// Gets the total mesh area.
/// </summary>
public double MeshArea { get { return meshArea; } }
#endregion
#region Private methods
private void GetAspectHistogram(Mesh mesh)
{
int[] aspecttable;
double[] ratiotable;
aspecttable = new int[16];
ratiotable = new double[]
{
1.5, 2.0, 2.5, 3.0, 4.0, 6.0, 10.0, 15.0, 25.0, 50.0,
100.0, 300.0, 1000.0, 10000.0, 100000.0, 0.0
};
Otri tri = default(Otri);
Vertex[] p = new Vertex[3];
double[] dx = new double[3], dy = new double[3];
double[] edgelength = new double[3];
double triarea;
double trilongest2;
double triminaltitude2;
double triaspect2;
int aspectindex;
int i, j, k;
tri.orient = 0;
foreach (var t in mesh.triangles)
{
tri.tri = t;
p[0] = tri.Org();
p[1] = tri.Dest();
p[2] = tri.Apex();
trilongest2 = 0.0;
for (i = 0; i < 3; i++)
{
j = plus1Mod3[i];
k = minus1Mod3[i];
dx[i] = p[j].x - p[k].x;
dy[i] = p[j].y - p[k].y;
edgelength[i] = dx[i] * dx[i] + dy[i] * dy[i];
if (edgelength[i] > trilongest2)
{
trilongest2 = edgelength[i];
}
}
//triarea = Primitives.CounterClockwise(p[0], p[1], p[2]);
triarea = Math.Abs((p[2].x - p[0].x) * (p[1].y - p[0].y) -
(p[1].x - p[0].x) * (p[2].y - p[0].y)) / 2.0;
triminaltitude2 = triarea * triarea / trilongest2;
triaspect2 = trilongest2 / triminaltitude2;
aspectindex = 0;
while ((triaspect2 > ratiotable[aspectindex] * ratiotable[aspectindex]) && (aspectindex < 15))
{
aspectindex++;
}
aspecttable[aspectindex]++;
}
}
#endregion
static readonly int[] plus1Mod3 = { 1, 2, 0 };
static readonly int[] minus1Mod3 = { 2, 0, 1 };
/// <summary>
/// Update statistics about the quality of the mesh.
/// </summary>
/// <param name="mesh"></param>
public void Update(Mesh mesh, int sampleDegrees)
{
Point[] p = new Point[3];
int k1, k2;
int degreeStep;
//sampleDegrees = 36; // sample every 5 degrees
//sampleDegrees = 45; // sample every 4 degrees
sampleDegrees = 60; // sample every 3 degrees
double[] cosSquareTable = new double[sampleDegrees / 2 - 1];
double[] dx = new double[3];
double[] dy = new double[3];
double[] edgeLength = new double[3];
double dotProduct;
double cosSquare;
double triArea;
double triLongest2;
double triMinAltitude2;
double triAspect2;
double radconst = Math.PI / sampleDegrees;
double degconst = 180.0 / Math.PI;
// New angle table
angleTable = new int[sampleDegrees];
minAngles = new int[sampleDegrees];
maxAngles = new int[sampleDegrees];
for (int i = 0; i < sampleDegrees / 2 - 1; i++)
{
cosSquareTable[i] = Math.Cos(radconst * (i + 1));
cosSquareTable[i] = cosSquareTable[i] * cosSquareTable[i];
}
for (int i = 0; i < sampleDegrees; i++)
{
angleTable[i] = 0;
}
minAspect = mesh.bounds.Width + mesh.bounds.Height;
minAspect = minAspect * minAspect;
maxAspect = 0.0;
minEdge = minAspect;
maxEdge = 0.0;
minArea = minAspect;
maxArea = 0.0;
minAngle = 0.0;
maxAngle = 2.0;
meshArea = 0.0;
bool acuteBiggest = true;
bool acuteBiggestTri = true;
double triMinAngle, triMaxAngle = 1;
foreach (var tri in mesh.triangles)
{
triMinAngle = 0; // Min angle: 0 < a < 60 degress
triMaxAngle = 1; // Max angle: 60 < a < 180 degress
p[0] = tri.vertices[0];
p[1] = tri.vertices[1];
p[2] = tri.vertices[2];
triLongest2 = 0.0;
for (int i = 0; i < 3; i++)
{
k1 = plus1Mod3[i];
k2 = minus1Mod3[i];
dx[i] = p[k1].x - p[k2].x;
dy[i] = p[k1].y - p[k2].y;
edgeLength[i] = dx[i] * dx[i] + dy[i] * dy[i];
if (edgeLength[i] > triLongest2)
{
triLongest2 = edgeLength[i];
}
if (edgeLength[i] > maxEdge)
{
maxEdge = edgeLength[i];
}
if (edgeLength[i] < minEdge)
{
minEdge = edgeLength[i];
}
}
//triarea = Primitives.CounterClockwise(p[0], p[1], p[2]);
triArea = Math.Abs((p[2].x - p[0].x) * (p[1].y - p[0].y) -
(p[1].x - p[0].x) * (p[2].y - p[0].y));
meshArea += triArea;
if (triArea < minArea)
{
minArea = triArea;
}
if (triArea > maxArea)
{
maxArea = triArea;
}
triMinAltitude2 = triArea * triArea / triLongest2;
if (triMinAltitude2 < minAspect)
{
minAspect = triMinAltitude2;
}
triAspect2 = triLongest2 / triMinAltitude2;
if (triAspect2 > maxAspect)
{
maxAspect = triAspect2;
}
for (int i = 0; i < 3; i++)
{
k1 = plus1Mod3[i];
k2 = minus1Mod3[i];
dotProduct = dx[k1] * dx[k2] + dy[k1] * dy[k2];
cosSquare = dotProduct * dotProduct / (edgeLength[k1] * edgeLength[k2]);
degreeStep = sampleDegrees / 2 - 1;
for (int j = degreeStep - 1; j >= 0; j--)
{
if (cosSquare > cosSquareTable[j])
{
degreeStep = j;
}
}
if (dotProduct <= 0.0)
{
angleTable[degreeStep]++;
if (cosSquare > minAngle)
{
minAngle = cosSquare;
}
if (acuteBiggest && (cosSquare < maxAngle))
{
maxAngle = cosSquare;
}
// Update min/max angle per triangle
if (cosSquare > triMinAngle)
{
triMinAngle = cosSquare;
}
if (acuteBiggestTri && (cosSquare < triMaxAngle))
{
triMaxAngle = cosSquare;
}
}
else
{
angleTable[sampleDegrees - degreeStep - 1]++;
if (acuteBiggest || (cosSquare > maxAngle))
{
maxAngle = cosSquare;
acuteBiggest = false;
}
// Update max angle for (possibly non-acute) triangle
if (acuteBiggestTri || (cosSquare > triMaxAngle))
{
triMaxAngle = cosSquare;
acuteBiggestTri = false;
}
}
}
// Update min angle histogram
degreeStep = sampleDegrees / 2 - 1;
for (int j = degreeStep - 1; j >= 0; j--)
{
if (triMinAngle > cosSquareTable[j])
{
degreeStep = j;
}
}
minAngles[degreeStep]++;
// Update max angle histogram
degreeStep = sampleDegrees / 2 - 1;
for (int j = degreeStep - 1; j >= 0; j--)
{
if (triMaxAngle > cosSquareTable[j])
{
degreeStep = j;
}
}
if (acuteBiggestTri)
{
maxAngles[degreeStep]++;
}
else
{
maxAngles[sampleDegrees - degreeStep - 1]++;
}
acuteBiggestTri = true;
}
minEdge = Math.Sqrt(minEdge);
maxEdge = Math.Sqrt(maxEdge);
minAspect = Math.Sqrt(minAspect);
maxAspect = Math.Sqrt(maxAspect);
minArea *= 0.5;
maxArea *= 0.5;
if (minAngle >= 1.0)
{
minAngle = 0.0;
}
else
{
minAngle = degconst * Math.Acos(Math.Sqrt(minAngle));
}
if (maxAngle >= 1.0)
{
maxAngle = 180.0;
}
else
{
if (acuteBiggest)
{
maxAngle = degconst * Math.Acos(Math.Sqrt(maxAngle));
}
else
{
maxAngle = 180.0 - degconst * Math.Acos(Math.Sqrt(maxAngle));
}
}
}
/// <summary>
/// Compute angle information for given triangle.
/// </summary>
/// <param name="triangle">The triangle to check.</param>
/// <param name="data">Array of doubles (length 6).</param>
/// <remarks>
/// On return, the squared cosines of the minimum and maximum angle will
/// be stored at position data[0] and data[1] respectively.
/// If the triangle was obtuse, data[2] will be set to -1 and maximum angle
/// is computed as (pi - acos(sqrt(data[1]))).
/// </remarks>
internal static void ComputeAngles(ITriangle triangle, double[] data)
{
double min = 0.0;
double max = 1.0;
var va = triangle.GetVertex(0);
var vb = triangle.GetVertex(1);
var vc = triangle.GetVertex(2);
double dxa = vb.x - vc.x;
double dya = vb.y - vc.y;
double lena = dxa * dxa + dya * dya;
double dxb = vc.x - va.x;
double dyb = vc.y - va.y;
double lenb = dxb * dxb + dyb * dyb;
double dxc = va.x - vb.x;
double dyc = va.y - vb.y;
double lenc = dxc * dxc + dyc * dyc;
// Dot products.
double dota = data[0] = dxb * dxc + dyb * dyc;
double dotb = data[1] = dxc * dxa + dyc * dya;
double dotc = data[2] = dxa * dxb + dya * dyb;
// Squared cosines.
data[3] = (dota * dota) / (lenb * lenc);
data[4] = (dotb * dotb) / (lenc * lena);
data[5] = (dotc * dotc) / (lena * lenb);
// The sign of the dot product will tell us, if the angle is
// acute (value < 0) or obtuse (value > 0).
bool acute = true;
double cos, dot;
for (int i = 0; i < 3; i++)
{
dot = data[i];
cos = data[3 + i];
if (dot <= 0.0)
{
if (cos > min)
{
min = cos;
}
if (acute && (cos < max))
{
max = cos;
}
}
else
{
// Update max angle for (possibly non-acute) triangle
if (acute || (cos > max))
{
max = cos;
acute = false;
}
}
}
data[0] = min;
data[1] = max;
data[2] = acute ? 1.0 : -1.0;
}
}
}

View file

@ -0,0 +1,11 @@
fileFormatVersion: 2
guid: 1af64e13882d6445bb141ef2f2810b17
MonoImporter:
externalObjects: {}
serializedVersion: 2
defaultReferences: []
executionOrder: 0
icon: {instanceID: 0}
userData:
assetBundleName:
assetBundleVariant:

View file

@ -0,0 +1,427 @@
// -----------------------------------------------------------------------
// <copyright file="TriangleQuadTree.cs" company="">
// Original code by Frank Dockhorn, [not available anymore: http://sourceforge.net/projects/quadtreesim/]
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace UnityEngine.U2D.Animation.TriangleNet
.Tools
{
using System.Collections.Generic;
using System.Linq;
using Animation.TriangleNet.Geometry;
/// <summary>
/// A Quadtree implementation optimized for triangles.
/// </summary>
internal class TriangleQuadTree
{
QuadNode root;
internal ITriangle[] triangles;
internal int sizeBound;
internal int maxDepth;
/// <summary>
/// Initializes a new instance of the <see cref="TriangleQuadTree" /> class.
/// </summary>
/// <param name="mesh">Mesh containing triangles.</param>
/// <param name="maxDepth">The maximum depth of the tree.</param>
/// <param name="sizeBound">The maximum number of triangles contained in a leaf.</param>
/// <remarks>
/// The quadtree does not track changes of the mesh. If a mesh is refined or
/// changed in any other way, a new quadtree has to be built to make the point
/// location work.
///
/// A node of the tree will be split, if its level if less than the max depth parameter
/// AND the number of triangles in the node is greater than the size bound.
/// </remarks>
public TriangleQuadTree(Mesh mesh, int maxDepth = 10, int sizeBound = 10)
{
this.maxDepth = maxDepth;
this.sizeBound = sizeBound;
triangles = mesh.Triangles.ToArray();
int currentDepth = 0;
root = new QuadNode(mesh.Bounds, this, true);
root.CreateSubRegion(++currentDepth);
}
public ITriangle Query(double x, double y)
{
var point = new Point(x, y);
var indices = root.FindTriangles(point);
foreach (var i in indices)
{
var tri = this.triangles[i];
if (IsPointInTriangle(point, tri.GetVertex(0), tri.GetVertex(1), tri.GetVertex(2)))
{
return tri;
}
}
return null;
}
/// <summary>
/// Test, if a given point lies inside a triangle.
/// </summary>
/// <param name="p">Point to locate.</param>
/// <param name="t0">Corner point of triangle.</param>
/// <param name="t1">Corner point of triangle.</param>
/// <param name="t2">Corner point of triangle.</param>
/// <returns>True, if point is inside or on the edge of this triangle.</returns>
internal static bool IsPointInTriangle(Point p, Point t0, Point t1, Point t2)
{
// TODO: no need to create new Point instances here
Point d0 = new Point(t1.x - t0.x, t1.y - t0.y);
Point d1 = new Point(t2.x - t0.x, t2.y - t0.y);
Point d2 = new Point(p.x - t0.x, p.y - t0.y);
// crossproduct of (0, 0, 1) and d0
Point c0 = new Point(-d0.y, d0.x);
// crossproduct of (0, 0, 1) and d1
Point c1 = new Point(-d1.y, d1.x);
// Linear combination d2 = s * d0 + v * d1.
//
// Multiply both sides of the equation with c0 and c1
// and solve for s and v respectively
//
// s = d2 * c1 / d0 * c1
// v = d2 * c0 / d1 * c0
double s = DotProduct(d2, c1) / DotProduct(d0, c1);
double v = DotProduct(d2, c0) / DotProduct(d1, c0);
if (s >= 0 && v >= 0 && ((s + v) <= 1))
{
// Point is inside or on the edge of this triangle.
return true;
}
return false;
}
internal static double DotProduct(Point p, Point q)
{
return p.x * q.x + p.y * q.y;
}
/// <summary>
/// A node of the quadtree.
/// </summary>
class QuadNode
{
const int SW = 0;
const int SE = 1;
const int NW = 2;
const int NE = 3;
const double EPS = 1e-6;
static readonly byte[] BITVECTOR = { 0x1, 0x2, 0x4, 0x8 };
Rectangle bounds;
Point pivot;
TriangleQuadTree tree;
QuadNode[] regions;
List<int> triangles;
byte bitRegions;
public QuadNode(Rectangle box, TriangleQuadTree tree)
: this(box, tree, false)
{
}
public QuadNode(Rectangle box, TriangleQuadTree tree, bool init)
{
this.tree = tree;
this.bounds = new Rectangle(box.Left, box.Bottom, box.Width, box.Height);
this.pivot = new Point((box.Left + box.Right) / 2, (box.Bottom + box.Top) / 2);
this.bitRegions = 0;
this.regions = new QuadNode[4];
this.triangles = new List<int>();
if (init)
{
int count = tree.triangles.Length;
// Allocate memory upfront
triangles.Capacity = count;
for (int i = 0; i < count; i++)
{
triangles.Add(i);
}
}
}
public List<int> FindTriangles(Point searchPoint)
{
int region = FindRegion(searchPoint);
if (regions[region] == null)
{
return triangles;
}
return regions[region].FindTriangles(searchPoint);
}
public void CreateSubRegion(int currentDepth)
{
// The four sub regions of the quad tree
// +--------------+
// | nw 2 | ne 3 |
// |------+pivot--|
// | sw 0 | se 1 |
// +--------------+
Rectangle box;
var width = bounds.Right - pivot.x;
var height = bounds.Top - pivot.y;
// 1. region south west
box = new Rectangle(bounds.Left, bounds.Bottom, width, height);
regions[0] = new QuadNode(box, tree);
// 2. region south east
box = new Rectangle(pivot.x, bounds.Bottom, width, height);
regions[1] = new QuadNode(box, tree);
// 3. region north west
box = new Rectangle(bounds.Left, pivot.y, width, height);
regions[2] = new QuadNode(box, tree);
// 4. region north east
box = new Rectangle(pivot.x, pivot.y, width, height);
regions[3] = new QuadNode(box, tree);
Point[] triangle = new Point[3];
// Find region for every triangle vertex
foreach (var index in triangles)
{
ITriangle tri = tree.triangles[index];
triangle[0] = tri.GetVertex(0);
triangle[1] = tri.GetVertex(1);
triangle[2] = tri.GetVertex(2);
AddTriangleToRegion(triangle, index);
}
for (int i = 0; i < 4; i++)
{
if (regions[i].triangles.Count > tree.sizeBound && currentDepth < tree.maxDepth)
{
regions[i].CreateSubRegion(currentDepth + 1);
}
}
}
void AddTriangleToRegion(Point[] triangle, int index)
{
bitRegions = 0;
if (TriangleQuadTree.IsPointInTriangle(pivot, triangle[0], triangle[1], triangle[2]))
{
AddToRegion(index, SW);
AddToRegion(index, SE);
AddToRegion(index, NW);
AddToRegion(index, NE);
return;
}
FindTriangleIntersections(triangle, index);
if (bitRegions == 0)
{
// we didn't find any intersection so we add this triangle to a point's region
int region = FindRegion(triangle[0]);
regions[region].triangles.Add(index);
}
}
void FindTriangleIntersections(Point[] triangle, int index)
{
// PLEASE NOTE:
// Handling of component comparison is tightly associated with the implementation
// of the findRegion() function. That means when the point to be compared equals
// the pivot point the triangle must be put at least into region 2.
//
// Linear equations are in parametric form.
// pivot.x = triangle[0].x + t * (triangle[1].x - triangle[0].x)
// pivot.y = triangle[0].y + t * (triangle[1].y - triangle[0].y)
int k = 2;
double dx, dy;
// Iterate through all triangle laterals and find bounding box intersections
for (int i = 0; i < 3; k = i++)
{
dx = triangle[i].x - triangle[k].x;
dy = triangle[i].y - triangle[k].y;
if (dx != 0.0)
{
FindIntersectionsWithX(dx, dy, triangle, index, k);
}
if (dy != 0.0)
{
FindIntersectionsWithY(dx, dy, triangle, index, k);
}
}
}
void FindIntersectionsWithX(double dx, double dy, Point[] triangle, int index, int k)
{
double t;
// find intersection with plane x = m_pivot.dX
t = (pivot.x - triangle[k].x) / dx;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
double yComponent = triangle[k].y + t * dy;
if (yComponent < pivot.y && yComponent >= bounds.Bottom)
{
AddToRegion(index, SW);
AddToRegion(index, SE);
}
else if (yComponent <= bounds.Top)
{
AddToRegion(index, NW);
AddToRegion(index, NE);
}
}
// find intersection with plane x = m_boundingBox[0].dX
t = (bounds.Left - triangle[k].x) / dx;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
double yComponent = triangle[k].y + t * dy;
if (yComponent < pivot.y && yComponent >= bounds.Bottom)
{
AddToRegion(index, SW);
}
else if (yComponent <= bounds.Top) // TODO: check && yComponent >= pivot.Y
{
AddToRegion(index, NW);
}
}
// find intersection with plane x = m_boundingBox[1].dX
t = (bounds.Right - triangle[k].x) / dx;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
double yComponent = triangle[k].y + t * dy;
if (yComponent < pivot.y && yComponent >= bounds.Bottom)
{
AddToRegion(index, SE);
}
else if (yComponent <= bounds.Top)
{
AddToRegion(index, NE);
}
}
}
void FindIntersectionsWithY(double dx, double dy, Point[] triangle, int index, int k)
{
double t, xComponent;
// find intersection with plane y = m_pivot.dY
t = (pivot.y - triangle[k].y) / dy;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
xComponent = triangle[k].x + t * dx;
if (xComponent > pivot.x && xComponent <= bounds.Right)
{
AddToRegion(index, SE);
AddToRegion(index, NE);
}
else if (xComponent >= bounds.Left)
{
AddToRegion(index, SW);
AddToRegion(index, NW);
}
}
// find intersection with plane y = m_boundingBox[0].dY
t = (bounds.Bottom - triangle[k].y) / dy;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
xComponent = triangle[k].x + t * dx;
if (xComponent > pivot.x && xComponent <= bounds.Right)
{
AddToRegion(index, SE);
}
else if (xComponent >= bounds.Left)
{
AddToRegion(index, SW);
}
}
// find intersection with plane y = m_boundingBox[1].dY
t = (bounds.Top - triangle[k].y) / dy;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
xComponent = triangle[k].x + t * dx;
if (xComponent > pivot.x && xComponent <= bounds.Right)
{
AddToRegion(index, NE);
}
else if (xComponent >= bounds.Left)
{
AddToRegion(index, NW);
}
}
}
int FindRegion(Point point)
{
int b = 2;
if (point.y < pivot.y)
{
b = 0;
}
if (point.x > pivot.x)
{
b++;
}
return b;
}
void AddToRegion(int index, int region)
{
//if (!(m_bitRegions & BITVECTOR[region]))
if ((bitRegions & BITVECTOR[region]) == 0)
{
regions[region].triangles.Add(index);
bitRegions |= BITVECTOR[region];
}
}
}
}
}

View file

@ -0,0 +1,11 @@
fileFormatVersion: 2
guid: 1a6b4c973006a4e1eb35cb87a8d6d385
MonoImporter:
externalObjects: {}
serializedVersion: 2
defaultReferences: []
executionOrder: 0
icon: {instanceID: 0}
userData:
assetBundleName:
assetBundleVariant:

View file

@ -0,0 +1,372 @@
// -----------------------------------------------------------------------
// <copyright file="VertexSorter.cs" company="">
// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace UnityEngine.U2D.Animation.TriangleNet
.Tools
{
using System;
using Animation.TriangleNet.Geometry;
/// <summary>
/// Sort an array of points using quicksort.
/// </summary>
internal class VertexSorter
{
private const int RANDOM_SEED = 57113;
Random rand;
Vertex[] points;
VertexSorter(Vertex[] points, int seed)
{
this.points = points;
this.rand = new Random(seed);
}
/// <summary>
/// Sorts the given vertex array by x-coordinate.
/// </summary>
/// <param name="array">The vertex array.</param>
/// <param name="seed">Random seed used for pivoting.</param>
internal static void Sort(Vertex[] array, int seed = RANDOM_SEED)
{
var qs = new VertexSorter(array, seed);
qs.QuickSort(0, array.Length - 1);
}
/// <summary>
/// Impose alternating cuts on given vertex array.
/// </summary>
/// <param name="array">The vertex array.</param>
/// <param name="length">The number of vertices to sort.</param>
/// <param name="seed">Random seed used for pivoting.</param>
internal static void Alternate(Vertex[] array, int length, int seed = RANDOM_SEED)
{
var qs = new VertexSorter(array, seed);
int divider = length >> 1;
// Re-sort the array of vertices to accommodate alternating cuts.
if (length - divider >= 2)
{
if (divider >= 2)
{
qs.AlternateAxes(0, divider - 1, 1);
}
qs.AlternateAxes(divider, length - 1, 1);
}
}
#region Quicksort
/// <summary>
/// Sort an array of vertices by x-coordinate, using the y-coordinate as a secondary key.
/// </summary>
/// <param name="left"></param>
/// <param name="right"></param>
/// <remarks>
/// Uses quicksort. Randomized O(n log n) time. No, I did not make any of
/// the usual quicksort mistakes.
/// </remarks>
private void QuickSort(int left, int right)
{
int oleft = left;
int oright = right;
int arraysize = right - left + 1;
int pivot;
double pivotx, pivoty;
Vertex temp;
var array = this.points;
if (arraysize < 32)
{
// Insertion sort
for (int i = left + 1; i <= right; i++)
{
var a = array[i];
int j = i - 1;
while (j >= left && (array[j].x > a.x || (array[j].x == a.x && array[j].y > a.y)))
{
array[j + 1] = array[j];
j--;
}
array[j + 1] = a;
}
return;
}
// Choose a random pivot to split the array.
pivot = rand.Next(left, right);
pivotx = array[pivot].x;
pivoty = array[pivot].y;
// Split the array.
left--;
right++;
while (left < right)
{
// Search for a vertex whose x-coordinate is too large for the left.
do
{
left++;
}
while ((left <= right) && ((array[left].x < pivotx) ||
((array[left].x == pivotx) && (array[left].y < pivoty))));
// Search for a vertex whose x-coordinate is too small for the right.
do
{
right--;
}
while ((left <= right) && ((array[right].x > pivotx) ||
((array[right].x == pivotx) && (array[right].y > pivoty))));
if (left < right)
{
// Swap the left and right vertices.
temp = array[left];
array[left] = array[right];
array[right] = temp;
}
}
if (left > oleft)
{
// Recursively sort the left subset.
QuickSort(oleft, left);
}
if (oright > right + 1)
{
// Recursively sort the right subset.
QuickSort(right + 1, oright);
}
}
#endregion
#region Alternate axes
/// <summary>
/// Sorts the vertices as appropriate for the divide-and-conquer algorithm with
/// alternating cuts.
/// </summary>
/// <param name="left"></param>
/// <param name="right"></param>
/// <param name="axis"></param>
/// <remarks>
/// Partitions by x-coordinate if axis == 0; by y-coordinate if axis == 1.
/// For the base case, subsets containing only two or three vertices are
/// always sorted by x-coordinate.
/// </remarks>
private void AlternateAxes(int left, int right, int axis)
{
int size = right - left + 1;
int divider = size >> 1;
if (size <= 3)
{
// Recursive base case: subsets of two or three vertices will be
// handled specially, and should always be sorted by x-coordinate.
axis = 0;
}
// Partition with a horizontal or vertical cut.
if (axis == 0)
{
VertexMedianX(left, right, left + divider);
}
else
{
VertexMedianY(left, right, left + divider);
}
// Recursively partition the subsets with a cross cut.
if (size - divider >= 2)
{
if (divider >= 2)
{
AlternateAxes(left, left + divider - 1, 1 - axis);
}
AlternateAxes(left + divider, right, 1 - axis);
}
}
/// <summary>
/// An order statistic algorithm, almost. Shuffles an array of vertices so that the
/// first 'median' vertices occur lexicographically before the remaining vertices.
/// </summary>
/// <param name="left"></param>
/// <param name="right"></param>
/// <param name="median"></param>
/// <remarks>
/// Uses the x-coordinate as the primary key. Very similar to the QuickSort()
/// procedure, but runs in randomized linear time.
/// </remarks>
private void VertexMedianX(int left, int right, int median)
{
int arraysize = right - left + 1;
int oleft = left, oright = right;
int pivot;
double pivot1, pivot2;
Vertex temp;
var array = this.points;
if (arraysize == 2)
{
// Recursive base case.
if ((array[left].x > array[right].x) ||
((array[left].x == array[right].x) &&
(array[left].y > array[right].y)))
{
temp = array[right];
array[right] = array[left];
array[left] = temp;
}
return;
}
// Choose a random pivot to split the array.
pivot = rand.Next(left, right);
pivot1 = array[pivot].x;
pivot2 = array[pivot].y;
left--;
right++;
while (left < right)
{
// Search for a vertex whose x-coordinate is too large for the left.
do
{
left++;
}
while ((left <= right) && ((array[left].x < pivot1) ||
((array[left].x == pivot1) && (array[left].y < pivot2))));
// Search for a vertex whose x-coordinate is too small for the right.
do
{
right--;
}
while ((left <= right) && ((array[right].x > pivot1) ||
((array[right].x == pivot1) && (array[right].y > pivot2))));
if (left < right)
{
// Swap the left and right vertices.
temp = array[left];
array[left] = array[right];
array[right] = temp;
}
}
// Unlike in vertexsort(), at most one of the following conditionals is true.
if (left > median)
{
// Recursively shuffle the left subset.
VertexMedianX(oleft, left - 1, median);
}
if (right < median - 1)
{
// Recursively shuffle the right subset.
VertexMedianX(right + 1, oright, median);
}
}
/// <summary>
/// An order statistic algorithm, almost. Shuffles an array of vertices so that
/// the first 'median' vertices occur lexicographically before the remaining vertices.
/// </summary>
/// <param name="left"></param>
/// <param name="right"></param>
/// <param name="median"></param>
/// <remarks>
/// Uses the y-coordinate as the primary key. Very similar to the QuickSort()
/// procedure, but runs in randomized linear time.
/// </remarks>
private void VertexMedianY(int left, int right, int median)
{
int arraysize = right - left + 1;
int oleft = left, oright = right;
int pivot;
double pivot1, pivot2;
Vertex temp;
var array = this.points;
if (arraysize == 2)
{
// Recursive base case.
if ((array[left].y > array[right].y) ||
((array[left].y == array[right].y) &&
(array[left].x > array[right].x)))
{
temp = array[right];
array[right] = array[left];
array[left] = temp;
}
return;
}
// Choose a random pivot to split the array.
pivot = rand.Next(left, right);
pivot1 = array[pivot].y;
pivot2 = array[pivot].x;
left--;
right++;
while (left < right)
{
// Search for a vertex whose x-coordinate is too large for the left.
do
{
left++;
}
while ((left <= right) && ((array[left].y < pivot1) ||
((array[left].y == pivot1) && (array[left].x < pivot2))));
// Search for a vertex whose x-coordinate is too small for the right.
do
{
right--;
}
while ((left <= right) && ((array[right].y > pivot1) ||
((array[right].y == pivot1) && (array[right].x > pivot2))));
if (left < right)
{
// Swap the left and right vertices.
temp = array[left];
array[left] = array[right];
array[right] = temp;
}
}
// Unlike in QuickSort(), at most one of the following conditionals is true.
if (left > median)
{
// Recursively shuffle the left subset.
VertexMedianY(oleft, left - 1, median);
}
if (right < median - 1)
{
// Recursively shuffle the right subset.
VertexMedianY(right + 1, oright, median);
}
}
#endregion
}
}

View file

@ -0,0 +1,11 @@
fileFormatVersion: 2
guid: 5f02453d23ed9451aa6e2ccb3efc0597
MonoImporter:
externalObjects: {}
serializedVersion: 2
defaultReferences: []
executionOrder: 0
icon: {instanceID: 0}
userData:
assetBundleName:
assetBundleVariant: