我有一个三角形列表,每个三角形节点都有 XYZ 坐标。我需要一个节点列表和一个三角形节点索引列表。关闭的节点应该获得相同的索引。 我的 C# 代码是:
// ensureNodeAt
// return node index. Avoid duplicates.
public int ensureNodeAt(XYZ p, double eps_dist = 1e-9)
{
// backwards - most likely the one we seek is at the back
for (int i = nodes.Count - 1; i >= 0; --i)
{
XYZ v = nodes[i].v - p;
// dx = v.X - p.X;
// dy = v.Y - p.Y;
// dz = v.Z - p.Z;
// dx = dx * dx + dy * dy + dz * dz;
if (v.X * v.X + v.Y * v.Y + v.Z * v.Z < eps_dist)
{ return i; }
}
// Must add a new node
nodes.Add(p);
return nodes.Count - 1;
}
void work(List<XYZ> triangles){
tri = new List<Tri>(triangles.Count / 3); // triangle indices
nodes = new List<XYZ>(triangles.Count); // triangle nodes
for (int i = 0; i < triangles.Count; i += 3)
{
var a = triangles[i];
var b = triangles[i + 1];
var c = triangles[i + 2];
tri.Add(new Tri(ensureNodeAt(a, eps_dist), ensureNodeAt(b, eps_dist), ensureNodeAt(c, eps_dist)));
}
}
有什么办法可以让它更快或并行运行吗?我有大约 900k 个节点。
为了将节点标记为重复项,它们必须在
eps_dist
内包含 X。因此,如果您将点按 X 排序,那么您就不需要检查这么多重复项。
您可以扫描已排序的集合,直到
nodes[i].v.X >= p.X-eps_dist
,然后开始测试重复项的距离,然后在 nodes[i].v.X > p+eps_dist
时停止。
Y或Z也是如此;因此,如果您有理由认为点在这些方向之一上会更加分离,则可以使用其中之一而不是 X。
我是这样解决的。我的数据主要是 X、Y 数据,带有轻微的 Z 偏移,因此按 x 范围对点进行分组对我来说很有效:
namespace Whatever
{
using System;
using System.Collections.Generic;
using WhateverXYZ; // class XYZ
public class SpacePartitioning
{
// PUBLIC
public double eps = 1e-3; // max error for neighbour nodes, so they are "the same"
public List<XYZ> nodes;
// PRIVATE
// A Space is a list of node indices within a range xmin/xmax.
// nodes very close to xmin/xmax are in both neighbour spaces,
// so finding a close node is ensured for both neighbour spaces.
class Space
{
public Space(int perSpaceNodes) { xmin = 0; xmax = 0; indices = new List<int>(perSpaceNodes); }
public double xmin, xmax; // the bounds of the space are extendes by +/- eps to check if a point is inside
public List<int> indices;
public double Width() { return xmax - xmin; }
}
private List<Space> spaces;
int perSpaceNodes = 128; // desired masimum number of nodes per space
const double minSpaceWidth = 2.0; // [ft] the minimum width of a space, even if it has lots of nodes
// C'tor
public SpacePartitioning(int nodeCountPreallocation = 1000000)
{
Clear(nodeCountPreallocation);
}
// Clear
public void Clear(int nodeCountPreallocation = 1000000)
{
// keep the number of nodes per spaces lower than the number of nodes
int sideLength = (int)Math.Sqrt(nodeCountPreallocation);
perSpaceNodes = Math.Max(128, Math.Min(4096, sideLength / 8));
nodes = new List<XYZ>(nodeCountPreallocation);
spaces = new List<Space>(64);
}
// PUBLIC
// Find or add a node to the list of nodes.
// Specify this.eps before starting to add nodes.
// Returns the index in m_nodes[]
public int FindOrAddNode(XYZ p)
{
double eps2 = eps * eps;
EnsureSpacesFor(p);
foreach (var space in spaces)
{
if (IsPointInSpace(space, p))
{
foreach (var ix in space.indices)
{
if (dist2(nodes[ix], p) < eps2) { return ix; }
}
// Point must go in this space
return InsertPointToSpaces(p);
}
}
int can_t_be = 1;
return -1;
}
// squared distance of two points
private static double dist2(XYZ p1, XYZ p2)
{
XYZ d = p2 - p1;
return d.X * d.X + d.Y * d.Y + d.Z * d.Z;
}
// Split a space with many nodes into two spaces with only half of the nodes, each.
private void SplitLargeSpace(Space spc)
{
if (spc.Width() < minSpaceWidth) { return; }
// find average point center
double xavg = 0.0;
foreach (int ix in spc.indices) { xavg += nodes[ix].X; }
xavg /= (double)spc.indices.Count;
// split space
Space sp2 = new Space(perSpaceNodes);
sp2.xmin = xavg;
sp2.xmax = spc.xmax;
spc.xmax = xavg;
// reasign nodes to the corresponding spaces
for (int i = 0; i < spc.indices.Count; i++)
{
int ix = spc.indices[i];
XYZ p = nodes[ix];
// both can be possible!
// remove from old space
if (!IsPointInSpace(spc, p))
{
spc.indices.RemoveAt(i);
--i;
}
// add in new space
if (IsPointInSpace(sp2, p))
{
sp2.indices.Add(ix);
}
}
spaces.Add(sp2);
}
// Adds point to nodes list and then puts the reference in any space
private int InsertPointToSpaces(XYZ p)
{
// add new point to the point list
int ix = nodes.Count;
nodes.Add(p);
// insert new point to spaces
foreach (var space in spaces)
{
if (IsPointInSpace(space, p))
{
space.indices.Add(ix);
}
}
// split large spaces
for (int i = 0; i < spaces.Count; i++)
{
var space = spaces[i];
if (space.indices.Count >= perSpaceNodes)
{
SplitLargeSpace(space);
}
}
return ix;
}
// Esnure, the spaces containers are wide enough to capture the point
private void EnsureSpacesFor(XYZ p)
{
const double largeWidth = 10000.0; // "infinite" range for new spaces
// empty space
if (spaces.Count == 0)
{
var ns = new Space(perSpaceNodes);
ns.xmin = p.X - largeWidth;
ns.xmax = p.X + largeWidth;
spaces.Add(ns);
}
// left of all spaces
if (spaces[0].xmin > p.X + eps)
{
var ns = new Space(perSpaceNodes);
ns.xmin = p.X - largeWidth;
ns.xmax = spaces[0].xmin;
}
// right of all spaces
if (spaces[spaces.Count - 1].xmax < p.X - eps)
{
var ns = new Space(perSpaceNodes);
ns.xmin = spaces[spaces.Count - 1].xmax;
ns.xmax = p.X + largeWidth;
}
}
// Is a point within the space (including some extra space for nodes close to the boundary of a space)
bool IsPointInSpace(Space space, XYZ p) { return p.X >= space.xmin - eps && p.X <= space.xmax + eps; }
}
}