在二维点集中寻找洞? [英] Finding holes in 2d point sets?
问题描述
我有一组 2d 点.它们是标准笛卡尔网格系统(在本例中为 UTM 区域)上的 X,Y 坐标.我需要找到该点集中的漏洞，最好具有设置找到这些漏洞的算法的灵敏度的能力.通常，这些点集非常密集，但有些点集的密度要小得多.
I have a set of 2d points. They are X,Y coordinates on a standard Cartesian grid system(in this case a UTM zone). I need to find the holes in that point set preferably with some ability to set the sensitivity of the algorithm that finds these holes. Typically these point sets are very dense but some can be much less dense.
它们是什么，如果有帮助的话，是对田间土壤进行采样的点，以获得农业人员显然认为有用的各种特性.有时在这些点样本中有巨大的岩石或沼泽地或充满小湖泊和池塘.
What they are, if it helps any, are points at which the soil in the field has been sampled for various properties that people in agriculture apparently find useful. Sometimes in these point samples there are giant rocks or swampy places or full on little lakes and ponds.
他们希望我从点集中找到大致定义这些孔的凹多边形.
From the point set they want me to find the concave polygon that roughly defines these holes.
我已经编写了算法来查找外部凹边界多边形，并允许他们为算法形成的多边形的粗糙程度或平滑程度设置敏感度.在那次运行之后，他们想找到洞并将这些洞作为凹多边形给他们，我猜在某些情况下可能是凸多边形，但这将是边缘情况而不是常态.
I already wrote the algorithm that finds the exterior concave boundary polygon and allows them to set sensitivity for how rough or smooth the polygon is that is formed by the algorithm. After that runs they would like to find holes and have those holes given to them as a concave polygon which I guess in some cases might be convex but that will be the edge case not the norm.
问题是我听说过的关于这个主题的唯一论文是天文学家写的，他们想在太空中找到大片的空洞，而我从来没有找到他们的一篇论文，其中显示了实际的算法将它们视为粗略的概念描述之外的任何内容.
The problem is the only papers I have ever heard of on this subject are ones done by astronomers who want to find big pockets of emptiness out in space and I have never been able to find one of their papers with an actual algorithm shown in them as anything other than a rough conceptual description.
我尝试过谷歌和各种学术论文搜索等，但到目前为止我还没有找到太多有用的东西.这让我想知道这种问题是否有一个名称而我不知道，所以我正在寻找错误的东西或什么?
I have tried Google and various scholarly paper searches etc. but I haven’t found much that is useful so far. Which makes me wonder if maybe there is a name for this kind of problem and I don’t know it so I am searching for the wrong thing or something?
所以经过冗长的解释之后，我的问题是:有没有人知道我应该搜索什么来找到最好使用定义明确的算法的论文，或者是否有人知道他们可以指出我的算法可以做到这一点?
So after that long winded explanation, my question is: Does anyone know what I should be searching for to find papers on this preferably with well defined algorithms or does somebody know an algorithm that will do this that they can point me to?
任何能帮助我解决这个问题的方法都会非常有用并且非常感谢，即使是正确的搜索词组也会出现潜在的算法或论文.
Anything to help me solve this problem would be very useful and greatly appreciated, even just correct search phrases that will turn up the potential algorithms or papers.
这将使用 C# 语言实现，但从 Mathematica 包到 MATLAB 或 ASM、C、C++、Python、Java 或 MathCAD 等任何示例都可以，只要示例中没有一些调用类似于 FindTheHole 等.其中 FindTheHole 未定义或专用于实现软件，例如MathCAD 示例通常有很多.
The language this will be implemented in will be C# but examples in anything from Mathematica packages to MATLAB or ASM, C, C++, Python, Java or MathCAD etc. would be fine so long as the example doesn’t have some calls in it that go to things like FindTheHole etc. Where FindTheHole isn’t defined or is proprietary to the implementing software e.g. MathCAD examples typically have a lot of that.
以下是实际点集的两个示例，一个密集的和一个稀疏的以及我需要找到的区域:
Below are two examples of actual point sets, one dense and one sparse and the areas in them I would need to find:
推荐答案
像这样的位图+矢量方法怎么样:
获取点云区域覆盖的边界框
如果未知，请执行此操作.它应该很简单 O(N) 循环遍历所有点.
Do this if it is not already known. It should be simple O(N) cycle through all points.
创建区域的map[N][N]
它是区域的位图"，便于数据密度计算.只需从 area(x,y) -> 创建投影map[i][j] 例如使用简单的比例尺.网格大小N也是输出的准确度，必须大于平均点距离！！！所以map内的每个单元格[][]覆盖区域至少有一个点(如果不在洞区域).
It is a 'bitmap' of the area for easy data density computation. Just create projection from area(x,y) -> map[i][j] for example with simple scale. Grid size N is also the accuracy of the output and must be bigger then average point distance !!! so each cell inside map[][] covers area with at least one point (if not in hole area).
计算map[][]
像馅饼一样简单，只需将 map[][].cnt(点的计数器)清除为 zero 并通过简单的 O(N) 计算> 循环 where do map[i][j].cnt++ for all points(x,y)
Easy as pie just clear map[][].cnt (counter of points) to zero and compute by simple O(N) cycle where do map[i][j].cnt++ for all points(x,y)
创建未使用区域列表(map[][].cnt==0)或(map[][].cnt<=treshold)
为了简单起见，我用水平线和垂直线来做
I do it by Horizontal and Vertical lines for simplicity
分割输出
只需将同一孔的线组合在一起(相交的...矢量方法)，也可以通过洪水填充(位图方法)在项目符号#4中完成
Just group lines of the same hole together (intersecting ones ... vector approach) and also can be done in bullet #4 by flood fill (bitmap approach)
多边形化输出
取同一孔/组的H、V线的所有边缘点并创建多边形(对它们进行排序，使它们的连接不相交).有很多关于这个的库、算法和源代码.
Take all edge points of H,V lines of the same hole/group and create polygon (sort them so their connection does not intersect anything). There are lot of libs,algorithms and source code about this.
我的这种方法的源代码:
void main_compute(int N)
{
// cell storage for density computation
struct _cell
{
double x0,x1,y0,y1; // bounding area of points inside cell
int cnt; // points inside cell
_cell(){}; _cell(_cell& a){ *this=a; }; ~_cell(){}; _cell* operator = (const _cell *a) { *this=*a; return this; }; /*_cell* operator = (const _cell &a) { ...copy... return this; };*/
};
// line storage for hole area
struct _line
{
double x0,y0,x1,y1; // line edge points
int id; // id of hole for segmentation/polygonize
int i0,i1,j0,j1; // index in map[][]
_line(){}; _line(_line& a){ *this=a; }; ~_line(){}; _line* operator = (const _line *a) { *this=*a; return this; }; /*_line* operator = (const _line &a) { ...copy... return this; };*/
};
int i,j,k,M=N*N; // M = max N^2 but usualy is much much less so dynamic list will be better
double mx,my; // scale to map
_cell *m; // cell ptr
glview2D::_pnt *p; // point ptr
double x0,x1,y0,y1; // used area (bounding box)
_cell **map=NULL; // cell grid
_line *lin=NULL; // temp line list for hole segmentation
int lins=0; // actual usage/size of lin[M]
// scan point cloud for bounding box (if it is known then skip it)
p=&view.pnt[0];
x0=p->p[0]; x1=x0;
y0=p->p[1]; y1=y0;
for (i=0;i<view.pnt.num;i++)
{
p=&view.pnt[i];
if (x0>p->p[0]) x0=p->p[0];
if (x1<p->p[0]) x1=p->p[0];
if (y0>p->p[1]) y0=p->p[1];
if (y1<p->p[1]) y1=p->p[1];
}
// compute scale for coordinate to map index conversion
mx=double(N)/(x1-x0); // add avoidance of division by zero if empty point cloud !!!
my=double(N)/(y1-y0);
// dynamic allocation of map[N][N],lin[M]
lin=new _line[M];
map=new _cell*[N];
for (i=0;i<N;i++) map[i]=new _cell[N];
// reset map[N][N]
for (i=0;i<N;i++)
for (j=0;j<N;j++)
map[i][j].cnt=0;
// compute point cloud density
for (k=0;k<view.pnt.num;k++)
{
p=&view.pnt[k];
i=double((p->p[0]-x0)*mx); if (i<0) i=0; if (i>=N) i=N-1;
j=double((p->p[1]-y0)*my); if (j<0) j=0; if (j>=N) j=N-1;
m=&map[i][j];
if (!m->cnt)
{
m->x0=p->p[0];
m->x1=p->p[0];
m->y0=p->p[1];
m->y1=p->p[1];
}
if (m->cnt<0x7FFFFFFF) m->cnt++; // avoid overflow
if (m->x0>p->p[0]) m->x0=p->p[0];
if (m->x1<p->p[0]) m->x1=p->p[0];
if (m->y0>p->p[1]) m->y0=p->p[1];
if (m->y1<p->p[1]) m->y1=p->p[1];
}
// find holes (map[i][j].cnt==0) or (map[i][j].cnt<=treshold)
// and create lin[] list of H,V lines covering holes
for (j=0;j<N;j++) // search lines
{
for (i=0;i<N;)
{
int i0,i1;
for (;i<N;i++) if (map[i][j].cnt==0) break; i0=i-1; // find start of hole
for (;i<N;i++) if (map[i][j].cnt!=0) break; i1=i; // find end of hole
if (i0< 0) continue; // skip bad circumstances (edges or no hole found)
if (i1>=N) continue;
if (map[i0][j].cnt==0) continue;
if (map[i1][j].cnt==0) continue;
_line l;
l.i0=i0; l.x0=map[i0][j].x1;
l.i1=i1; l.x1=map[i1][j].x0;
l.j0=j ; l.y0=0.25*(map[i0][j].y0+map[i0][j].y1+map[i1][j].y0+map[i1][j].y1);
l.j1=j ; l.y1=l.y0;
lin[lins]=l; lins++;
}
}
for (i=0;i<N;i++) // search columns
{
for (j=0;j<N;)
{
int j0,j1;
for (;j<N;j++) if (map[i][j].cnt==0) break; j0=j-1; // find start of hole
for (;j<N;j++) if (map[i][j].cnt!=0) break; j1=j; // find end of hole
if (j0< 0) continue; // skip bad circumstances (edges or no hole found)
if (j1>=N) continue;
if (map[i][j0].cnt==0) continue;
if (map[i][j1].cnt==0) continue;
_line l;
l.i0=i ; l.y0=map[i][j0].y1;
l.i1=i ; l.y1=map[i][j1].y0;
l.j0=j0; l.x0=0.25*(map[i][j0].x0+map[i][j0].x1+map[i][j1].x0+map[i][j1].x1);
l.j1=j1; l.x1=l.x0;
lin[lins]=l; lins++;
}
}
// segmentate lin[] ... group lines of the same hole together by lin[].id
// segmentation based on vector lines data
// you can also segmentate the map[][] directly as bitmap during hole detection
for (i=0;i<lins;i++) lin[i].id=i; // all lines are separate
for (;;) // join what you can
{
int e=0,i0,i1;
_line *a,*b;
for (a=lin,i=0;i<lins;i++,a++)
{
for (b=a,j=i;j<lins;j++,b++)
if (a->id!=b->id)
{
// do 2D lines a,b intersect ?
double xx0,yy0,xx1,yy1;
double kx0,ky0,dx0,dy0,t0;
double kx1,ky1,dx1,dy1,t1;
double x0=a->x0,y0=a->y0;
double x1=a->x1,y1=a->y1;
double x2=b->x0,y2=b->y0;
double x3=b->x1,y3=b->y1;
// discart lines with non intersecting bound rectangles
double a0,a1,b0,b1;
if (x0<x1) { a0=x0; a1=x1; } else { a0=x1; a1=x0; }
if (x2<x3) { b0=x2; b1=x3; } else { b0=x3; b1=x2; }
if (a1<b0) continue;
if (a0>b1) continue;
if (y0<y1) { a0=y0; a1=y1; } else { a0=y1; a1=y0; }
if (y2<y3) { b0=y2; b1=y3; } else { b0=y3; b1=y2; }
if (a1<b0) continue;
if (a0>b1) continue;
// compute intersection
kx0=x0; ky0=y0; dx0=x1-x0; dy0=y1-y0;
kx1=x2; ky1=y2; dx1=x3-x2; dy1=y3-y2;
t1=divide(dx0*(ky0-ky1)+dy0*(kx1-kx0),(dx0*dy1)-(dx1*dy0));
xx1=kx1+(dx1*t1);
yy1=ky1+(dy1*t1);
if (fabs(dx0)>=fabs(dy0)) t0=divide(kx1-kx0+(dx1*t1),dx0);
else t0=divide(ky1-ky0+(dy1*t1),dy0);
xx0=kx0+(dx0*t0);
yy0=ky0+(dy0*t0);
// check if intersection exists
if (fabs(xx1-xx0)>1e-6) continue;
if (fabs(yy1-yy0)>1e-6) continue;
if ((t0<0.0)||(t0>1.0)) continue;
if ((t1<0.0)||(t1>1.0)) continue;
// if yes ... intersection point = xx0,yy0
e=1; break;
}
if (e) break; // join found ... stop searching
}
if (!e) break; // no join found ... stop segmentation
i0=a->id; // joid ids ... rename i1 to i0
i1=b->id;
for (a=lin,i=0;i<lins;i++,a++)
if (a->id==i1)
a->id=i0;
}
// visualize lin[]
for (i=0;i<lins;i++)
{
glview2D::_lin l;
l.p0.p[0]=lin[i].x0;
l.p0.p[1]=lin[i].y0;
l.p1.p[0]=lin[i].x1;
l.p1.p[1]=lin[i].y1;
// l.col=0x0000FF00;
l.col=(lin[i].id*0x00D00C10A)+0x00800000; // color is any function of ID
view.lin.add(l);
}
// dynamic deallocation of map[N][N],lin[M]
for (i=0;i<N;i++) delete[] map[i];
delete[] map;
delete[] lin;
}
//---------------------------------------------------------------------------
忽略我的 glview2D 东西(它是我的几何图形渲染引擎)
Just ignore my glview2D stuff (it is my gfx render engine for geometry)
view.pnt[]是点的动态列表(随机生成)view.lin[]是动态列表输出H,V 线，仅用于可视化lin[]是你的行输出
view.pnt[]is dynamic list of your points (generated by random)view.lin[]is dynamic list output H,V lines for visualization onlylin[]is your lines output
这是输出:
我懒得添加多边形化了，现在你可以看到分割有效(着色).如果您还需要多边形化的帮助，请评论我，但我认为这应该没有任何问题.
I am too lazy to add polygonize for now you can see that segmentation works (coloring). If you need also help with polygonize then comment me but I think that should not be any problem.
复杂度估计取决于整体空洞覆盖率
Complexity estimation depends on the overall hole coverage
但是对于大部分代码来说，它是 O(N) 并且对于孔搜索/分割 ~O((M^2)+(U^2))其中:
but for most of the code it is O(N) and for hole search/segmentation ~O((M^2)+(U^2)) where:
N是点数M是地图网格大小U是 H,V 线 依赖于孔的计数...M <<
Nis point countMis map grid sizeUis H,V lines count dependent on holes ...M << N, U << M*M
正如您在上图中看到的 3783 点 30x30 网格，我的设置花费了几乎 9ms
as you can see for 3783 points 30x30 grid on the image above it took almost 9ms on my setup
[Edit1] 稍微玩弄矢量多边形
对于简单的孔很好，但对于更复杂的孔还有一些问题
for simple holes is fine but for more complicated ones there are some hick-ups yet
[Edit2] 终于有时间做这个了，就是这样:
这是一个以更愉快/易于管理的形式进行孔/多边形搜索的简单类:
This is simple class for hole/polygon search in more pleasant/manageable form:
//---------------------------------------------------------------------------
class holes
{
public:
int xs,ys,n; // cell grid x,y - size and points count
int **map; // points density map[xs][ys]
// i=(x-x0)*g2l; x=x0+(i*l2g);
// j=(y-y0)*g2l; y=y0+(j*l2g);
double mg2l,ml2g; // scale to/from global/map space (x,y) <-> map[i][j]
double x0,x1,y0,y1; // used area (bounding box)
struct _line
{
int id; // id of hole for segmentation/polygonize
int i0,i1,j0,j1; // index in map[][]
_line(){}; _line(_line& a){ *this=a; }; ~_line(){}; _line* operator = (const _line *a) { *this=*a; return this; }; /*_line* operator = (const _line &a) { ...copy... return this; };*/
};
List<_line> lin;
int lin_i0; // start index for perimeter lines (smaller indexes are the H,V lines inside hole)
struct _point
{
int i,j; // index in map[][]
int p0,p1; // previous next point
int used;
_point(){}; _point(_point& a){ *this=a; }; ~_point(){}; _point* operator = (const _point *a) { *this=*a; return this; }; /*_point* operator = (const _point &a) { ...copy... return this; };*/
};
List<_point> pnt;
// class init and internal stuff
holes() { xs=0; ys=0; n=0; map=NULL; mg2l=1.0; ml2g=1.0; x0=0.0; y0=0.0; x1=0.0; y1=0.0; lin_i0=0; };
holes(holes& a){ *this=a; };
~holes() { _free(); };
holes* operator = (const holes *a) { *this=*a; return this; };
holes* operator = (const holes &a)
{
xs=0; ys=0; n=a.n; map=NULL;
mg2l=a.mg2l; x0=a.x0; x1=a.x1;
ml2g=a.ml2g; y0=a.y0; y1=a.y1;
_alloc(a.xs,a.ys);
for (int i=0;i<xs;i++)
for (int j=0;j<ys;j++) map[i][j]=a.map[i][j];
return this;
}
void _free() { if (map) { for (int i=0;i<xs;i++) if (map[i]) delete[] map[i]; delete[] map; } xs=0; ys=0; }
void _alloc(int _xs,int _ys) { int i=0; _free(); xs=_xs; ys=_ys; map=new int*[xs]; if (map) for (i=0;i<xs;i++) { map[i]=new int[ys]; if (map[i]==NULL) { i=-1; break; } } else i=-1; if (i<0) _free(); }
// scann boundary box interface
void scann_beg();
void scann_pnt(double x,double y);
void scann_end();
// dynamic allocations
void cell_size(double sz); // compute/allocate grid from grid cell size = sz x sz
// scann holes interface
void holes_beg();
void holes_pnt(double x,double y);
void holes_end();
// global(x,y) <- local map[i][j] + half cell offset
inline void l2g(double &x,double &y,int i,int j) { x=x0+((double(i)+0.5)*ml2g); y=y0+((double(j)+0.5)*ml2g); }
// local map[i][j] <- global(x,y)
inline void g2l(int &i,int &j,double x,double y) { i= double((x-x0) *mg2l); j= double((y-y0) *mg2l); }
};
//---------------------------------------------------------------------------
void holes::scann_beg()
{
x0=0.0; y0=0.0; x1=0.0; y1=0.0; n=0;
}
//---------------------------------------------------------------------------
void holes::scann_pnt(double x,double y)
{
if (!n) { x0=x; y0=y; x1=x; y1=y; }
if (n<0x7FFFFFFF) n++; // avoid overflow
if (x0>x) x0=x; if (x1<x) x1=x;
if (y0>y) y0=y; if (y1<y) y1=y;
}
//---------------------------------------------------------------------------
void holes::scann_end()
{
}
//---------------------------------------------------------------------------
void holes::cell_size(double sz)
{
int x,y;
if (sz<1e-6) sz=1e-6;
x=ceil((x1-x0)/sz);
y=ceil((y1-y0)/sz);
_alloc(x,y);
ml2g=sz; mg2l=1.0/sz;
}
//---------------------------------------------------------------------------
void holes::holes_beg()
{
int i,j;
for (i=0;i<xs;i++)
for (j=0;j<ys;j++)
map[i][j]=0;
}
//---------------------------------------------------------------------------
void holes::holes_pnt(double x,double y)
{
int i,j;
g2l(i,j,x,y);
if ((i>=0)&&(i<xs))
if ((j>=0)&&(j<ys))
if (map[i][j]<0x7FFFFFFF) map[i][j]++; // avoid overflow
}
//---------------------------------------------------------------------------
void holes::holes_end()
{
int i,j,e,i0,i1;
List<int> ix; // hole lines start/stop indexes for speed up the polygonization
_line *a,*b,l;
_point *aa,*bb,p;
lin.num=0; lin_i0=0;// clear lines
ix.num=0; // clear indexes
// find holes (map[i][j].cnt==0) or (map[i][j].cnt<=treshold)
// and create lin[] list of H,V lines covering holes
for (j=0;j<ys;j++) // search lines
for (i=0;i<xs;)
{
int i0,i1;
for (;i<xs;i++) if (map[i][j]==0) break; i0=i-1; // find start of hole
for (;i<xs;i++) if (map[i][j]!=0) break; i1=i; // find end of hole
if (i0< 0) continue; // skip bad circumstances (edges or no hole found)
if (i1>=xs) continue;
if (map[i0][j]==0) continue;
if (map[i1][j]==0) continue;
l.i0=i0;
l.i1=i1;
l.j0=j ;
l.j1=j ;
l.id=-1;
lin.add(l);
}
for (i=0;i<xs;i++) // search columns
for (j=0;j<ys;)
{
int j0,j1;
for (;j<ys;j++) if (map[i][j]==0) break; j0=j-1; // find start of hole
for (;j<ys;j++) if (map[i][j]!=0) break; j1=j ; // find end of hole
if (j0< 0) continue; // skip bad circumstances (edges or no hole found)
if (j1>=ys) continue;
if (map[i][j0]==0) continue;
if (map[i][j1]==0) continue;
l.i0=i ;
l.i1=i ;
l.j0=j0;
l.j1=j1;
l.id=-1;
lin.add(l);
}
// segmentate lin[] ... group lines of the same hole together by lin[].id
// segmentation based on vector lines data
// you can also segmentate the map[][] directly as bitmap during hole detection
for (i=0;i<lin.num;i++) lin[i].id=i; // all lines are separate
for (;;) // join what you can
{
for (e=0,a=lin.dat,i=0;i<lin.num;i++,a++)
{
for (b=a,j=i;j<lin.num;j++,b++)
if (a->id!=b->id)
{
// if a,b not intersecting or neighbouring
if (a->i0>b->i1) continue;
if (b->i0>a->i1) continue;
if (a->j0>b->j1) continue;
if (b->j0>a->j1) continue;
// if they do mark e for join groups
e=1; break;
}
if (e) break; // join found ... stop searching
}
if (!e) break; // no join found ... stop segmentation
i0=a->id; // joid ids ... rename i1 to i0
i1=b->id;
for (a=lin.dat,i=0;i<lin.num;i++,a++)
if (a->id==i1)
a->id=i0;
}
// sort lin[] by id
for (e=1;e;) for (e=0,a=&lin[0],b=&lin[1],i=1;i<lin.num;i++,a++,b++)
if (a->id>b->id) { l=*a; *a=*b; *b=l; e=1; }
// re id lin[] and prepare start/stop indexes
for (i0=-1,i1=-1,a=&lin[0],i=0;i<lin.num;i++,a++)
if (a->id==i1) a->id=i0;
else { i0++; i1=a->id; a->id=i0; ix.add(i); }
ix.add(lin.num);
// polygonize
lin_i0=lin.num;
for (j=1;j<ix.num;j++) // process hole
{
i0=ix[j-1]; i1=ix[j];
// create border pnt[] list (unique points only)
pnt.num=0; p.used=0; p.p0=-1; p.p1=-1;
for (a=&lin[i0],i=i0;i<i1;i++,a++)
{
p.i=a->i0;
p.j=a->j0;
map[p.i][p.j]=0;
for (aa=&pnt[0],e=0;e<pnt.num;e++,aa++)
if ((aa->i==p.i)&&(aa->j==p.j)) { e=-1; break; }
if (e>=0) pnt.add(p);
p.i=a->i1;
p.j=a->j1;
map[p.i][p.j]=0;
for (aa=&pnt[0],e=0;e<pnt.num;e++,aa++)
if ((aa->i==p.i)&&(aa->j==p.j)) { e=-1; break; }
if (e>=0) pnt.add(p);
}
// mark not border points
for (aa=&pnt[0],i=0;i<pnt.num;i++,aa++)
if (!aa->used) // ignore marked points
if ((aa->i>0)&&(aa->i<xs-1)) // ignore map[][] border points
if ((aa->j>0)&&(aa->j<ys-1))
{ // ignore if any non hole cell around
if (map[aa->i-1][aa->j-1]>0) continue;
if (map[aa->i-1][aa->j ]>0) continue;
if (map[aa->i-1][aa->j+1]>0) continue;
if (map[aa->i ][aa->j-1]>0) continue;
if (map[aa->i ][aa->j+1]>0) continue;
if (map[aa->i+1][aa->j-1]>0) continue;
if (map[aa->i+1][aa->j ]>0) continue;
if (map[aa->i+1][aa->j+1]>0) continue;
aa->used=1;
}
// delete marked points
for (aa=&pnt[0],e=0,i=0;i<pnt.num;i++,aa++)
if (!aa->used) { pnt[e]=*aa; e++; } pnt.num=e;
// connect neighbouring points distance=1
for (i0= 0,aa=&pnt[i0];i0<pnt.num;i0++,aa++)
if (aa->used<2)
for (i1=i0+1,bb=&pnt[i1];i1<pnt.num;i1++,bb++)
if (bb->used<2)
{
i=aa->i-bb->i; if (i<0) i=-i; e =i;
i=aa->j-bb->j; if (i<0) i=-i; e+=i;
if (e!=1) continue;
aa->used++; if (aa->p0<0) aa->p0=i1; else aa->p1=i1;
bb->used++; if (bb->p0<0) bb->p0=i0; else bb->p1=i0;
}
// try to connect neighbouring points distance=sqrt(2)
for (i0= 0,aa=&pnt[i0];i0<pnt.num;i0++,aa++)
if (aa->used<2)
for (i1=i0+1,bb=&pnt[i1];i1<pnt.num;i1++,bb++)
if (bb->used<2)
if ((aa->p0!=i1)&&(aa->p1!=i1))
if ((bb->p0!=i0)&&(bb->p1!=i0))
{
if ((aa->used)&&(aa->p0==bb->p0)) continue; // avoid small closed loops
i=aa->i-bb->i; if (i<0) i=-i; e =i*i;
i=aa->j-bb->j; if (i<0) i=-i; e+=i*i;
if (e!=2) continue;
aa->used++; if (aa->p0<0) aa->p0=i1; else aa->p1=i1;
bb->used++; if (bb->p0<0) bb->p0=i0; else bb->p1=i0;
}
// try to connect to closest point
int ii,dd;
for (i0= 0,aa=&pnt[i0];i0<pnt.num;i0++,aa++)
if (aa->used<2)
{
for (ii=-1,i1=i0+1,bb=&pnt[i1];i1<pnt.num;i1++,bb++)
if (bb->used<2)
if ((aa->p0!=i1)&&(aa->p1!=i1))
if ((bb->p0!=i0)&&(bb->p1!=i0))
{
i=aa->i-bb->i; if (i<0) i=-i; e =i*i;
i=aa->j-bb->j; if (i<0) i=-i; e+=i*i;
if ((ii<0)||(e<dd)) { ii=i1; dd=e; }
}
if (ii<0) continue;
i1=ii; bb=&pnt[i1];
aa->used++; if (aa->p0<0) aa->p0=i1; else aa->p1=i1;
bb->used++; if (bb->p0<0) bb->p0=i0; else bb->p1=i0;
}
// add connected points to lin[] ... this is hole perimeter !!!
// lines are 2 x duplicated so some additional code for sort the order of line swill be good idea
l.id=lin[ix[j-1]].id;
for (i0=0,aa=&pnt[i0];i0<pnt.num;i0++,aa++)
{
l.i0=aa->i;
l.j0=aa->j;
// [edit3] this avoid duplicating lines
if (aa->p0>i0) { bb=&pnt[aa->p0]; l.i1=bb->i; l.j1=bb->j; lin.add(l); }
if (aa->p1>i0) { bb=&pnt[aa->p1]; l.i1=bb->i; l.j1=bb->j; lin.add(l); }
//if (aa->p0>=0) { bb=&pnt[aa->p0]; l.i1=bb->i; l.j1=bb->j; lin.add(l); }
//if (aa->p1>=0) { bb=&pnt[aa->p1]; l.i1=bb->i; l.j1=bb->j; lin.add(l); }
}
}
}
//---------------------------------------------------------------------------
您只需要将我的 List 模板替换为 std::list 或其他(我无法共享的模板).它是 T ...
You just need to replace my List<T> template with std::list or whatever (that template I cannot share). It is an dynamic 1D array of T ...
List与x; int x[]; 相同x.add();向 x 添加空项x.add(a);向 x 添加一项x.reset()清空数组x.allocate(size)预分配空间以避免在运行时重新分配很慢x.num是 x[] 中的项目数 ... 项目中使用的大小
List<int> x;is the same asint x[];x.add();add empty item to xx.add(a);add a item to xx.reset()clears the arrayx.allocate(size)preallocate space to avoid reallocations on the run which is slowx.numis number of items in x[] ... used size in items
在原始代码中只是静态数组，所以如果您感到困惑，请改为使用它.
in the original code are only static arrays so if you are confused check with it instead.
现在如何使用它:
h.scann_beg(); for (i=0;i<view.pnt.num;i++) { p=view.pnt[i].p0.p; h.scann_pnt(p[0],p[1]); } h.scann_end();
h.cell_size(2.5);
h.holes_beg(); for (i=0;i<view.pnt.num;i++) { p=view.pnt[i].p0.p; h.holes_pnt(p[0],p[1]); } h.holes_end();
其中 view.pnt[] 是输入点列表及其内部: view.pnt[i].p0.p[ 2 ]= { x,y }
where view.pnt[] is list of input points and inside it: view.pnt[i].p0.p[ 2 ]= { x,y }
输出在 h.lin[] 和 lin_i0 中，其中:
Output is in h.lin[] and lin_i0 where:
h.lin[i] i= <0,lin_i0 )是里面的 H,V 线h.lin[i] i= <lin_i0,h.lin.num )是周长
h.lin[i] i= < 0,lin_i0 )are the inside H,V linesh.lin[i] i= < lin_i0,h.lin.num )are the perimeter
周边线没有排序并且重复两次，所以只需对它们进行排序并删除重复项(太懒了).lin[] 里面是 id .. id 线所属的孔 0,1,2,3,... 和 i,j 地图内的坐标.所以为了正确输出到您的世界坐标中，请执行以下操作:
The perimeter lines are not ordered and are duplicated twice so just order them and remove duplicates (too lazy for that). Inside lin[] are id .. id of hole 0,1,2,3,... to which the line belongs and i,j coordinates inside map. so for proper output into your world coordinates do something like this:
int i,j;
holes h; // holes class
double *p; // input point list ptr
h.scann_beg(); for (i=0;i<view.pnt.num;i++) { p=view.pnt[i].p0.p; h.scann_pnt(p[0],p[1]); } h.scann_end();
h.cell_size(2.5);
h.holes_beg(); for (i=0;i<view.pnt.num;i++) { p=view.pnt[i].p0.p; h.holes_pnt(p[0],p[1]); } h.holes_end();
DWORD coltab[]=
{
0x000000FF,
0x0000FF00,
0x00FF0000,
0x0000FFFF,
0x00FFFF00,
0x00FF00FF,
0x00FFFFFF,
0x00000088,
0x00008800,
0x00880000,
0x00008888,
0x00888800,
0x00880088,
0x00888888,
};
for (i=0;i<h.lin.num;i++) // draw lin[]
{
glview2D::_lin a;
holes::_line *b=&h.lin[i];
h.l2g(a.p0.p[0],a.p0.p[1],b->i0,b->j0);
h.l2g(a.p1.p[0],a.p1.p[1],b->i1,b->j1);
if (i<h.lin_i0) // H,V lines inside hole(b->id) .. gray [edit3] was <= which is wrong and miss-color first perimeter line
{
a.col=0x00808080;
}
else{ // hole(b->id) perimeter lines ... each hole different collor
if ((b->id>=0)&&(b->id<14)) a.col=coltab[b->id];
if (b->id==-1) a.col=0x00FFFFFF; // special debug lines
if (b->id==-2) a.col=0x00AA8040; // special debug lines
}
view.lin.add(a); // here draw your line or add it to your polygon instead
}
- 我的
view.lin[]有成员:p0,p1,，它们是view.pnt[]和col即颜色 - my
view.lin[]has members:p0,p1,which are points asview.pnt[]andcolwhich is color
当孔太小时，我只看到一个问题(直径 <3 个单元) 否则没问题
I saw only one issue with this when holes are too small (diameter < 3 cells) otherwise is OK
[edit4] 重新排列边界线
这样做而不是这样:
/* add connected points to lin[] ... this is hole perimeter !!!
// lines are 2 x duplicated so some additional code for sort the order of line swill be good idea
l.id=lin[ix[j-1]].id;
for (i0=0,aa=&pnt[i0];i0<pnt.num;i0++,aa++)
{
l.i0=aa->i;
l.j0=aa->j;
// [edit3] this avoid duplicating lines
if (aa->p0>i0) { bb=&pnt[aa->p0]; l.i1=bb->i; l.j1=bb->j; lin.add(l); }
if (aa->p1>i0) { bb=&pnt[aa->p1]; l.i1=bb->i; l.j1=bb->j; lin.add(l); }
//if (aa->p0>=0) { bb=&pnt[aa->p0]; l.i1=bb->i; l.j1=bb->j; lin.add(l); }
//if (aa->p1>=0) { bb=&pnt[aa->p1]; l.i1=bb->i; l.j1=bb->j; lin.add(l); }
} */
这样做:
// add connected points to lin[] ... this is hole perimeter !!!
l.id=lin[ix[j-1]].id;
// add index of points instead points
int lin_i1=lin.num;
for (i0=0,aa=&pnt[i0];i0<pnt.num;i0++,aa++)
{
l.i0=i0;
if (aa->p0>i0) { l.i1=aa->p0; lin.add(l); }
if (aa->p1>i0) { l.i1=aa->p1; lin.add(l); }
}
// reorder perimeter lines
for (i0=lin_i1,a=&lin[i0];i0<lin.num-1;i0++,a++)
for (i1=i0+1 ,b=&lin[i1];i1<lin.num ;i1++,b++)
{
if (a->i1==b->i0) { a++; l=*a; *a=*b; *b=l; a--; break; }
if (a->i1==b->i1) { a++; l=*a; *a=*b; *b=l; i=a->i0; a->i0=a->i1; a->i1=i; a--; break; }
}
// convert point indexes to points
for (i0=lin_i1,a=&lin[i0];i0<lin.num;i0++,a++)
{
bb=&pnt[a->i0]; a->i0=bb->i; a->j0=bb->j;
bb=&pnt[a->i1]; a->i1=bb->i; a->j1=bb->j;
}
[Edit5] holes::holes_end 内部的多边形化是如何工作的
How polygonize inside holes::holes_end works
作为输入，您需要所有 H,V 线 lin[] 按孔分段/分组/排序的列表和密度图 map[][].
As input for this you need the list of all H,V lines lin[] segmentated/grouped/sorted by hole and the density map map[][].
遍历所有孔
遍历加工孔的所有H、V线
创建所有唯一线路端点的列表pnt[](无重复).因此，为每条线取 2 个端点，然后查看每个点是否已在列表中.如果不添加它，则忽略它.
Create list of all unique line endpoints pnt[] (no duplicates). So take 2 endpoints for each line and look if each point is already in the list. If not add it there else ignore it.
从列表中删除所有非边界点
因此通过查看密度中的 4 个相邻点来删除所有与非孔区域没有接触的点 map[][]
So remove all points that have no contact with non-hole area by looking into 4 neighbors in the density map[][]
对点进行连通分量分析
- set
used=0;p0=-1;p1=-1;对于pnt[]列表中的所有点 用
distance=1
- set
used=0; p0=-1; p1=-1;for all points inpnt[]list connect points with
distance=1
使用 used<2 遍历所有点 pnt[] 这意味着它们尚未完全使用，并且对于每个这样的点搜索 pnt[] 再次为另一个这样的点也有 distance = 1 到它.这意味着它是它的 4 个邻居并且应该被连接，所以添加连接信息到它们的摊位(使用 p0 或 p1 索引，无论未使用 (-1)) 并增加两个点的使用量.
loop through all points pnt[] with used<2 which means they are not fully used yet and for each such point search pnt[] again for another such point that has also distance = 1 to it. It means it is its 4-neighbors and should be connected so add the connection info to booth of them (use p0 or p1 index which ever is unused (-1)) and increment usage of both points.
尝试用distance=sqrt(2)
几乎与 #2 相同，除了现在选择 8 个邻居的对角线的距离.这次也避免了闭环，所以不要连接已经连接到它的点.
is almost the same as #2 except the distance which now selects diagonals of 8-neighbors. This time also avoid closed loops so do not connect point that is already connected to it.
尝试连接最近的点
再次与 #2,#3 几乎相同，但选择最近的点，并避免闭环.
again is almost the same as #2,#3 but select the closest point instead and also avoid closed loops.
从pnt[]
所以选择列表中的第一个点并将其添加到多边形中.然后将连接点添加到它(与您以哪种方式启动 p0 或 p1 无关).然后添加它的连接点(不同于之前添加到多边形的点以避免前后循环).添加与 pnt[] 中的点一样多的点.
so pick first point in list and add it to the polygon. then add the connected point to it (does not matter which way you start p0 or p1). Then add its connected point (different then previous added point to polygon to avoid back and forward loops). Add so many points as you have points in a pnt[].
这篇关于在二维点集中寻找洞?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
