¿Encontrar agujeros en conjuntos de 2d puntos?

Tengo un conjunto de 2d points . Son X,Y coordinates en un sistema de grilla cartesiana estándar (en este caso, una UTM zone ). Necesito encontrar los agujeros en ese punto establecidos preferiblemente con alguna habilidad para establecer la sensibilidad del algoritmo que encuentra estos agujeros. Normalmente, estos conjuntos de puntos son muy densos, pero algunos pueden ser mucho menos densos.

Lo que son, si es que lo ayuda, puntos en los que se ha muestreado el suelo en el campo para varias propiedades que las personas en la agricultura aparentemente encuentran útiles. Algunas veces en estas muestras puntuales hay rocas gigantes o lugares pantanosos o llenos de pequeños lagos y estanques.

Desde el punto establecido, quieren que encuentre el polígono cóncavo que define estos agujeros aproximadamente.

Ya escribí el algoritmo que encuentra el polígono de límite cóncavo exterior y les permite establecer la sensibilidad para qué tan áspero o liso sea el polígono formado por el algoritmo. Después de eso, les gustaría encontrar los agujeros y tener esos agujeros como un polígono cóncavo, que supongo que en algunos casos puede ser convexo, pero ese será el borde de la caja, no la norma.

El problema es que los únicos documentos que he escuchado sobre este tema son los realizados por astrónomos que quieren encontrar grandes bolsas de vacío en el espacio y nunca he podido encontrar uno de sus artículos con un algoritmo real mostrado en ellos como cualquier cosa que no sea una descripción conceptual aproximada.

He intentado con Google y varias búsquedas de documentos académicos, etc. pero no he encontrado mucho que sea útil hasta ahora. Lo que me hace preguntarme si tal vez hay un nombre para este tipo de problema y no lo sé, ¿entonces estoy buscando algo equivocado o algo así?

Entonces, después de esa larga explicación, mi pregunta es: ¿Alguien sabe lo que debería buscar para encontrar artículos sobre esto preferiblemente con algoritmos bien definidos o alguien sabe un algoritmo que hará esto para que me lo indiquen?

Cualquier cosa que me ayude a resolver este problema sería muy útil y muy apreciada, incluso las frases de búsqueda correctas que harán aparecer los posibles algoritmos o documentos.

El lenguaje en el que se implementará será C # pero los ejemplos en cualquier cosa, desde paquetes de Mathematica hasta MATLAB or ASM, C, C++, Python, Java or MathCAD , estarán bien siempre y cuando el ejemplo no tenga algunas llamadas que vaya a cosas como FindTheHole etc. Donde FindTheHole no está definido o es propiedad del software de implementación, por ejemplo, los ejemplos de MathCAD suelen tener mucho de eso.

A continuación hay dos ejemplos de conjuntos de puntos reales, uno denso y uno disperso y las áreas en ellos que necesitaría encontrar: Puntos dispersos, por ejemploPuntos densos, por ejemplo

¿Qué pasa con un enfoque vectorial de bitmap como este?

  1. obtener un cuadro delimitador de la cobertura del área de nubes de puntos

    Haz esto si aún no se sabe. Debe ser simple O(N) ciclo a través de todos los puntos.

  2. crear el map[N][N] del área

    Es un “bitmap” del área para un fácil cálculo de la densidad de datos. Simplemente crea una proyección desde el area(x,y) -> map[i][j] por ejemplo con escala simple. El tamaño de cuadrícula N también es la precisión de la salida y debe ser mayor que la distancia promedio del punto. de modo que cada celda dentro del map[][] cubre el área con al menos un punto (si no está en el área del agujero).

  3. calcular la densidad de datos para cada celda del map[][]

    Fácil como un pastel, simplemente borre el map[][].cnt (contador de puntos) a zero y calcule con un simple ciclo O(N) donde haga el map[i][j].cnt++ para todos los points(x,y)

  4. crear una lista de área no utilizada (map[][].cnt==0) o (map[][].cnt<=treshold)

    Lo hago por líneas horizontales y verticales por simplicidad

  5. salida segmentada

    Solo agrupe líneas del mismo agujero juntas (intersecando ... enfoque vectorial) y también se puede hacer en la viñeta # 4 por relleno de inundación (enfoque de bitmap)

  6. salida poligonalizada

    Tome todos los puntos de borde de las líneas H, V del mismo agujero / grupo y cree un polígono (ordénelos para que su conexión no intersecte con nada). Hay muchas librerías, algoritmos y código fuente sobre esto.

Mi código fuente para este enfoque:

 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;ip->p[0]) x0=p->p[0]; if (x1p[0]) x1=p->p[0]; if (y0>p->p[1]) y0=p->p[1]; if (y1p[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;ip[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->x1p[0]) m->x1=p->p[0]; if (m->y0>p->p[1]) m->y0=p->p[1]; if (m->y1p[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) 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) 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;iid!=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 (x0b1) continue; if (y0b1) 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;iid==i1) a->id=i0; } // visualize lin[] for (i=0;i 

Solo ignore mis cosas glview2D (es mi motor de renderizado gfx para geometría)

  • view.pnt[] es una lista dinámica de tus puntos (generada por azar)
  • view.lin[] es una salida de lista dinámica H, líneas V para visualización solamente
  • lin[] es la salida de tus líneas

Esto es resultado:

vista previa de agujeros

Soy demasiado perezoso para agregar Polygonize por ahora puedes ver que la segmentación funciona (colorear). Si también necesita ayuda con polygonize, coménteme, pero creo que no debería haber ningún problema.

La estimación de la complejidad depende de la cobertura total del pozo

pero para la mayoría del código es O(N) y para la búsqueda / segmentación de agujeros ~O((M^2)+(U^2)) donde:

  • N es el conteo de puntos
  • M es el tamaño de cuadrícula del mapa
  • U es H, las líneas V cuentan en función de los agujeros ...
  • M << N, U << M*M

como se puede ver en 3783 puntos, la cuadrícula de 30x30 en la imagen de arriba tomó casi 9ms en mi configuración

[Edit1] jugado con vector polígono un poco

agujeros confinados

para los agujeros simples está bien, pero para los más complicados hay algunos hick-ups todavía

[Edit2] finalmente tiene un poco de tiempo para esto, así que aquí está:

Esta es una clase simple para la búsqueda de agujeros / polígonos en forma más agradable / manejable:

 //--------------------------------------------------------------------------- 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=an; 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;ix) x0=x; if (x1y) y0=y; if (y1=0)&&(i=0)&&(j 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=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=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;iid!=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;iid==i1) a->id=i0; } // sort lin[] by id for (e=1;e;) for (e=0,a=&lin[0],b=&lin[1],i=1;iid>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;iid==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;ji0; pj=a->j0; map[pi][pj]=0; for (aa=&pnt[0],e=0;ei==pi)&&(aa->j==pj)) { e=-1; break; } if (e>=0) pnt.add(p); pi=a->i1; pj=a->j1; map[pi][pj]=0; for (aa=&pnt[0],e=0;ei==pi)&&(aa->j==pj)) { e=-1; break; } if (e>=0) pnt.add(p); } // mark not border points for (aa=&pnt[0],i=0;iused) // ignore marked points if ((aa->i>0)&&(aa->ij>0)&&(aa->ji-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;iused) { pnt[e]=*aa; e++; } pnt.num=e; // connect neighbouring points distance=1 for (i0= 0,aa=&pnt[i0];i0used<2) for (i1=i0+1,bb=&pnt[i1];i1used<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];i0used<2) for (i1=i0+1,bb=&pnt[i1];i1used<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];i0used<2) { for (ii=-1,i1=i0+1,bb=&pnt[i1];i1used<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)||(eused++; 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];i0i; 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); } } } } //--------------------------------------------------------------------------- 

Solo necesita reemplazar mi plantilla List con std::list o lo que sea (esa plantilla que no puedo compartir). Es una matriz 1D dinámica de T ...

  • List x; es lo mismo que int x[];
  • x.add(); añadir un elemento vacío a x
  • x.add(a); agregar un elemento a x
  • x.reset() borra la matriz
  • x.allocate(size) preasignar espacio para evitar reasignaciones en la ejecución que es lenta
  • x.num es el número de elementos en x [] ... tamaño utilizado en los elementos

en el código original solo hay matrices estáticas, por lo que si está confundido, consulte con él.

Ahora cómo usarlo:

 h.scann_beg(); for (i=0;i 

donde view.pnt[] es una lista de puntos de entrada y dentro de él: view.pnt[i].p0.p[ 2 ]= { x,y }

La salida está en h.lin[] y lin_i0 donde:

  • h.lin[i] i= < 0,lin_i0 ) son las líneas interiores H, V
  • h.lin[i] i= < lin_i0,h.lin.num ) son el perímetro

Las líneas del perímetro no están ordenadas y se duplican dos veces, así que simplemente pídalas y elimine duplicados (demasiado flojo para eso). Dentro de lin[] encuentran id .. id del agujero 0,1,2,3,... al que pertenece la línea i,j coordenadas dentro del mapa. entonces, para una salida correcta en las coordenadas de tu mundo, haz algo como esto:

 int i,j; holes h; // holes class double *p; // input point list ptr h.scann_beg(); for (i=0;ii0,b->j0); h.l2g(a.p1.p[0],a.p1.p[1],b->i1,b->j1); if (iid) .. 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 } 
  • mi view.lin[] tiene miembros: p0,p1, que son puntos como view.pnt[] y col que es color

Solo vi un problema cuando los agujeros son demasiado pequeños (diameter < 3 cells) contrario, está bien

[edit4] reordenando líneas de perímetro

para hacer eso solo en vez de esto:

  /* 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];i0i; 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); } } */ 

hacer esto:

  // 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];i0p0>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];i0i1==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];i0i0]; a->i0=bb->i; a->j0=bb->j; bb=&pnt[a->i1]; a->i1=bb->i; a->j1=bb->j; } 

[Edit5] Cómo holes::holes_end funciona

Como entrada para esto necesita la lista de todas las líneas H, V lin[] segmentadas / agrupadas / ordenadas por agujero y el mapa de map[][] densidad map[][] .

  1. recorrer todos los agujeros

    1. recorrer todas las líneas H, V del agujero procesado

      Cree una lista de todos los puntos finales de línea únicos pnt[] (sin duplicados). Así que tome 2 puntos finales para cada línea y observe si cada punto ya está en la lista. Si no lo agregas allí, ignóralo.

    2. eliminar todos los puntos no fronterizos de la lista

      Por lo tanto, elimine todos los puntos que no tienen contacto con el área sin orificios mirando en 4 vecinos en el map[][] densidad map[][]

    3. hacer análisis de componentes conectados en los puntos

      1. set used=0; p0=-1; p1=-1; used=0; p0=-1; p1=-1; para todos los puntos en la lista pnt[]
      2. conecta puntos con distance=1

        recorra todos los puntos pnt[] con used<2 que significa que aún no se utilizan por completo y para cada uno de estos puntos busque pnt[] nuevamente para otro punto que también tenga distance = 1 en él. Significa que son sus 4 vecinos y deben estar conectados, así que agregue la información de conexión a la cabina de ellos (use los índices p0 o p1 que no se usen nunca (-1) ) e incremente el uso de ambos puntos.

      3. intente conectar puntos con distance=sqrt(2)

        es casi lo mismo que # 2, excepto la distancia que ahora selecciona las diagonales de 8 vecinos. Esta vez también evite los bucles cerrados, así que no conecte el punto que ya está conectado a él.

      4. tratar de conectar los puntos más cercanos

        nuevamente es casi lo mismo que # 2, # 3, pero seleccione el punto más cercano y también evite los bucles cerrados.

      5. form polygon from pnt[]

        así que elija el primer punto en la lista y agréguelo al polígono. luego agregue el punto conectado a él (no importa de qué manera inicie p0 o p1 ). A continuación, agregue su punto conectado (diferente del punto agregado anterior al polígono para evitar bucles hacia atrás y hacia adelante). Agrega tantos puntos como puntos en un pnt[] .

La triangulación de Delauney puede ayudar. Tiene propiedad de que no hay ningún punto de entrada dentro del circunferencia circunscrita de ningún triángulo en la triangulación. Debido a eso, los puntos del límite del agujero estarán conectados por triangularjs más grandes / más amplios que cubren ese agujero. En sus casos, la triangulación tendrá muchos triangularjs de tamaño similar y algunos triangularjs de mayor tamaño que cubren los agujeros. Probablemente sea suficiente para filtrar los más grandes y conectarlos para encontrar un agujero.

Esta es mi entusiasta solución no científica:

1 – Escanee todo el área 2D con el paso mínimo predefinido (dx, dy). Para cada paso coord, encuentre el círculo más grande que podría caber sin ningún punto dentro. Deseche todos los círculos con un radio inferior a un tamaño predefinido.

enter image description here

2 – Ahora encuentre todos los grupos de círculos en colisión, prueba fácil de distancia y radio, almacene y agrupe en listas separadas. (Pregunte, si quiere más detalles sobre cómo agruparlos, es realmente fácil)

enter image description here

3 – Encuentra el polígono delimitador cóncavo para cada grupo de círculos, muy similar al algoritmo para encontrar el polígono convexo alrededor de un grupo de puntos que ya escribiste, y tus últimos angularjs de pregunta entre vectores estuvieron relacionados.

enter image description here

Notas

Sugerencias de optimización: antes del paso 1, puede almacenar todos los puntos en una cuadrícula para que el cálculo de la distancia se simplifique y limite a cuadrados cercanos a la cuadrícula del radio del círculo dado.

Precisión: Obtiene más precisión para valores más pequeños del paso de escaneo y el radio mínimo del círculo permitido.

No probado por mí mismo, pero estoy seguro de que funciona. ¡Buena suerte!

Probablemente sería mejor usar la triangulación de Delaunay para encontrar el gráfico de Gabriel . A continuación, ordena en ángulo el gráfico de Gabriel y realiza recorridos circulares para generar una lista de polígonos convexos. Luego puede ordenar esos polígonos por área. Estarías interesado en los que tienen el área más grande.

También será más eficiente modificar el gráfico ordenado por ángulo para que pueda seguir el camino de A a B, y ver qué sigue a continuación, ya sea en el sentido de las agujas del reloj o en el sentido contrario a las agujas del reloj (desde el ángulo). Un dictonario de diccionarios podría ser útil, que se define como algo así como “gráfico [A] [B] = (en el sentido de las agujas del reloj, en sentido antihorario)”. Un algoritmo de ejemplo (python) que utiliza el diccionario de enfoques de diccionarios.

 pre_graph = gabriel_graph(point) graph = {} for a in pre_graph: graph[a] = {} angle_sorted = sort(pre_graph[a], key=calc_angle_from(a)) for i,b in enumerate(angle_sorted): clockwise = angle_sorted[(i - 1) % len(angle_sorted)] counterclockwise = angle_sorted[(i + 1) % len(angle_sorted)] graph[a][b] = (clockwise, counterclockwise) polygons = [] for A in points: for B in graph[A]: for direction in [0,1]: polygon = [A] next_point = B: while next != A: polygon.append(next) next_point = graph[A][B][direction] if polygon[0] == min(polygon): # This should avoid duplicates polygons.add(polygon) 

También puede ser útil combinarlo con la sugerencia de Ante.

No sé de ningún algoritmo en la cabeza, pero una cosa que podría intentar (y este es mi primer impulso aquí) es algo similar a cómo se calcula la densidad en métodos libres de malla como la hidrodinámica de partículas suavizadas .

Si pudieras calcular un valor de densidad para cualquier punto en el espacio (y no solo en los puntos de muestra que tienes) podrías definir los límites de los agujeros como curvas de nivel de la función de densidad. Es decir, los agujeros son donde la densidad cae por debajo de un umbral (que probablemente le permita al usuario configurar). Puede encontrar los límites usando algo así como marchar cuadrados .

Si desea alguna aclaración sobre el funcionamiento de este tipo de funciones de interpolación de densidad, puedo proporcionarlo (a mi leal saber y entender) si lo desea.

No sé qué tan bien funcionaría realmente, pero espero que te dé una dirección.

Aquí hay un pensamiento:

  • For each point x find the distance d(x,y) (where y is the closest neighbor to x ). Define f(x)=d(x,y) as above.
  • Find the mean and variance of f(x) .
  • Find the ‘outliers’ – the points that their f values are very far from the mean, far by at least \alpha standard deviations. (\alpha is a parameter for the algorithm).
  • This will detect the ‘holes’ – all you have to do now is set the outlier of each hole.

It seems that you can address this by means of (binary) mathematical morphology on images.

Create a white image and plot all your points. Then “inflate” them into rectangles which are larger than the normal horizontal and vertical spacing. You do it by means of a so called erosion operation with a rectangular structuring element.

Doing this you will fill the plane, except at places where the points are too sparse.

The unfilled areas that you detect this way are smaller than the actual voids. You will restre to the full size by applying a dilation, using the same structuring element.

Both transforms combined are called an opening.

http://en.wikipedia.org/wiki/Mathematical_morphology